more fixes

Author: xtalax

Manifest.toml

@@ -31,6 +31,12 @@ git-tree-sha1 = "265b06e2b1f6a216e0e8f183d28e4d354eab3220"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 version = "3.2.1"
 
+[[ArrayLayouts]]
+deps = ["FillArrays", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "e1ba79094cae97b688fb42d31cbbfd63a69706e4"
+uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
+version = "0.7.8"
+
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
@@ -39,6 +45,12 @@ git-tree-sha1 = "45bb6705d93be619b81451bb2006b7ee5d4e4453"
 uuid = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
 version = "0.2.0"
 
+[[BandedMatrices]]
+deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "ce68f8c2162062733f9b4c9e3700d5efc4a8ec47"
+uuid = "aae01518-5342-5314-be14-df237901396f"
+version = "0.16.11"
+
 [[BangBang]]
 deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"]
 git-tree-sha1 = "0ad226aa72d8671f20d0316e03028f0ba1624307"
@@ -64,6 +76,18 @@ git-tree-sha1 = "bc1317f71de8dce26ea67fcdf7eccc0d0693b75b"
 uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
 version = "0.1.1"
 
+[[BlockArrays]]
+deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
+git-tree-sha1 = "5524e27323cf4c4505699c3fb008c3f772269945"
+uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+version = "0.16.9"
+
+[[BlockBandedMatrices]]
+deps = ["ArrayLayouts", "BandedMatrices", "BlockArrays", "FillArrays", "LinearAlgebra", "MatrixFactorizations", "SparseArrays", "Statistics"]
+git-tree-sha1 = "b1db5b5daca19070297580e1c5b5095e7ada6792"
+uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
+version = "0.11.1"
+
 [[CPUSummary]]
 deps = ["Hwloc", "IfElse", "Static"]
 git-tree-sha1 = "87b0c9c6ee0124d6c1f4ce8cb035dcaf9f90b803"
@@ -206,6 +230,12 @@ git-tree-sha1 = "0aa2d003ec9efe2a93f93ae722de05a870ffc0b2"
 uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
 version = "8.0.0"
 
+[[DiffEqOperators]]
+deps = ["BandedMatrices", "BlockBandedMatrices", "DiffEqBase", "DomainSets", "ForwardDiff", "LazyArrays", "LazyBandedMatrices", "LinearAlgebra", "LoopVectorization", "NNlib", "NonlinearSolve", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "SparseArrays", "SparseDiffTools", "StaticArrays", "SuiteSparse"]
+git-tree-sha1 = "8e7e2b9da87481294bcfcb7c7cf9e46832f843d7"
+uuid = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
+version = "4.38.0"
+
 [[DiffResults]]
 deps = ["StaticArrays"]
 git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
@@ -429,6 +459,18 @@ git-tree-sha1 = "83b56449c39342a47f3fcdb3bc782bd6d66e1d97"
 uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
 version = "0.1.4"
 
+[[LazyArrays]]
+deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "MatrixFactorizations", "SparseArrays", "StaticArrays"]
+git-tree-sha1 = "6dfb5dc9426e0cb7e237a71aa78c6b8c3e78a7fc"
+uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+version = "0.22.4"
+
+[[LazyBandedMatrices]]
+deps = ["ArrayLayouts", "BandedMatrices", "BlockArrays", "BlockBandedMatrices", "FillArrays", "LazyArrays", "LinearAlgebra", "MatrixFactorizations", "SparseArrays", "StaticArrays"]
+git-tree-sha1 = "8f4c427f6b505294dfefd343c8373b8c74e09973"
+uuid = "d7e5e226-e90b-4449-9968-0f923699bf6f"
+version = "0.7.6"
+
 [[LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
@@ -488,6 +530,12 @@ version = "0.1.6"
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
+[[MatrixFactorizations]]
+deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"]
+git-tree-sha1 = "1a0358d0283b84c3ccf9537843e3583c3b896c59"
+uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
+version = "0.8.5"
+
 [[MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
@@ -551,6 +599,12 @@ git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1"
 uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 version = "4.5.1"
 
+[[NNlib]]
+deps = ["Adapt", "ChainRulesCore", "Compat", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
+git-tree-sha1 = "afcfd88042cc5a62743c98ca311668ff5c03148a"
+uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+version = "0.7.32"
+
 [[NaNMath]]
 git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"

src/MethodOfLines.jl

@@ -4,6 +4,7 @@ module MethodOfLines
     using DiffEqBase
     using DiffEqOperators
     using ModelingToolkit
+    using ModelingToolkit: operation, istree, arguments, variable
     using SymbolicUtils, Symbolics
     import DomainSets
 

src/discretization/MOL_discretization.jl

@@ -12,10 +12,10 @@ end
 # Constructors. If no order is specified, both upwind and centered differences will be 2nd order
 function MOLFiniteDifference(dxs, time=nothing; upwind_order = 1, centered_order = 2, grid_align=center_align)
 
-    if discretization.centered_order % 2 != 0
-        throw(ArgumentError("Discretization centered_order must be even, given $(discretization.centered_order)"))
+    if centered_order % 2 != 0
+        warn("Discretization centered_order must be even, rounding up to $(centered_order+1)")
     end
-    return MOLFiniteDifference(dxs, time, approx_order, grid_align)
+    return MOLFiniteDifference(dxs, time, centered_order, grid_align)
 end
 
 function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::MethodOfLines.MOLFiniteDifference)
@@ -66,10 +66,10 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
             nottime = first(allnottime)
             @assert length(allindvars) == 1
             indvars = first(allindvars)
-            nspace = length(nottime)
-
-            s = DiscreteSpace(pdesys.domain, depvars, indvars, nottime, discretization)
-
+            
+            # Get the grid
+            s = DiscreteSpace(domain, depvars, indvars, nottime, discretization)
+            
             derivweights = DifferentialDiscretizer(pde, s, discretization)
 
             interior = s.Igrid[[2:(length(grid[x])-1) for x in s.nottime]]

src/discretization/discretize_vars.jl

@@ -73,8 +73,6 @@ function DiscreteSpace(domain, depvars, indvars, nottime, discretization)
     if grid_align == center_align
         grid = axies
     elseif grid_align == edge_align
-        # boundary conditions implementation assumes centered_order=2
-        @assert discretization.centered_order == 2
         # construct grid including ghost nodes beyond outer edges
         # e.g. space 0:dx:1 goes to grid -dx/2:dx:1+dx/2
         grid =  map(nottime) do x

src/discretization/generate_finite_difference_rules.jl

@@ -120,7 +120,7 @@ function spherical_diffusion(innerexpr, II, derivweights, s, r, u)
     if substitute(r, _rsubs(r, II)) ≈ 0
         D_2_u = central_difference(D_2, II, s, (s.x2i(r), r), u, ufunc_u)
         return 3exprhere*D_2_u # See appendix B of the paper
-
+    end
     D_1_u = central_difference(D_1, II, s, (s.x2i[r], r), u, ufunc_u)
     # See scheme 1 in appendix A of the paper
 
@@ -186,20 +186,22 @@ function generate_finite_difference_rules(II, s, pde, derivweights)
 
     rhs_arg = istree(pde.rhs) && (SymbolicUtils.operation(pde.rhs) == +) ? SymbolicUtils.arguments(pde.rhs) : [pde.rhs]
     lhs_arg = istree(pde.lhs) && (SymbolicUtils.operation(pde.lhs) == +) ? SymbolicUtils.arguments(pde.lhs) : [pde.lhs]
+    #These have to be done seperately to ensure that the spherical rules are applied first.
     nonlinlap_rules = []
     for t in vcat(lhs_arg,rhs_arg)
-        for r in cartesian_deriv_rules
+        for r in spherical_deriv_rules
             if r(t) !== nothing
                 push!(nonlinlap_rules, t => r(t))
             end
         end
-        for r in spherical_deriv_rules
+    end
+    for t in vcat(lhs_arg,rhs_arg)
+        for r in cartesian_deriv_rules
             if r(t) !== nothing
                 push!(nonlinlap_rules, t => r(t))
             end
         end
     end
-
     rules = vcat(vec(nonlinlap_rules),
                 vec(central_deriv_rules_cartesian),
                 vec(central_deriv_rules_spherical),