Merge pull request #414 from SciML/straydiffs

Handle D(0.001) type of post rule subsitutions

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "MethodOfLines"
 uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 authors = ["Alex Jones, <[email protected]>"]
-version = "0.11.3"
+version = "0.11.4"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/MethodOfLines.jl

@@ -7,7 +7,7 @@ using ModelingToolkit: operation, iscall, arguments, variable, get_metadata, get
 parameters, defaults, varmap_to_vars
 using SymbolicIndexingInterface
 using SymbolicUtils, Symbolics
-using Symbolics: unwrap, solve_for, expand_derivatives, diff2term, setname, rename
+using Symbolics: unwrap, symbolic_linear_solve, expand_derivatives, diff2term, setname, rename
 using SymbolicUtils: operation, arguments
 using IfElse
 using StaticArrays

src/discretization/generate_ic_defaults.jl

@@ -11,7 +11,7 @@ function PDEBase.generate_ic_defaults(tconds, s::DiscreteSpace,
             indexmap = Dict([args[i] => i for i in 1:length(args)])
             D = ic.order == 0 ? identity : (Differential(t)^ic.order)
             defaultvars = D.(s.discvars[depvar(ic.u, s)])
-            broadcastable_rhs = [solve_for(ic.eq, D(ic.u))]
+            broadcastable_rhs = [symbolic_linear_solve(ic.eq, D(ic.u))]
             out = substitute.(
                 broadcastable_rhs, valmaps(s, depvar(ic.u, s), ic.depvars, indexmap))
             vec(defaultvars .=>

src/scalar_discretization.jl

@@ -43,7 +43,7 @@ function discretize_equation_at_point(
         boundaryrules,
         valmaps(s, eqvar, depvars, II, indexmap))
     try
-        return substitute(pde.lhs, rules) ~ substitute(pde.rhs, rules)
+        return expand_derivatives(substitute(pde.lhs, rules)) ~ substitute(pde.rhs, rules)
     catch e
         println("A scheme has been incorrectly applied to the following equation: $pde.\n")
         println("The following rules were constructed at index $II:")

src/system_parsing/pde_system_transformation.jl

@@ -271,7 +271,7 @@ function generate_bc_rules(bcs, v)
         deriv = bc.order == 0 ? identity : (Differential(bc.x)^bc.order)
         bcrule_lhs = deriv(operation(bc.u)(v.args[operation(bc.u)]...))
         bcterm = deriv(bc.u)
-        rhs = solve_for(bc.eq, bcterm)
+        rhs = symbolic_linear_solve(bc.eq, bcterm)
         bcrule_lhs => rhs
     end
 end

test/pde_systems/MOL_Mixed_Deriv.jl

@@ -34,16 +34,15 @@ using ModelingToolkit: Differential
     dx = (xmax - xmin) / 20
     discretization = MOLFiniteDifference([x => dx], t, advection_scheme = WENOScheme())
 
-    @test_broken (discretize(pdesys, discretization) isa ODEProblem)
-    # prob = discretize(pdesys, discretization)
-    # sol = solve(prob, FBDF())
+    prob = discretize(pdesys, discretization)
+    sol = solve(prob, FBDF())
 
-    # xdisc = sol[x]
-    # tdisc = sol[t]
-    # usol = sol[u(t,x)]
+    xdisc = sol[x]
+    tdisc = sol[t]
+    usol = sol[u(t,x)]
 
-    # asol = [assf(t, x) for t in tdisc, x in xdisc]
-    # @test_broken usol ≈ asol atol = 1e-3
+    asol = [assf(t, x) for t in tdisc, x in xdisc]
+    @test_broken usol ≈ asol atol = 1e-3
 end
 
 @testset "Test 01: Dt(u) ~ Dxy(u)" begin

test/pde_systems/MOLtest2.jl

@@ -40,7 +40,7 @@ using Test
         diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
         @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
         discretization = MOLFiniteDifference([x => Δx], t)
-        @test_broken (discretize(pdesys, discretization) isa ODEProblem)
+        @test (discretize(pdesys, discretization) isa ODEProblem)
     end
 
     @testset "Test 02: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin
@@ -63,7 +63,7 @@ using Test
         diff_eq = ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))
         @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)])
         discretization = MOLFiniteDifference([x => Δx], t)
-        @test_broken (discretize(pdesys, discretization) isa ODEProblem)
+        @test (discretize(pdesys, discretization) isa ODEProblem)
     end
 
     @testset "Test 05: ∂t(c(x, t)) ~ ∂x(1/x * ∂x(c(x, t)))" begin