Tests passing - Release

Author: xtalax

Manifest.toml

@@ -455,7 +455,7 @@ uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 version = "7.1.1"
 
 [[LinearAlgebra]]
-deps = ["Libdl"]
+deps = ["Libdl", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[LogExpFunctions]]
@@ -571,6 +571,10 @@ git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 version = "1.10.8"
 
+[[OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
 [[OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
@@ -659,7 +663,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[Random]]
-deps = ["Serialization"]
+deps = ["SHA", "Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[RandomNumbers]]
@@ -972,6 +976,10 @@ git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
 uuid = "700de1a5-db45-46bc-99cf-38207098b444"
 version = "0.2.2"
 
+[[libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
 [[nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

Project.toml

@@ -8,6 +8,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"

src/bcs/generate_bc_eqs.jl

@@ -18,12 +18,12 @@ function generate_maps(bcs, t, s, tspan, grid_align)
         depvarbcmaps = reduce(vcat,[subvar(depvar, edgevals(s)) .=> edgevar for (depvar, edgevar) in zip(s.vars, edgevars(s))])
     end
     # depvarderivbcmaps will dictate what to replace the Differential terms with in the bcs
-    if length(s.params) == 1
+    if length(s.nottime) == 1
         # 1D system
-        bc_der_orders = filter(!iszero,sort(unique([count_differentials(bc.lhs, s.params[1]) for bc in bcs])))
+        bc_der_orders = filter(!iszero,sort(unique([count_differentials(bc.lhs, s.nottime[1]) for bc in bcs])))
         left_weights(d_order,j) = calculate_weights(d_order, s.grid[j][1], s.grid[j][1:1+d_order])
         right_weights(d_order,j) = calculate_weights(d_order, s.grid[j][end], s.grid[j][end-d_order:end])
-        # central_neighbor_idxs(II,j) = [II-CartesianIndex((1:length(s.params).==j)...),II,II+CartesianIndex((1:length(s.params).==j)...)]
+        # central_neighbor_idxs(II,j) = [II-CartesianIndex((1:length(s.nottime).==j)...),II,II+CartesianIndex((1:length(s.nottime).==j)...)]
         left_idxs(d_order) = CartesianIndices(s.grid[1])[1:1+d_order]
         right_idxs(d_order,j) = CartesianIndices(s.grid[j])[end-d_order:end]
         # Constructs symbolic spatially discretized terms of the form e.g. au₂ - bu₁
@@ -32,7 +32,7 @@ function generate_maps(bcs, t, s, tspan, grid_align)
         subderivar(depvar,d_order,x) = substitute.(((Differential(x)^d_order)(depvar),), edgevals(s))
         # Create map of symbolic Differential terms with symbolic spatially discretized terms
         depvarderivbcmaps = []
-        for x in s.params
+        for x in s.nottime
             rs = (subderivar(depvar,bc_der_orders[j],x) .=> derivars[i][j] for j in 1:length(bc_der_orders), (i,depvar) in enumerate(s.vars))
             for r in rs
                 push!(depvarderivbcmaps, r)
@@ -45,7 +45,7 @@ function generate_maps(bcs, t, s, tspan, grid_align)
                     for depvar in s.discvars]
             # replace u(t,0) with (u₁ + u₂) / 2, etc
             depvarbcmaps = reduce(vcat,[subvar(depvar, edgevals(s)) .=> bcvars[i]
-                                for (i, depvar) in enumerate(s.vars) for param in s.params])
+                                for (i, depvar) in enumerate(s.vars) for param in s.nottime])
 
         end
     else
@@ -62,7 +62,7 @@ function generate_maps(bcs, t, s, tspan, grid_align)
     end
 
     # All unique order of derivates in BCs
-    bc_der_orders = filter(!iszero,sort(unique([count_differentials(bc.lhs, s.params[1]) for bc in bcs])))
+    bc_der_orders = filter(!iszero,sort(unique([count_differentials(bc.lhs, s.nottime[1]) for bc in bcs])))
     # no. of different orders in BCs
     n = length(bc_der_orders)
     return BoundaryMap(depvarbcmaps, DerivativeMap{n}(depvarderivbcmaps, bc_der_orders), initmaps)
@@ -72,8 +72,9 @@ end
 """
 Mutates bceqs and u0 by finding relevant equations and discretizing them
 """
-function generate_u0_and_bceqs!!(u0, bceqs, bcs, t, s, depvar_ops, tspan, discretization)
+function generate_u0_and_bceqs!!(u0, bceqs, bcs, s, depvar_ops, tspan, discretization)
     grid_align = discretization.grid_align
+    t=discretization.time
     boundarymap = generate_maps(bcs, t, s, tspan, grid_align)
       # Generate initial conditions and bc equations
     for bc in bcs
@@ -94,19 +95,19 @@ function generate_u0_and_bceqs!!(u0, bceqs, bcs, t, s, depvar_ops, tspan, discre
                     # Replace Differential terms in the bc lhs with the symbolic spatially discretized terms
                     # TODO: Fix Neumann and Robin on higher dimension
                     # Update: Fixed for 1D systems
-                    derivativecount = count_differentials(bc.lhs, s.params[1])
+                    derivativecount = count_differentials(bc.lhs, s.nottime[1])
                     derivativeindex = findfirst(isequal(derivativecount), boundarymap.derivatives.orders)
                     k = initindex%2 == 0 ? 2 : 1
 
-                    if (length(s.params) == 1) && (derivativeindex !== nothing)
+                    if (length(s.nottime) == 1) && (derivativeindex !== nothing)
                         lhs = substitute(bc.lhs, boundarymap.derivatives.map[derivativeindex + count(boundarymap.derivatives)*Int(floor((initindex-1)/2))][k]) 
                     else
                         lhs = bc.lhs
                     end
 
                     # Replace symbol in the bc lhs with the spatial discretized term
                     lhs = substitute(lhs,boundarymap.vars[initindex])
-                    rhs = substitute.((bc.rhs,), edgemaps(t, s)[bcargs])
+                    rhs = substitute.((bc.rhs,), edgemaps(s)[bcargs])
                     lhs = lhs isa Vector ? lhs : [lhs] # handle 1D
                     push!(bceqs,lhs .~ rhs)
                 end

src/discretization/MOL_discretization.jl

@@ -11,7 +11,7 @@ struct MOLFiniteDifference{T,T2} <: DiffEqBase.AbstractDiscretization
 end
 
 # Constructors. If no order is specified, both upwind and centered differences will be 2nd order
-MOLFiniteDifference(dxs, time; upwind_order = 1, centered_order = 2, grid_align=center_align) =
+MOLFiniteDifference(dxs, time=nothing; upwind_order = 1, centered_order = 2, grid_align=center_align) =
     MOLFiniteDifference(dxs, time, upwind_order, centered_order, grid_align)
 
 function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::MethodOfLines.MOLFiniteDifference)
@@ -63,7 +63,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
             indvars = first(allindvars)
             nspace = length(nottime)
 
-            s = DiscreteSpace(pdesys.domain, depvars, nottime, grid_align, discretization)
+            s = DiscreteSpace(pdesys.domain, depvars, indvars, nottime, grid_align, discretization)
 
             #---- Count Boundary Equations --------------------
             # Count the number of boundary equations that lie at the spatial boundary on
@@ -97,7 +97,6 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
 
 
             # Calculate buffers
-            I1 = oneunit(first(s.Igrid))
             #TODO: Update this when allowing higher order derivatives to correctly deal with boundary upwinding
             interior = s.Igrid[[let bcs = get_bc_counts(i)
                                     (1 + first(bcs)):length(g)-last(bcs)
@@ -112,11 +111,11 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
             end
             # For all cartesian indices in the interior, generate finite difference rules
             eqs = vec(map(interior) do II
-                rules = generate_finite_difference_rules(II, I1, s, pde, discretization)
+                rules = generate_finite_difference_rules(II, s, pde, discretization)
                 substitute(pde.lhs,rules) ~ substitute(pde.rhs,rules)
             end)
 
-            generate_u0_and_bceqs!!(u0, bceqs, pdesys.bcs, t, s, depvar_ops, tspan, discretization)
+            generate_u0_and_bceqs!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, discretization)
             push!(alleqs,eqs)
             push!(alldepvarsdisc, reduce(vcat, s.discvars))
         end

src/discretization/discretize_vars.jl

@@ -6,7 +6,7 @@ end
 
 
 @inline function get_edgevals(s, i)
-    return [s.params[i] => first(s.axies[i]), s.params[i] => last(s.axies[i])]
+    return [s.nottime[i] => first(s.axies[i]), s.nottime[i] => last(s.axies[i])]
 end
 
 @inline function subvar(depvar, edge_vals)
@@ -15,7 +15,8 @@ end
 struct DiscreteSpace{N,M}
     vars
     discvars
-    params # Note that these aren't necessarily @params
+    nottime # Note that these aren't necessarily @parameters
+    params
     axies
     grid
     dxs
@@ -24,7 +25,7 @@ struct DiscreteSpace{N,M}
     Iedge
 end
 
-function DiscreteSpace(domain, depvars, nottime, grid_align, discretization)
+function DiscreteSpace(domain, depvars, indvars, nottime, grid_align, discretization)
     t = discretization.time
     nspace = length(nottime)
     # Discretize space
@@ -72,7 +73,7 @@ function DiscreteSpace(domain, depvars, nottime, grid_align, discretization)
     Iedge = reduce(vcat, [[vcat([Colon() for j = 1:i-1], 1, [Colon() for j = i+1:nspace]),
         vcat([Colon() for j = 1:i-1], length(axies[i]), [Colon() for j = i+1:nspace])] for i = 1:nspace])
 
-    return DiscreteSpace{nspace,length(depvars)}(depvars, depvarsdisc, nottime, axies, grid, dxs, Iaxies, Igrid, Iedge)
+    return DiscreteSpace{nspace,length(depvars)}(depvars, depvarsdisc, nottime, indvars, axies, grid, dxs, Iaxies, Igrid, Iedge)
 end
 
 nparams(::DiscreteSpace{N,M}) where {N,M} = N
@@ -82,19 +83,19 @@ grid_idxs(s::DiscreteSpace) = CartesianIndices(((axes(g)[1] for g in s.grid)...,
 edge_idxs(s::DiscreteSpace{N}) where {N} = reduce(vcat, [[vcat([Colon() for j = 1:i-1], 1, [Colon() for j = i+1:N]),
     vcat([Colon() for j = 1:i-1], length(s.axies[i]), [Colon() for j = i+1:N])] for i = 1:N])
 
-axiesvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.params[i] => s.axies[i][y.I[i]] for i = 1:N], s.Iaxies)
-gridvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.params[i] => s.grid[i][y.I[i]] for i = 1:N], s.Igrid)
+axiesvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.nottime[i] => s.axies[i][y.I[i]] for i = 1:N], s.Iaxies)
+gridvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.nottime[i] => s.grid[i][y.I[i]] for i = 1:N], s.Igrid)
 
 ## Boundary methods ##
