Keep array equations through the compiled ODEProblem path

[ChrisRackauckas-Claude]

Summary

discretize preserves the array (slice-form) equations by building a DAEProblem without mtkcompile, but the ODEProblem path needed by explicit Runge–Kutta methods (Tsit5, SSPRK54, ...) still went through the scalarizing mtkcompile, so its compile time grew with the grid.

ModelingToolkit can now keep array equations through compilation with mtkcompile(sys; scalarize_arrays = false) (SciML/ModelingToolkit.jl#5100 and SciML/ModelingToolkit.jl#5101). This PR uses that for the compiled ODEProblem path, so the same O(1) examples work for ODEProblem.

What changes

• ode_compile(sys) (exported): mtkcompile(sys; scalarize_arrays = false) when the installed ModelingToolkit supports it and the result is an explicit ODE (every remaining equation is D(x) ~ f); otherwise the default mtkcompile. The fallback covers older ModelingToolkit versions (detected with isdefined(ModelingToolkitBase, :arrays_scalarized)) and systems the tearing-free compiler cannot make explicit (second order in time, coupled algebraic equations).
• _ode_problem, ODEFunction and ODEFunctionExpr go through ode_compile.
• New ODEProblem(pdesys, discretization; kwargs...) constructor, so the compiled path no longer requires symbolic_discretize + mtkcompile by hand.
• Docs (MOLFiniteDifference.md, index.md, howitworks.md, FAQ, tutorials) no longer say the ODE path always scalarizes, and use ode_compile / ODEProblem(pdesys, disc). NEWS entry for 1.5; version bumped to 1.5.0.

Tests

test/Discretization/problem_construction.jl:

• ODEProblem(pdesys, disc) for the 1D heat equation solves with Tsit5 to the analytic solution;
• "array-preserving ODE path": for 11, 21, 41 and 81 grid points the compiled system has n - 2 unknowns, mass_matrix === I, second-order accuracy, and — when ModelingToolkit preserves array equations — exactly one equation with a symbolic tree size independent of n. With an older ModelingToolkit the test instead asserts the scalarized shape, so CI passes with the current release too;
• the wave equation (second order in time) falls back to the scalarizing compilation and still solves.

Run locally against the ModelingToolkit branches above: problem_construction.jl 81/81 pass (with the released ModelingToolkit the fallback branch of the tests is exercised).

How to verify the O(1) scaling

using MethodOfLines, ModelingToolkit, DomainSets, OrdinaryDiffEq
@parameters t x
@variables u(..)
Dt = Differential(t); Dxx = Differential(x)^2
@named pdesys = PDESystem(Dt(u(t, x)) ~ Dxx(u(t, x)),
    [u(0, x) ~ sinpi(x), u(t, 0) ~ 0.0, u(t, 1) ~ 0.0],
    [t ∈ Interval(0.0, 0.1), x ∈ Interval(0.0, 1.0)], [t, x], [u(t, x)])
for n in (11, 81, 641)
    sys, _ = symbolic_discretize(pdesys, MOLFiniteDifference([x => n], t))
    csys = ode_compile(sys)
    @show n length(equations(csys))   # 1 for every n with MTK ≥ 11.43
    solve(ODEProblem(pdesys, MOLFiniteDifference([x => n], t)), Tsit5())
end

Depends on SciML/ModelingToolkit.jl#5100 and SciML/ModelingToolkit.jl#5101 for the preserved path; degrades to the previous behaviour without them.

Made with Cursor

[ChrisRackauckas-Claude]

Closing: this PR is built on mtkcompile(sys; scalarize_arrays = false) from SciML/ModelingToolkit.jl#5100, and that approach has been withdrawn (see the closing comment there). It compiled the system with tearing and index reduction switched off rather than making those algorithms keep array equations, and we do not want MethodOfLines (or its docs) to depend on that.

What replaces it, in two steps:

1. Build ODEProblem from array differential equations with complete only ModelingToolkit.jl#5101 (rewritten from master) makes ODEProblem(complete(sys), ...) accept array differential equations directly, like DAEProblem already does, with generated code independent of the grid size. The residual form MethodOfLines emits for the interior (-1 .* lap .+ D(u[2:n-1]) ~ 0 .+ 0 .* u[2:n-1]) is rewritten to D(u[2:n-1]) ~ ... automatically. Checked on the 1D heat equation from symbolic_discretize for 11, 41 and 81 points: complete → ODEProblem → Rodas5P matches mtkcompile → ODEProblem to 1e-9, with the same generated function for 41 and 81 points. The one thing standing between that and a MethodOfLines ODEProblem path is the form of the boundary conditions: complete → ODEProblem requires algebraic equations as 0 ~ g (as it always has for scalar systems), while MethodOfLines emits u[1] ~ 0.0. Either MethodOfLines writes its algebraic residuals as 0 ~ lhs - rhs when building this problem, or ModelingToolkit accepts the other form; that will be a small follow-up PR on whichever side, giving an O(1) mass-matrix ODEProblem (usable with Rodas/FBDF-type solvers; explicit RK solvers still need the boundary unknowns eliminated).

2. Getting the same O(1) behaviour through mtkcompile itself (so that the boundary conditions become observed and Tsit5 works) requires tearing and index reduction to preserve array equations. That is separate ModelingToolkit work and is not something MethodOfLines should work around in the meantime.

The DAEProblem path already on master is unaffected.

[ChrisRackauckas-Claude]

Checkpoint note: this PR still depends on the rejected scalarize_arrays=false / skip-tearing intermediate (MTK #5100, closed).

Preferred stack going forward:

1. MTK #5101 — complete → ODEProblem for array DEs
2. MTK #5102 + StateSelection Mixed Space and Time Derivatives #158 — draft checkpoint for preserving array DEs through real structural simplification (preserve_array_equations / ArrayEquationGroup), not turning tearing off

Please treat this MOL PR as needing a rewrite once those land (or close and reopen). Status doc: https://github.com/ChrisRackauckas-Claude/ModelingToolkit.jl/blob/array-aware-tearing/docs/src/internals/ARRAY_EQUATIONS_CHECKPOINT.md