Build ODEProblem from array differential equations with complete only

[ChrisRackauckas-Claude]

Summary

Rewritten (no longer stacked on #5100, which is closed): this PR makes ODEProblem(complete(sys), ...) work for systems whose differential equations are written over array slices, the same way DAEProblem already does, without mtkcompile and without any new keyword.

For a finite-difference discretization such as D(u[2:(n - 1)]) ~ lap(u) (or the residual form D(u[2:(n - 1)]) .- lap(u) ~ 0 that MethodOfLines emits), ODEProblem used to throw "The system has array equations. Call mtkcompile", and mtkcompile scalarizes, so the generated code and compile time grow with n. Now the array equation stays one equation and the generated function is independent of n.

Making mtkcompile (tearing / index reduction) keep array equations intact is separate work; this PR stops at complete → ODEProblem.

What changes (all in ModelingToolkitBase)

• generate_rhs (explicit ODE branch): array equations are rewritten to D(x) ~ f (explicit_array_derivative_form, which handles D(x) .- f ~ 0 and D(x) .+ g ~ 0), then the right-hand sides are packed into du with array_residual_maker — the same contiguous-region assembly the implicit-DAE path uses. check_lhs accepts D(x) where x is an array or slice whose elements are all unknowns; check_operator_variables checks uniqueness of the derivatives element-wise, so D(u[1:3]) together with D(u[3]) is still rejected.
• calculate_massmatrix gives an array equation one row per element (identity rows for differential equations, zero rows for zeros(n) ~ g), laid out like the generated code.
• __process_SciMLProblem skips the array-equation check for ODEFunction constructors (array unknowns are still rejected: pass collect(u) as unknowns, as for DAEProblem).
• full_equations expands array equations into scalar rows (after the residual → explicit rewrite), so the symbolic jacobian, sparsity pattern, tgrad and the initialization system see one equation per row. jac = true, sparse = true and the default initialization problem work from such systems.
• ODEFunction only computes W_sparsity when sparse or sparsity is set; it was computed unconditionally and is an O(n) symbolic computation that nothing consumed otherwise.
• build_explicit_observed_function records the element derivatives of a slice equation instead of calling default_toterm on D(u[2:4]) (which has no toterm name and threw).

Example

using ModelingToolkit, OrdinaryDiffEqRosenbrock
using ModelingToolkit: t_nounits as t, D_nounits as D

n = 96
@variables u(t)[1:n]
dx = 1 / (n - 1)
lap = (u[1:(n - 2)] .- 2 .* u[2:(n - 1)] .+ u[3:n]) ./ dx^2
@named heat = System([0 ~ u[1], D(u[2:(n - 1)]) ~ lap, 0 ~ u[n]], t, collect(u), [])
sys = complete(heat)

prob = ODEProblem(sys, [u => sinpi.(range(0, 1, length = n))], (0.0, 0.1))
prob.f.mass_matrix   # Diagonal([0, 1, ..., 1, 0])
sol = solve(prob, Rodas5P())

D(u) ~ -u over a whole array gives mass_matrix === I and solves with Tsit5.

Tests

lib/ModelingToolkitBase/test/array_equation_ode.jl (new; array_equation_dae.jl is also now wired into runtests.jl):

• D(u) ~ -u from complete: mass_matrix === I, in-place and out-of-place functions, solution with Tsit5, with and without the default initialization problem;
• heat equation over a slice in explicit and residual form: diagonal mass matrix, du equals the discrete Laplacian, solution matches the analytic solution, the mtkcompile path and a hand-scalarized complete path to 1e-10, default initialization works;
• the generated in-place and out-of-place Exprs have the same node count for n = 24, 48, 96 and contain no array_literal;
• a 2D slice D(w[2:n-1, 2:n-1]) ~ lap matches the scalarized system row for row;
• jac = true, sparse = true gives the tridiagonal-plus-boundary pattern and solves;
• full_equations returns n scalar rows for both forms;
• duplicated derivatives, slices with non-unknown elements and array unknowns are rejected.

Checked separately against the real (tearing) mtkcompile from ModelingToolkit for n = 21, 41: complete → ODEProblem → Rodas5P agrees with mtkcompile → ODEProblem → Tsit5 to 2e-11, with and without the initialization problem.

Notes

• Algebraic equations from complete still have to be written 0 ~ g (as before for scalar systems); u[1] ~ 0 keeps raising "The LHS cannot contain nondifferentiated variables". MethodOfLines emits its boundary conditions in the latter form, so its ODEProblem path needs either that rewrite on the MethodOfLines side or a follow-up here.
• Initialization from array equations is built from the scalarized rows, so it is correct but not O(1) in n; build_initializeprob = false gives the O(1) construction.

Made with Cursor

[ChrisRackauckas-Claude]

Checkpoint before usage cutoff

• This PR (Build ODEProblem from array differential equations with complete only #5101) remains the first useful slice: complete → ODEProblem for array DEs.
• Hard-path WIP saved as draft checkpoints (not the closed Preserve array differential equations with mtkcompile(sys; scalarize_arrays = false) #5100 skip-tearing approach):
  • CHECKPOINT: preserve array DEs through mtkcompile (ArrayEquationGroup / not skip-tearing) #5102
  • MTKTearing: keep array differential equations intact through structural simplification JuliaComputing/StateSelection.jl#158
• Full done/next writeup: https://github.com/ChrisRackauckas-Claude/ModelingToolkit.jl/blob/array-aware-tearing/docs/src/internals/ARRAY_EQUATIONS_CHECKPOINT.md

[ChrisRackauckas-Claude]

CI status after 232d7a4 (the array-equation DAE tests now solve from a consistent du0 with NoInit(); BrownFullBasicInit trips a ParameterDespecializationWrapper assertion in the downgrade environment's DiffEqBase 7.15, which aborted the group before the new ODE tests ran):

• Array Equation ODEProblem Test (69) and Array Equation DAEProblem Test (15) pass on Julia 1 / lts / pre and in the downgrade InterfaceII job.
• Root tests / InterfaceI, InterfaceII, Initialization, SymbolicIndexingInterface, Downgrade InterfaceI, and the MethodOfLines / Catalyst / MTKStdlib downstream jobs pass.
• Remaining red jobs are the same as on master (9b458b0): MTKBase InterfaceII fails only in JumpSystem Test (abstol field error), QA fails on the same 234 JET may be undefined reports from @match expansions, x86 InterfaceI can't precompile QuasiMonteCarlo, and StructuralIdentifiability hits a power_series_solution MethodError inside its own Nemo code.

[ChrisRackauckas-Claude]

Superseded: the ODEProblem array-DE support here was folded into #5148 ("Unify array-equation support behind a capability trait; fold in #5101 ODE support"), merged upstream and released in ModelingToolkitBase 1.75 / MTK 11.45. Verified on the released version: ODEProblem(complete(sys), ...) builds and solves a whole-array differential equation without mtkcompile. Closing to keep the queue honest.