Use MethodOfLines v1 array-form DAE path

[ChrisRackauckas]

Summary

• require MethodOfLines v1 in the test and documentation environments
• add DEIM reduction directly from a symbolic SciMLBase DAEProblem and its solution
• keep the full-order MethodOfLines model on the array-form DAE compilation path, then compile only the reduced ODE
• migrate the README, tutorial, and end-to-end FitzHugh-Nagumo test away from symbolic_discretize and the scalarized full-order ODE path

This implements the MethodOfLines v1 path referenced by SciMLDocs#350 and supersedes the dependency-only update in #174.

Implementation

MethodOfLines v1 now constructs the full problem with discretize(...; fallback = false) and solves that DAE directly. For POD-DEIM, the new overload:

1. symbolically tears the stored DAE system into its explicit differential subsystem without constructing a full-order ODEProblem
2. maps the torn differential unknowns back to the corresponding saved DAE state rows by symbolic identity
3. evaluates nonlinear snapshot terms by symbolic substitution, avoiding full-order scalar function generation
4. clears stale full-order initialization data and compiles only the reduced ODE
5. preserves and extends the observed equations so MethodOfLines PDE variable indexing still works on the reduced solution

Projecting all raw DAE residuals would change the existing explicit-ODE DEIM semantics, so this keeps the same algebraic elimination semantics as the old compiled path.

The ModelingToolkit and SciMLBase lower bounds match the MethodOfLines v1 requirements.

Tests

• GROUP=Core Pkg.test(): 71/71 passed
  • verifies fallback=false returns 4 array equations for 12 DAE states
  • solves the full DAE and reduced ODE successfully
  • reconstructs both PDE fields from reduced-solution metadata
  • retains coverage of the existing explicit ODESystem overload
• full Documenter build, doctests, rendering, and link checks: passed
• Runic check: passed
• typos: passed
• focused JET tests: 5/5 passed
• local package-wide JET reporting encountered the existing TSVD constant-redefinition issue; the remaining Aqua checks passed and current remote QA is green

References

• SciMLDocs#350
• feat!: make array-form DAE discretization the v1 default MethodOfLines.jl#650
• Support MethodOfLines v1 array-form reductions #174

[ChrisRackauckas]

Implementation scratchpad and commit map:

1. 1c1bf50 adds the symbolic DAE-to-DEIM bridge, SciMLBase dependency and v1 compatibility floors, array-DAE structural regression, end-to-end full/reduced solves, initialization cleanup, and legacy explicit-system coverage.
2. 3ef906c migrates the README and tutorial to discretize(...; fallback = false), solve(full_prob), and deim(full_prob, sol, ...), and restricts docs to MethodOfLines v1.

Validated after both commits:

• Core: 71/71 passed
• full Documenter build and doctests: passed
• Runic check: passed
• typos: passed
• diff check: clean

The full-order path never constructs an ODEProblem or generated scalar nonlinear function. Tearing is used only to recover the explicit differential equations and the saved DAE states are reordered by symbolic identity before reduction.

[ChrisRackauckas]

CI follow-up:

• The original Core / Julia 1 job remained in its test step beyond the repository's one-hour stuck threshold, so I cancelled that workflow attempt and reran only the cancelled job.
• The rerun passed in 28m45s, consistent with recent successful current-Julia runs (roughly 26-31 minutes).
• All 10 PR checks now pass, including current Julia, Julia LTS, Julia prerelease, QA, downgrade, documentation, Runic, and typos.

Local verification on 3ef906cbfde00685147341c4bcaebca2d94b7911 also passed: Core tests 71/71, full docs build and doctests, Runic check, typos, and git diff --check.

[ChrisRackauckas-Claude]

Adversarial review

Verdict

Do not merge this as the “fully array-form MOR” solution.

It successfully moves the full-order MethodOfLines solve onto the array-form DAEProblem path, but the MOR transformation and returned ROM remain substantially scalarized. PR 174 is much closer to the correct architecture.

Findings

