feat!: make array-form DAE discretization the v1 default

Project.toml

@@ -1,6 +1,6 @@
 name = "MethodOfLines"
 uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
-version = "0.12.0"
+version = "1.0.0"
 authors = ["Alex Jones, <alex.jones@juliahub.com>"]
 
 [deps]
@@ -11,6 +11,7 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PDEBase = "a7812802-0625-4b9e-961c-d332478797e5"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -23,12 +24,13 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Combinatorics = "1"
-DiffEqBase = "7.3"
+DiffEqBase = "7.18.1"
 DomainSets = "0.7, 0.8"
 ForwardDiff = "0.10, 1"
 Interpolations = "0.15, 0.16"
 Latexify = "0.16"
-ModelingToolkit = "11.31.2"
+ModelingToolkit = "11.39.1"
+ModelingToolkitBase = "1.67"
 NonlinearSolve = "4"
 OrdinaryDiffEq = "7"
 OrdinaryDiffEqLowOrderRK = "2"
@@ -39,7 +41,7 @@ PDEBase = "0.1.30"
 PrecompileTools = "1.2.1"
 RuntimeGeneratedFunctions = "0.5.12"
 SafeTestsets = "0.1"
-SciMLBase = "3.1"
+SciMLBase = "3.48"
 SciMLTesting = "2.8"
 StableRNGs = "1"
 StaticArrays = "1.9.18"

benchmark/weno/suite.jl

@@ -53,6 +53,11 @@ function solve_benchmark(prob, dt)
     ) evals = 1
 end
 
+function ode_problem(pdesys, disc)
+    sys, tspan = symbolic_discretize(pdesys, disc)
+    return ODEProblem(mtkcompile(sys), nothing, tspan)
+end
+
 function discretize_benchmark(pdesys, disc)
     return @benchmarkable discretize($pdesys, $disc) seconds = 60 samples = 3 evals = 1
 end