-edgevals(s::DiscreteSpace{N}) where {N} = reduce(vcat, [get_edgevals(s.params, s.axies, i) for i = 1:N])
+edgevals(s::DiscreteSpace{N}) where {N} = reduce(vcat, [get_edgevals(s.nottime, s.axies, i) for i = 1:N])
 edgevars(s::DiscreteSpace) = [[d[e...] for e in s.Iedge] for d in s.discvars]
 
-@inline function edgemaps(t, s::DiscreteSpace)
-    bclocs(s::DiscreteSpace) = map(e -> substitute.(vcat(t, s.params), e), edgevals(s))
+@inline function edgemaps(s::DiscreteSpace)
+    bclocs(s::DiscreteSpace) = map(e -> substitute.(s.params, e), edgevals(s))
     return Dict(bclocs(s) .=> [axiesvals(s)[e...] for e in s.Iedge])
 end
 
-map_symbolic_to_discrete(II::CartesianIndex, s::DiscreteSpace{N,M}) where {N,M} = vcat([s.vars[k] => s.discvars[k][II] for k = 1:M], [s.params[j] => s.grid[j][II[j]] for j = 1:N])
+map_symbolic_to_discrete(II::CartesianIndex, s::DiscreteSpace{N,M}) where {N,M} = vcat([s.vars[k] => s.discvars[k][II] for k = 1:M], [s.nottime[j] => s.grid[j][II[j]] for j = 1:N])
 
 
 #varmap(s::DiscreteSpace{N,M}) where {N,M} = [s.vars[i] => i for i = 1:M]