[ChrisRackauckas-Claude]

Ignore this PR until reviewed by @ChrisRackauckas.

Follow-up to #178. The package is now organized around the symbolic array form of the reduced model: the only symbolic outputs are one array differential equation for the reduced state and one array (or scalar) reconstruction equation per source field. Everything that was left over from the scalarized representation is gone.

What changed

• The reduced equation is array linear algebra. It now reads D(ŷ) ~ A*ŷ + g + G*forcing + C*F where A = V'AV, g = V'g_numeric, G, C = V'U(P'U)⁻¹, and the stencil rows V_S of the state basis are array parameters of the reduced system, forcing is the array of distinct symbolic forcing expressions, and F is the array of DEIM-sampled nonlinear terms written in (V_S*ŷ)[i]. Before, V'*g and projector*sampled were evaluated eagerly into k scalar sums, so only the samples themselves (which are inherently scalar) remain scalar. Numeric matrices are parameters rather than literals because a sum of array terms holding literal matrices breaks mtkcompile (isless(::Matrix, ::Matrix) in SymbolicUtils argument sorting).
• Removed the scalar low-level deim(full_vars, A, g, F, ŷ, V, U) method. It returned scalarized projection equations y_i ~ (Vŷ)_i, which the high-level path discarded. The numeric projection now lives in the internal Projection struct and the symbolic assembly in _reduced_rhs; the initial derivative is evaluated numerically from the same projection.
• New src/full_order.jl. FullOrderModel holds the explicit first-order form of a source system (compiled as needed, linear * y + forcing + nonlinear split, source fields, snapshot row map, trained parameter values, preserved observed equations). separate_terms is reimplemented on top of the public Symbolics.semilinear_form instead of reaching into SymbolicUtils .dict/.coeff internals; it no longer takes the independent variable. The src/utils.jl helpers (get_deqs, only_dvs) and _sort_observed_equations are removed; the reduced observed equations are ordered by construction.
• deim(prob, sol, pod_dim) now accepts ODEProblems (Union{AbstractODEProblem, AbstractDAEProblem}) in addition to DAEProblems. ODEProblem(mtkcompile(rom), ...; build_initializeprob = false) is tested as a supported path alongside the array DAEProblem path.
• Dispatch and docstrings use ModelingToolkit.System; ODESystem is a deprecated alias in MTK v11. The Setfield dependency and the @set! sys.complete = false hack are gone (mtkcompile only rejects already-compiled systems). The unused AbstractMORProblem type is removed from the code and the developer docs. Parameter defaults are stored as plain numbers.
• Developer docs gain a "Reduction pipeline" section documenting FullOrderModel, SourceField, and separate_terms. AGENTS.md records the term-classification rule, the scalar-unknown requirement, and the two supported problem constructors.

Not changed: the numerical algorithm (POD via TSVD, DEIM indices, Galerkin projection, least-squares reconstruction of eliminated states), the D(ŷ) initial condition (still required by DAEProblem with build_initializeprob = false; guesses do not supply du0), and the POD API.

A new regression in test/array_structure.jl checks that a forcing term present on every grid row (.+ sin(t)) yields the same reduced graph for 4 and 12 states; a naive lazy G*forcing with one entry per row would fail it.

Verification

Julia 1.12.7, ModelingToolkit 11.41.0, MethodOfLines 1.3.0, Symbolics 7.39.0.

GROUP=Core julia +1 --project -e 'using Pkg; Pkg.test()'

Core/DataReduction.jl   |   15     15  15.9s
Core/array_structure.jl |    7      7  34.0s
Core/deim.jl            |   58     58  3m58.3s
Core/full_order.jl      |   32     32  4.1s
Core/interface.jl       |   30     30  0.4s
     Testing ModelOrderReduction tests passed

GROUP=QA julia +1 --project -e 'using Pkg; Pkg.test()'

QA/jet_tests.jl |    5      5  13.4s
QA/qa.jl        |   21     21  3m57.4s
     Testing ModelOrderReduction tests passed

GKSwstype=100 julia +1 --project=docs docs/make.jl

[ Info: Doctest: running doctests.
┌ Warning: Documenter could not auto-detect the building environment. Skipping deployment.
exit status 0 (full build with the executed FitzHugh-Nagumo tutorial, doctests, and link checks)

julia +1 --project=@runic -m Runic --check ., typos ., and git diff --check all exit 0.

New capability: failing before, passing after

deim(::ODEProblem, sol, k) on unmodified main (benchmark-style cubic array system from @mtkcompile, ODEProblem + Tsit5):

ERROR: LoadError: MethodError: no method matching deim(::ODEProblem{Vector{Float64}, ...

Same script on this branch:

reduced equations: 1

The new test/deim.jl testsets "ODEProblem sources and scalar code generation" and "nonautonomous linear terms" cover this path, the mtkcompile + ODEProblem construction, and exact recovery (relative error < 1e-6) of a rank-two nonautonomous linear system with two POD/DEIM modes. The nonautonomous case also passes on main (relative error 4.2e-11 there); it is regression coverage for the term classification, not a bug fix.

Benchmark

benchmark/array_compilation.jl, fresh process per size, two POD/DEIM modes:

n	offline main	offline branch	construction main	construction branch	rhs main	rhs branch	reconstruction
32	12.05 s	9.24 s	4.33 s	4.81 s	0.23 s	0.63 s	0.24 s
512	12.28 s	9.58 s	4.32 s	4.90 s	0.23 s	0.71 s	0.23 s
2048	13.39 s	11.88 s	4.50 s	4.36 s	0.22 s	0.59 s	0.23 s

Offline reduction is faster than main because the initial derivative is evaluated numerically from the projection instead of by symbolic substitution into the scalarized right-hand side. The first residual evaluation is about 0.4 s slower: the generated code now performs the A*ŷ, C*F, and V_S*ŷ array products with parameter matrices instead of inlined scalar sums. All numbers are grid-independent, and both the FitzHugh-Nagumo graph-size regression (reduced_system_size(5) == reduced_system_size(20)) and the new uniform-forcing regression pass. Branch numbers were measured while the Core/QA/docs runs shared the machine.

Not verified

• Julia LTS/prerelease and the downgrade configuration were not run locally.
• MethodOfLines sol[field] convenience indexing still goes through Array observed state is misclassified as parameter-only by getu ModelingToolkit.jl#5072; tests use SymbolicIndexingInterface.observed for non-PDE fields, as before.
• DAEProblem(rom, nothing, tspan) with the default initialization still fails inside MTK (discover_maybe_zeros on array initial conditions); not addressed here.

Review notes

• Version stays at 1.0.0: it is not registered yet (General has 0.1.11), so removing the scalar low-level method is not a post-release breaking change. If custom projection matrices should stay a supported entry point, a deim(sys, snapshot, V, U) method is a small follow-up on top of _reduced_system.
• Term classification is unchanged in effect (numeric-coefficient linear terms projected exactly; everything else state-dependent DEIM-sampled) but now goes through semilinear_form, so symbolic-coefficient linear terms come back as symbolic matrix entries and are moved to the nonlinear part explicitly.

Links

• Keep POD-DEIM ROMs array-based for v1 #178
• Array observed state is misclassified as parameter-only by getu ModelingToolkit.jl#5072

🤖 Generated with Claude Code (model: claude-fable-5-1)

https://claude.ai/code/session_01H5RuLDtEazwTKPkNbPeZVu