1. Blocker: the reduced model is not array-form.

   This PR scalarizes the reduced variables and differential equations, then adds one scalar reconstruction equation per full-order dynamic state. Consequently, the observed/reconstruction graph grows at least as O(n·k) with grid size. It only preserves array form for the full-order solve, not the generated ROM. See

   ModelOrderReduction.jl/src/deim.jl

   Lines 236 to 281 in 3ef906c

   	var_name = gensym(:ŷ)
   	ŷ = (@variables $var_name(iv)[1:pod_dim])[1]
   	@set! sys.unknowns = Symbolics.value.(Symbolics.scalarize(ŷ)) # new variables from POD
   	ModelingToolkit.get_var_to_name(sys)[SymbolicIndexingInterface.getname(ŷ)] = Symbolics.unwrap(ŷ)
   	deqs, eqs = get_deqs(sys) # split eqs into differential and non-differential equations
   	rhs = [eq.rhs for eq in deqs]
   	# a sparse matrix of coefficients for the linear part,
   	# a vector of constant terms and a vector of nonlinear terms about dvs
   	A, g, F = separate_terms(rhs, dvs, iv)
   	nonlinear_snapshot = if isnothing(snapshot_times)
   	# Generate an in-place function from the symbolic nonlinear expressions for the
   	# existing explicit-system API.
   	F_func! = build_function(F, dvs; expression = Val{false}, kwargs...)[2]
   	values = similar(snapshot)
   	for i in axes(snapshot, 2)
   	F_func!(view(values, :, i), view(snapshot, :, i))
   	end
   	values
   	else
   	# The DAE entry point deliberately avoids generating a full-order scalar function.
   	_evaluate_symbolic_snapshot(F, dvs, snapshot, iv, snapshot_times)
   	end
   	deim_reducer = POD(nonlinear_snapshot, deim_dim)
   	reduce!(deim_reducer, TSVD())
   	U = deim_reducer.rbasis # DEIM projection basis
   	reduced_rhss, linear_projection_eqs = deim(dvs, A, g, F, ŷ, V, U; kwargs...)
   	reduced_deqs = D.(ŷ) ~ reduced_rhss
   	@set! sys.eqs = [Symbolics.scalarize(reduced_deqs); eqs]
   	old_observed = ModelingToolkit.get_observed(sys)
   	new_observed = [old_observed; linear_projection_eqs]
   	@set! sys.observed = _sort_observed_equations(new_observed)
   	# Replace full-order initialization data with the projected initial state. Array-form
   	# systems can retain parent-array guesses that are not keyed by the scalarized `dvs`.
   	@set! sys.guesses = Dict{Any, Any}()
   	@set! sys.initialization_eqs = Equation[]
   	reduced_initial = V' * view(snapshot, :, 1)
   	@set! sys.initial_conditions = Dict(Symbolics.unwrap(ŷ) => reduced_initial)
   	return complete(sys)

   .

   PR 174 instead emits one array differential equation and one array reconstruction observable per source field, storing the n×k lift as an array parameter. That is the appropriate online representation: offline analysis may scale with n, but the generated symbolic graph should not. See https://github.com/SciML/ModelOrderReduction.jl/blob/2e95959e2d1c15d7f142c59267cd47006cb3aa64/src/deim.jl#L322-L405.

2. Blocker: parameterized models are broken.

   The new snapshot evaluator substitutes states and time but never parameter values. I reproduced the complete DAE path with a parameterized reaction term; it fails with:

   ArgumentError: nonlinear snapshot expression ((u(t))[10]^2)*α
   evaluated to 0.0009549150281252625α after substituting states and time
   

   The cause is visible at

   ModelOrderReduction.jl/src/utils.jl

   Lines 28 to 57 in 3ef906c

   	function _evaluate_symbolic_snapshot(
   	expressions::AbstractVector, variables::AbstractVector,
   	snapshot::AbstractMatrix, iv, times::AbstractVector
   	)
   	size(snapshot, 1) == length(variables) ||
   	throw(DimensionMismatch("snapshot rows must match the differential unknowns"))
   	size(snapshot, 2) == length(times) ||
   	throw(DimensionMismatch("snapshot columns must match the saved times"))
   	values = similar(snapshot, length(expressions), size(snapshot, 2))
   	substitutions = Dict{Any, Any}()
   	for column in axes(snapshot, 2)
   	empty!(substitutions)
   	for (variable, value) in zip(variables, view(snapshot, :, column))
   	substitutions[variable] = value
   	end
   	substitutions[iv] = times[column]
   	for row in eachindex(expressions)
   	value = Symbolics.value(substitute(expressions[row], substitutions))
   	value isa Number ||
   	throw(
   	ArgumentError(
   	"nonlinear snapshot expression $(expressions[row]) evaluated to $(value) " *
   	"after substituting states and time"
   	)
   	)
   	values[row, column] = value
   	end
   	end
   	return values

   .

   Separately, the shared implementation replaces the entire defaults/initialization dictionary with only the reduced initial state. I reproduced a system retaining α in get_ps while losing its numeric default. This regresses the pre-existing deim(::ODESystem, ...) path too. PR 174 explicitly preserves and substitutes parameter defaults.

