refactor: update to ModelingToolkit@10

Author: AayushSabharwal

docs/src/howitworks.md

@@ -12,4 +12,4 @@ Next, the boundary conditions are discretized, creating an equation for each poi
 
 After that, the system of PDEs is discretized, creating a finite difference equation for each point in their interior. Specific terms are recognized, and the best implemented scheme for these terms dispatched. For example, advection terms are discretized with the upwind scheme. There are also special schemes for the nonlinear Laplacian and spherical Laplacian. See [here for how this term matching occurs](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/generate_finite_difference_rules.jl), note that the order the generated rules are applied is important, with more specific rules applied first to avoid their terms being matched incorrectly by more general rules.  The [SymbolicUtils.jl docs](https://symbolicutils.juliasymbolics.org/rewrite/) are a useful companion here. See [here for the practical implementation of the finite difference schemes](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/differential_discretizer.jl).
 
-Now we have a system of equations which are either ODEs, linear, or nonlinear equations and an equal number of unknowns. See [here for the system](@ref brusssys) that is generated for the Brusselator at low point count. The structure of the system is simplified with `ModelingToolkit.structural_simplify`, and then either an `ODEProblem` or `NonlinearProblem` is returned. Under the hood, the `ODEProblem` generates a fast semidiscretization, written in Julia with `RuntimeGeneratedFunctions`. See [here for an example of the generated code](@ref brusscode) for the Brusselator system at low point count.
+Now we have a system of equations which are either ODEs, linear, or nonlinear equations and an equal number of unknowns. See [here for the system](@ref brusssys) that is generated for the Brusselator at low point count. The structure of the system is simplified with `ModelingToolkit.mtkcompile`, and then either an `ODEProblem` or `NonlinearProblem` is returned. Under the hood, the `ODEProblem` generates a fast semidiscretization, written in Julia with `RuntimeGeneratedFunctions`. See [here for an example of the generated code](@ref brusscode) for the Brusselator system at low point count.

src/MOL_discretization.jl

@@ -66,15 +66,15 @@ function get_discrete(pdesys, discretization)
         [Num(x) => s.grid[x] for x in s.x̄], [Num(u) => s.discvars[u] for u in s.ū]))
 end
 
-function ModelingToolkit.ODEFunctionExpr(
+function ODEFunctionExpr(
         pdesys::PDESystem, discretization::MethodOfLines.MOLFiniteDifference)
     sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization)
     try
         if tspan === nothing
             @assert true "Codegen for NonlinearSystems is not yet implemented."
         else
-            simpsys = structural_simplify(sys)
-            return ODEFunctionExpr(simpsys)
+            simpsys = mtkcompile(sys)
+            return ODEFunction(simpsys; expression = Val{true})
         end
     catch e
         println("The system of equations is:")
@@ -94,7 +94,7 @@ function SciMLBase.ODEFunction(
         if tspan === nothing
             @assert true "Codegen for NonlinearSystems is not yet implemented."
         else
-            simpsys = structural_simplify(sys)
+            simpsys = mtkcompile(sys)
             if analytic !== nothing
                 analytic = analytic isa Dict ? analytic : Dict(analytic)
                 s = getfield(sys, :metadata).discretespace

src/discretization/staggered_discretize.jl

@@ -3,7 +3,7 @@ function SciMLBase.discretize(pdesys::PDESystem,
         analytic = nothing, kwargs...) where {G <: StaggeredGrid}
     sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization)
     try
-        simpsys = structural_simplify(sys)
+        simpsys = mtkcompile(sys)
         if tspan === nothing
             add_metadata!(get_metadata(sys), sys)
             return prob = NonlinearProblem(simpsys, ones(length(simpsys.unknowns));

test/pde_systems/MOL_NonlinearProblem.jl

@@ -14,7 +14,7 @@ using DomainSets
         0 ~ b + d - 2 * c,
         0 ~ d - 1]
     @named ns = NonlinearSystem(eqs, [a, b, c, d], [])
-    sns = structural_simplify(ns)
+    sns = mtkcompile(ns)
     prob = NonlinearProblem(sns, zeros(length(get_unknowns(sns))), [])
 
     sol = NonlinearSolve.solve(prob, NewtonRaphson())