@@ -87,7 +92,7 @@ function build_weno_suite(;
             slv[label][kstr] = BenchmarkGroup()
             dsc[label][kstr] = BenchmarkGroup()
             for n in resolutions
-                prob = discretize(sys.pdesys, mkdisc(sys, kind, n))
+                prob = ode_problem(sys.pdesys, mkdisc(sys, kind, n))
                 rhs[label][kstr]["N=$n"] = rhs_benchmark(prob)
                 dt = CFL_TARGET * min_spacing(kind, sys.xspan..., n) / wavespeed
                 slv[label][kstr]["N=$n"] = solve_benchmark(prob, dt)
@@ -108,7 +113,7 @@ function build_weno_suite(;
         slv["interface"][kstr] = BenchmarkGroup()
         dsc["interface"][kstr] = BenchmarkGroup()
         for n in interface_resolutions
-            prob = discretize(itf.pdesys, interface_discretization(itf, kind, n))
+            prob = ode_problem(itf.pdesys, interface_discretization(itf, kind, n))
             rhs["interface"][kstr]["N=$n"] = rhs_benchmark(prob)
             dt = CFL_TARGET * interface_min_spacing(kind, n)
             slv["interface"][kstr]["N=$n"] = solve_benchmark(prob, dt)

docs/Project.toml

@@ -9,7 +9,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
 OrdinaryDiffEqSSPRK = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
 PDEBase = "a7812802-0625-4b9e-961c-d332478797e5"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -30,10 +29,9 @@ Documenter = "1"
 DomainSets = "0.7, 0.8"
 Interpolations = "0.14, 0.15, 0.16"
 Latexify = "0.15, 0.16"
-MethodOfLines = "0.11, 0.12"
+MethodOfLines = "0.11, 0.12, 1"
 ModelingToolkit = "11"
 OrdinaryDiffEq = "6, 7"
-OrdinaryDiffEqBDF = "2"
 OrdinaryDiffEqSSPRK = "2"
 PDEBase = "0.1.20"
 Plots = "1"

docs/src/MOLFiniteDifference.md

@@ -1,14 +1,16 @@
 # [Discretization](@id molfd)
+
 ```julia
-struct MOLFiniteDifference{G} <: DiffEqBase.AbstractDiscretization
-    dxs
-    time
+struct MOLFiniteDifference{G} <: AbstractEquationSystemDiscretization
+    dxs::Any
+    time::Any
     approx_order::Int
-    advection_scheme
+    advection_scheme::Any
     grid_align::G
     should_transform::Bool
-    use_ODAE::Bool
-    kwargs
+    useIR::Bool
+    callbacks::Any
+    kwargs::Any
 end
 ```
 
@@ -20,13 +22,12 @@ domain = [your domain, a vector of Intervals i.e. x ∈ Interval(x_min, x_max)]
 
 @named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
 
-discretization = MOLFiniteDifference(dxs, 
-                                      <your choice of continuous variable, usually time>; 
-                                      advection_scheme = <UpwindScheme() or WENOScheme()>, 
-                                      approx_order = <Order of derivative approximation, starting from 2> 
+discretization = MOLFiniteDifference(dxs,
+                                      <your choice of continuous variable, usually time>;
+                                      advection_scheme = <UpwindScheme() or WENOScheme()>,
+                                      approx_order = <Order of derivative approximation, starting from 2>,
                                       grid_align = <your grid type choice>,
-                                      should_transform = <Whether to automatically transform the PDESystem (see below)>
-                                      use_ODAE = <Whether to use ODAEProblem>)
+                                      should_transform = <Whether to automatically transform the PDESystem (see below)>)
 prob = discretize(pdesys, discretization)
 ```
 Where `dxs` is a vector of pairs of parameters to the grid step in this dimension, i.e. `[x=>0.2, y=>0.1]`. If the value given for a dimension is a subtype of `Integer`, the domain for that variable will be discretized in to that integer number of equally spaced points.
@@ -45,8 +46,65 @@ Currently supported options are `grid_align`: `center_align` and `edge_align`. E
 
 `should_transform`: Whether to automatically transform the system to make it compatible with MethodOfLines where possible, defaults to true. If your system has no mixed derivatives, all derivatives are purely of a dependent variable i.e. `Dx(u_aux(t,x))` not `Dx(v(t,x)*u(t,x))`, excepting nonlinear and spherical Laplacians for which this holds for the innermost derivative argument, and no expandable derivatives, this can be set to false for better discretization performance at the cost of generality, if you perform these transformations yourself.
 
-`use_ODAE`: MethodOfLines will automatically make use of `ODAEProblem` where relevant, which improves performance for DAEs (as discretized PDEs are in general), if this is set to true. Defaults to false.
+MethodOfLines generates the interior of each PDE as a single symbolic array equation over slices of the discretized variables, e.g. `D(u[2:n-1]) - (u[1:n-2] .- 2u[2:n-1] .+ u[3:n]) ./ dx^2 ~ 0`. This keeps the number of symbolic equations independent of the grid resolution and scales much better during symbolic processing. Nonlinear Laplacians `Dx(a(u) * Dx(u))`, spherical Laplacians `r^-2 Dr(r^2 Dr(u))`, mixed first-order derivatives `Dx(Dy(u))`, WENO / functional advection (on uniform grids, and on nonuniform grids — periodic or not — for schemes that provide a coefficient split; WENO does; user schemes opt in by defining a method on `MethodOfLines.array_scheme_split`), staggered grids, self-periodic interfaces and boundary values appearing inside an interior equation (e.g. `u(t, 1)`) are included in the array form, including on wrap boxes near a periodic seam. Patterns without a slice representation (nonuniform advection schemes without a coefficient split, schemes that read the grid coordinate, integrals, higher mixed derivatives, two-domain interface BCs, derivatives of boundary values, time-literal references such as `u(0, x)`, boundary values on edge-aligned grids, stationary systems) automatically fall back to pointwise scalar equations for the affected equation.
+
+Any unrecognized keyword arguments are passed to the generated problem constructor; see the [ModelingToolkit problem documentation](https://docs.sciml.ai/ModelingToolkit/stable/API/problems/#Dynamical-systems) for available options.
+
+## Problem types
+
+`discretize` returns a `DAEProblem` for a time-dependent system:
+
+```julia
+
+disc = MOLFiniteDifference([x => n], t)
+prob = discretize(pdesys, disc)
+sol = solve(prob)
+```
+
+MethodOfLines emits residuals of the form `D(u) - f ~ 0`, which are already in
+implicit-DAE form. Building a `DAEProblem` therefore needs no `mtkcompile`, and the array
+equations reach the generated code intact — isolating the derivative for an `ODEProblem`
+is structural simplification, and it scalarizes them. Calling `solve(prob)` lets
+OrdinaryDiffEq select its default DAE algorithm.
+
+`initializealg` defaults to `BrownFullBasicInit()`, chosen only when the discretized
+system's initialization equations are ones that algorithm preserves.
+
+A few systems cannot be posed as a first-order DAE: those second order in time, and those
+whose initialization equations `BrownFullBasicInit` would not honour. They fall back to
+`mtkcompile` plus an `ODEProblem`, which scalarizes the array equations. Pass
+`fallback = false` to `discretize` to make that an error instead.
+
+Time-independent systems have no derivative to keep implicit, and discretize to a
+`NonlinearProblem`.
+
+The solution is a `PDETimeSeriesSolution` in every case, indexed and interpolated by the
+`PDESystem`'s own variables: `sol[u(t, x)]`, `sol(t, x)`.
+
+## Building a problem yourself: `symbolic_discretize`
+
+To construct a problem type `discretize` does not build — an `ODEProblem`, or anything
+else — start from `symbolic_discretize`, which returns the discretized system and the
+time span:
+
+```julia
+sys, tspan = symbolic_discretize(pdesys, disc)
+
+# an ODEProblem needs `D(x) = f(x)`, so compile first
+prob = ODEProblem(mtkcompile(sys), nothing, tspan)
+sol = solve(prob, Rodas5P())
+```
+
+Note that `mtkcompile` scalarizes the array equations, so this path gives up the scaling
+benefit of the array form. Prefer `discretize` unless you specifically need an
+`ODEProblem`.
 
-`discretization_strategy`: How the discretized equations are represented symbolically. `ScalarizedDiscretization()` (the default) generates one scalar equation per interior grid point. `ArrayDiscretization()` generates the interior of each PDE as a single symbolic array equation over slices of the discretized variables, e.g. `D(u[2:n-1]) - (u[1:n-2] .- 2u[2:n-1] .+ u[3:n]) ./ dx^2 ~ 0`, which scales much better to large systems during symbolic processing. Nonlinear Laplacians `Dx(a(u) * Dx(u))`, spherical Laplacians `r^-2 Dr(r^2 Dr(u))`, mixed first-order derivatives `Dx(Dy(u))`, WENO / functional advection (on uniform grids, and on nonuniform grids — periodic or not — for schemes that provide a coefficient split; WENO does; user schemes opt in by defining a method on `MethodOfLines.array_scheme_split`), staggered grids, self-periodic interfaces and boundary values appearing inside an interior equation (e.g. `u(t, 1)`) are included in the array form, including on wrap boxes near a periodic seam. Patterns without a slice representation (nonuniform advection schemes without a coefficient split, schemes that read the grid coordinate, integrals, two-domain interface BCs, derivatives of boundary values, time-literal references such as `u(0, x)`, boundary values on edge-aligned grids, stationary systems) automatically fall back to pointwise scalar equations for the affected equation, matching `ScalarizedDiscretization` there. Where the array form is used, numerics match the scalar path whenever the scalar path can express the same boundary-value substitutions; the array path also substitutes periodic-face and free-standing-corner references that scalar `boundaryvalfuncs` currently leave symbolic.
+## Migrating to v1
 
-Any unrecognized keyword arguments will be passed to the `ODEProblem` constructor, see [its documentation](https://docs.sciml.ai/ModelingToolkit/stable/API/problems/#Dynamical-systems) for available options.
+- `discretize` returns a `DAEProblem` rather than an `ODEProblem` for time-dependent
+  systems. Call `solve(prob)` to use the default DAE algorithm; an `ODEProblem` solver
+  like `Tsit5()` will no longer accept the result. Solution indexing is unchanged.
+- Discretization strategy options were removed. MethodOfLines always uses array-form
+  equations with automatic pointwise fallback for unsupported patterns.
+- To construct the pre-v1 compiled `ODEProblem`, use `symbolic_discretize` plus
+  `mtkcompile` as above.

docs/src/faq.md

@@ -71,7 +71,8 @@ domains = [
 discretization = MOLFiniteDifference(
     [x => dx, y => dy], t; advection_scheme = UpwindScheme())
 
-prob = discretize(pdesys, discretization)
+sys, tspan = symbolic_discretize(pdesys, discretization)
+prob = ODEProblem(mtkcompile(sys), nothing, tspan)
 
 sol = solve(prob, SSPRK54(), dt = 0.01, saveat = 0.1)
 

docs/src/generated/bruss_ode_eqs.md

@@ -64,4 +64,4 @@ v[1, 5](t) ~ v[5, 5](t),
 v[1, 1](t) ~ 0] # Invalid corner point set to 0
 ```
 
-On the call to ODEProblem, this [code](@ref brusscode) is generated.
+On the compiled `ODEProblem` path, this [code](@ref brusscode) is generated.

docs/src/howitworks.md

@@ -12,4 +12,10 @@ 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.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.
+The result is a system of differential-algebraic, linear, or nonlinear equations with an
+equal number of unknowns. See [the generated Brusselator system](@ref brusssys) at low
+point count. Time-dependent systems are normally converted directly to a `DAEProblem`,
+which preserves symbolic array equations. Stationary systems are compiled with
+`ModelingToolkit.mtkcompile` into a `NonlinearProblem`; the optional compiled
+`ODEProblem` path also uses `mtkcompile` and scalarizes the equations. The problem
+constructors generate executable Julia functions with `RuntimeGeneratedFunctions`.

docs/src/solutions.md

@@ -1,8 +1,9 @@
 # [Solution Retrieval - PDESolutions](@id sol)
 
-MethodOfLines automatically wraps the `ODESolution` that results from its generated `ODEProblem` in
-a `PDESolution` object, reshaping it and providing a convenient interface for accessing your results,
-as well as interpolations.
+MethodOfLines wraps the solution of its generated `DAEProblem`, `ODEProblem`, or
+`NonlinearProblem` in a PDE solution object, reshaping it and providing a convenient
+interface for accessing results and interpolations. Time-dependent problems use
+`PDETimeSeriesSolution`; stationary problems use `PDENoTimeSolution`.
 
 ## Solution Retrieval
 
@@ -57,17 +58,18 @@ uv_interp = sol(0.4, 1.7, 2.6)
 
 ## Original solution
 
-The original `ODESolution` is stored in `sol.original_sol`.
+The original solver solution is stored in `sol.original_sol`.
 
 To avoid wrapping entirely, use the `wrap` keyword argument to `solve`:
 
 ```julia
-sol = solve(prob, Tsit5(); wrap = Val(false))
+
+sol = solve(prob; wrap = Val(false))
 ```
 
 ```
 > typeof(sol)
-ODESolution
+DAESolution
 ```
 
 This is useful where speed is important, but the shape of the solution is not.

docs/src/staggered.md

@@ -10,7 +10,7 @@ Solvers should be chosen carefully, the only officially supported solver is `Spl
 
 Staggered grid functionality is still in its infancy. Please open issues if unexpected results occur or needed functionality is not present.
 
-Staggered grids are compatible with `discretization_strategy=ArrayDiscretization()`, which emits the interior of each PDE as a single symbolic array equation over slices, keeping the symbolic equation count independent of the grid resolution.
+Staggered grids emit the interior of each PDE as a single symbolic array equation over slices, keeping the symbolic equation count independent of the grid resolution.
 
 ```
 using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
@@ -39,7 +39,8 @@ domains = [t in Interval(0.0, tmax),
 @named pdesys = PDESystem(eq, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]);
 
 discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid(), edge_aligned_var=ϕ(t,x));
-prob = discretize(pdesys, discretization);
+sys, tspan = symbolic_discretize(pdesys, discretization)
+prob = ODEProblem(mtkcompile(sys), nothing, tspan);
 
 sol = solve(prob, SplitEuler(), dt=dt);
 ```

docs/src/tutorials/PIDE.md

@@ -26,7 +26,7 @@ With BCs and ICs:
 We can discretize such a system like this:
 
 ```@example pide
-using MethodOfLines, ModelingToolkit, OrdinaryDiffEq, OrdinaryDiffEqBDF, DomainSets, Plots
+using MethodOfLines, ModelingToolkit, OrdinaryDiffEq, DomainSets, Plots
 
 @parameters t, x
 @variables u(..) cumuSum(..)
@@ -52,7 +52,7 @@ order = 2
 discretization = MOLFiniteDifference([x => 30], t)
 
 prob = MethodOfLines.discretize(pde_system, discretization)
-sol = solve(prob, QNDF(), saveat = 0.1);
+sol = solve(prob; saveat = 0.1);
 
 solu = sol[u(t, x)]
 
@@ -91,7 +91,7 @@ disc = MOLFiniteDifference([x => 120], t)
 
 prob = discretize(pde_system, disc)
 
-sol = solve(prob, Tsit5())
+sol = solve(prob)
 
 xdisc = sol[x]
 tdisc = sol[t]