3. High: the tests do not establish numerical correctness.

   The end-to-end test checks only successful retcodes, dimensions, and reconstruction shapes. It never compares the ROM against the full DAE or the old scalarized reduction. Incorrect row ordering, forcing, algebraic elimination, or reconstruction values could all pass. See

   ModelOrderReduction.jl/test/deim.jl

   Lines 36 to 65 in 3ef906c

   	N = 5 # (minimum number of) equidistant discretization intervals
   	dx = (L - 0.0) / N
   	dxs = [x => dx]
   	order = 2
   	discretization = MOLFiniteDifference(dxs, t; approx_order = order)
   	dae_prob = discretize(pde_sys, discretization; fallback = false)
   	@test dae_prob isa DAEProblem
   	@test length(ModelingToolkit.get_eqs(dae_prob.f.sys)) == 4
   	@test length(ModelingToolkit.get_unknowns(dae_prob.f.sys)) == 12
   	sol = solve(dae_prob; saveat = 1.0)
   	@test successful_retcode(sol)
   	pod_dim = 3
   	deim_sys = @test_nowarn deim(dae_prob, sol, pod_dim)
   	# check the number of dependent variables in the new system
   	@test length(ModelingToolkit.get_unknowns(deim_sys)) == pod_dim
   	@test isempty(ModelingToolkit.initialization_equations(deim_sys))
   	deim_prob = ODEProblem(complete(deim_sys), nothing, dae_prob.tspan)
   	deim_sol = solve(deim_prob, Rodas5P(), saveat = 1.0)
   	@test successful_retcode(deim_sol)
   	nₓ = length(sol[x])
   	nₜ = length(sol[t])
   	# Test solution retrieval through the MethodOfLines metadata.
   	@test size(deim_sol[v(x, t)]) == (nₓ, nₜ)
   	@test size(deim_sol[w(x, t)]) == (nₓ, nₜ)

   .

   Required coverage should include deterministic RHS equivalence against the old path and a trajectory/reconstruction error assertion. PR 174 already adds parameter, multidimensional ordering, numerical-error, and grid-independent tree-size tests.

4. High: this is not general DAE reduction.

   The method dispatches on any symbolic DAEProblem, but immediately tears it into an explicit ODE and rejects any remaining algebraic residual. It is therefore ODE POD-DEIM after algebraic elimination, suitable only for DAEs reducible to one explicit differential equation per dynamic unknown. A genuine descriptor/PDAE requires residual or mass-matrix projection and should be a distinct algorithm/API. See

   ModelOrderReduction.jl/src/deim.jl

   Lines 333 to 375 in 3ef906c

   	function deim(
   	prob::SciMLBase.DAEProblem, sol::SciMLBase.AbstractODESolution,
   	pod_dim::Integer; deim_dim::Integer = pod_dim,
   	name::Union{Nothing, Symbol} = nothing, kwargs...
   	)::ODESystem
   	hasproperty(prob.f, :sys) && prob.f.sys isa ODESystem ||
   	throw(ArgumentError("the DAE problem must contain a symbolic ModelingToolkit system"))
   	raw_sys = prob.f.sys
   	hasproperty(sol, :prob) && hasproperty(sol.prob.f, :sys) && sol.prob.f.sys === raw_sys ||
   	throw(ArgumentError("the solution must have been obtained from the supplied DAE problem"))
   	first(sol.t) == first(prob.tspan) ||
   	throw(ArgumentError("the solution must save the state at the start of the DAE problem"))
   	raw_variables = ModelingToolkit.get_unknowns(raw_sys)
   	full_snapshot = Array(sol)
   	size(full_snapshot, 1) == length(raw_variables) ||
   	throw(DimensionMismatch("solution states must match the DAE system unknowns"))
   	sys = complete(tearing(raw_sys))
   	deqs, residual_equations = get_deqs(sys)
   	isempty(residual_equations) ||
   	throw(ArgumentError("tearing the DAE must produce an explicit ODE system"))
   	length(deqs) == length(ModelingToolkit.get_unknowns(sys)) ||
   	throw(ArgumentError("the torn DAE system must have one differential equation per unknown"))
   	raw_indices = Dict(
   	Symbolics.unwrap(variable) => index
   	for (index, variable) in enumerate(raw_variables)
   	)
   	snapshot_rows = map(ModelingToolkit.get_unknowns(sys)) do variable
   	index = get(raw_indices, Symbolics.unwrap(variable), nothing)
   	isnothing(index) &&
   	throw(
   	ArgumentError(
   	"a differential unknown produced by tearing is absent from the DAE state vector"
   	)
   	)
   	index
   	end
   	snapshot = full_snapshot[snapshot_rows, :]
   	reduced_name = isnothing(name) ? Symbol(nameof(raw_sys), :_deim) : name
   	return _deim_impl(
   	sys, snapshot, pod_dim;
   	deim_dim, name = reduced_name, snapshot_times = sol.t, kwargs...
   	)

   .

5. Medium: fallback=false does not guarantee array equations.

   It prevents MethodOfLines from falling back from DAEProblem to the compiled ODEProblem path. Unsupported discretization patterns can still fall back individually to pointwise scalar equations inside the DAE. MethodOfLines does this explicitly at https://github.com/SciML/MethodOfLines.jl/blob/v1.1.1/src/discretization/discretize_equations.jl#L121-L143. The single N=5, four-equation assertion does not establish grid independence generally.

O(1) compile-time result

For supported array stencils, yes: the full-order DAE’s symbolic/JIT compilation is consistent with O(1) in grid size. Total DAEProblem construction is not.

Using Julia 1.12.7, MethodOfLines 1.1.1, ModelingToolkit 11.40.0, and the same kind of two-field 1D stencil:

Unknowns	Equations	Fresh-process compiler time
12	4	105.56 s
102	4	105.36 s
1,002	4	107.48 s

After warming compilation, increasing from 102 to 1,002 to 10,002 unknowns gave compiler times of 0.50, 0.29, and 0.35 seconds. But total construction time grew from 0.50 to 0.54 to 20.44 seconds, and allocations from 37 MiB to 136 MiB to 7.99 GiB. MethodOfLines necessarily enumerates every scalar unknown when constructing its operating point: https://github.com/SciML/MethodOfLines.jl/blob/v1.1.1/src/dae_discretization.jl#L158-L243.

Therefore:

• Full-order array residual compilation: effectively grid-independent for supported patterns.
• Full DAEProblem construction: not O(1).
• This PR’s offline MOR transformation: not O(1) and heavily scalar/symbolic.
• This PR’s returned reconstruction graph: grid-dependent.
• PR 174’s array-valued online representation: directionally the right approach.

All ten remote checks are green. My local GROUP=Core rerun was stopped after an hour without reaching assertions because Julia 1.12.7 repeatedly rebuilt the temporary test environment and hit an internal irinterp is unable to handle heavy recursion correctly precompile error; this is not evidence of a PR test failure.

The PR also does not currently follow the repository workflow supplied for this review: it is not a draft, lacks the ignore-until-reviewed note, and its two commit messages have empty bodies without the required co-author/agent trailers.

Links

• PR 177: Use MethodOfLines v1 array-form DAE path #177
• Alternative PR 174: Support MethodOfLines v1 array-form reductions #174
• MethodOfLines v1 array DAE change: feat!: make array-form DAE discretization the v1 default MethodOfLines.jl#650
• PR 177 Core CI run: https://github.com/SciML/ModelOrderReduction.jl/actions/runs/32855718121

---

🤖 Review generated by an AI agent using Codex (harness version unavailable; model: GPT-5). Local session; no shareable conversation URL was exposed.