Keep POD-DEIM ROMs array-based for v1

AGENTS.md

@@ -0,0 +1,14 @@
+# ModelOrderReduction development
+
+- High-level `deim` must preserve array-valued reduced dynamics and array-valued field
+  reconstruction. Keep full-grid reconstruction coefficients in array parameters;
+  scalar reconstruction equations make the generated graph grow with the grid.
+- Compilation complexity claims apply at fixed reduced dimensions, field count, forcing
+  expression count, and nonlinear stencil size. Offline reduction, problem construction, reconstruction
+  storage, and field evaluation can still scale with the full-order dimension.
+- Validate changes with `GROUP=Core julia --project -e 'using Pkg; Pkg.test()'` and
+  `GROUP=QA julia --project -e 'using Pkg; Pkg.test()'`. Build documentation with
+  `julia --project=docs docs/make.jl` when changing public behavior or examples.
+- Array-form regression coverage must check both generated graph size across grids
+  and numerical dynamics/reconstruction. Equation counts alone do not establish
+  numerical correctness or constant compilation cost.

Project.toml

@@ -1,14 +1,15 @@
 name = "ModelOrderReduction"
 uuid = "207801d6-6cee-43a9-ad0c-f0c64933efa0"
 authors = ["Bowen S. Zhu <bowenzhu@mit.edu> and contributors"]
-version = "0.1.11"
+version = "1.0.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
@@ -17,9 +18,10 @@ TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 [compat]
 DocStringExtensions = "0.9.5"
 LinearAlgebra = "1"
-ModelingToolkit = "11.23.0"
+ModelingToolkit = "11.39.1"
 PrecompileTools = "1.2.1"
 RandomizedLinAlg = "0.1"
+SciMLBase = "3.48"
 Setfield = "1.1.2"
 SparseArrays = "1"
 SymbolicIndexingInterface = "0.3.49"

README.md

@@ -31,10 +31,11 @@ the documentation, which contains the unreleased features.
 
 ```julia
 using ModelingToolkit, MethodOfLines, DifferentialEquations, ModelOrderReduction
-using SciMLBase: symbolic_discretize
+using SciMLBase: discretize
 
 # firstly construct a ModelingToolkit.PDESystem for the FitzHugh-Nagumo model
-@variables x t v(..) w(..)
+@independent_variables x t
+@variables v(..) w(..)
 Dx = Differential(x)
 Dxx = Dx^2
 Dt = Differential(t)
@@ -53,25 +54,22 @@ ivs = [x, t]
 dvs = [v(x, t), w(x, t)]
 pde_sys = PDESystem(eqs, bcs, domains, ivs, dvs; name = Symbol("FitzHugh-Nagumo"))
 
-# transfer to a ModelingToolkit.ODESystem by automated discretization via MethodOfLines
+# discretize to MethodOfLines v1's array-form DAEProblem
 N = 15 # equidistant discretization intervals
 dx = (L - 0.0) / N
 dxs = [x => dx]
 discretization = MOLFiniteDifference(dxs, t)
-ode_sys, tspan = symbolic_discretize(pde_sys, discretization)
-simp_sys = structural_simplify(ode_sys)
-ode_prob = ODEProblem(simp_sys, nothing, tspan)
+full_prob = discretize(pde_sys, discretization; fallback = false)
 
 # solve the full-order model to get snapshots
-sol = solve(ode_prob, Tsit5())
-snapshot_simpsys = Array(sol.original_sol)
+sol = solve(full_prob)
 
 # set POD and DEIM dimensions
 # apply POD-DEIM to obtain the reduced-order model
 pod_dim = deim_dim = 5
-deim_sys = deim(simp_sys, snapshot_simpsys, pod_dim; deim_dim = deim_dim)
-deim_prob = ODEProblem(deim_sys, nothing, tspan)
-deim_sol = solve(deim_prob, Tsit5())
+deim_sys = deim(full_prob, sol, pod_dim; deim_dim = deim_dim)
+deim_prob = DAEProblem(deim_sys, nothing, full_prob.tspan; build_initializeprob = false)
+deim_sol = solve(deim_prob)
 
 # retrieve the approximate solution of the original full-order model
 sol_deim_x = deim_sol[x]
@@ -82,3 +80,8 @@ sol_deim_w = deim_sol[w(x, t)]
 The following figure shows the comparison of the solutions of the 32-dimension full-order model and the POD5-DEIM5 reduced-order model.
 
 ![comparison](https://user-images.githubusercontent.com/45696147/195765614-df9092a2-4fca-4602-bb15-81e65b2b572e.svg)
+
+For fixed reduced dimensions, field and forcing structure, and local nonlinear stencils,
+the reduced dynamics and field reconstruction use a grid-independent symbolic array graph. Offline reduction
+and full-field reconstruction still scale with the grid. See the tutorial for the
+algebraic reconstruction approximation and supported DAE scope.

benchmark/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+ModelOrderReduction = "207801d6-6cee-43a9-ad0c-f0c64933efa0"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
+
+[compat]
+ModelOrderReduction = "1"
+ModelingToolkit = "11"
+OrdinaryDiffEq = "7"
+SymbolicIndexingInterface = "0.3"

benchmark/README.md

@@ -0,0 +1,20 @@
+Run from the repository root:
+
+```sh
+julia --project=benchmark -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+julia --project=benchmark benchmark/array_compilation.jl 8 32 128
+```
+
+POD and DEIM dimensions stay fixed at two. The script reports offline reduction,
+problem construction, first residual evaluation, and first field reconstruction
+separately. Compiler time is Julia's `@timed` compilation counter; total elapsed time
+includes work other than compilation. Package loading occurs before the measurements.
+
+To compare fresh-process compilation across grids, invoke the second command once
+per size. Passing several sizes measures the first case cold and subsequent cases
+with reusable compilation already warm. Equation and observable counts describe the
+symbolic representation; timings are empirical evidence, not an asymptotic proof.
+
+Reconstruction uses `SymbolicIndexingInterface.observed`, retaining the whole array
+expression. The convenience `getu` path currently encounters
+https://github.com/SciML/ModelingToolkit.jl/issues/5072 for these observables.

benchmark/array_compilation.jl

@@ -0,0 +1,31 @@
+using ModelOrderReduction, ModelingToolkit, OrdinaryDiffEq
+import SymbolicIndexingInterface as SII
+@independent_variables t
+D = Differential(t)
+function trial(n)
+    @variables z(t)[1:n]
+    states = Symbolics.unwrap.(Symbolics.scalarize(z))
+    sys = System([D(z) ~ -z - z .^ 3], t, states, []; name = :scaling)
+    snapshot = [sin(0.13 * i * j) + cos(0.07 * i * (j + 1)) for i in 1:n, j in 1:8]
+    offline = @timed deim(sys, snapshot, 2)
+    rom = offline.value
+    construction = @timed DAEProblem(rom, nothing, (0.0, 0.1); build_initializeprob = false)
+    prob = construction.value
+    du = zeros(2); residual = zeros(2)
+    rhs = @timed prob.f(residual, du, prob.u0, prob.p, 0.0)
+    reconstruction = @timed SII.observed(prob, z)(prob.u0, prob.p, first(prob.tspan))
+    println(
+        (
+            n = n, equations = length(equations(rom)), observables = length(observed(rom)),
+            offline_seconds = offline.time, construction_seconds = construction.time,
+            construction_compile_seconds = construction.compile_time,
+            rhs_seconds = rhs.time, rhs_compile_seconds = rhs.compile_time,
+            reconstruction_seconds = reconstruction.time,
+            reconstruction_compile_seconds = reconstruction.compile_time,
+        )
+    )
+    return flush(stdout)
+end
+for n in parse.(Int, ARGS)
+    trial(n)
+end

docs/Project.toml

@@ -13,8 +13,8 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 DifferentialEquations = "8"
 Documenter = "1"
 LaTeXStrings = "1.3"
-MethodOfLines = "0.11, 0.12"
-ModelOrderReduction = "0.1"
+MethodOfLines = "1"
+ModelOrderReduction = "1"
 ModelingToolkit = "11"
 Plots = "1.40"
 SciMLBase = "3"

docs/src/tutorials/deim_FitzHugh-Nagumo.md

@@ -43,7 +43,8 @@ describing the original FitzHugh-Nagumo model.
 ```@example deim_FitzHugh_Nagumo
 using ModelingToolkit
 using LaTeXStrings
-@variables x t v(..) w(..)
+@independent_variables x t
+@variables v(..) w(..)
 Dx = Differential(x)
 Dxx = Dx^2
 Dt = Differential(t)
@@ -69,28 +70,27 @@ nothing # hide
 ```
 
 Next, we apply finite difference discretization using
-[MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/).
+[MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/). MethodOfLines v1 constructs
+an array-form `DAEProblem` directly. Setting `fallback = false` keeps this path explicit and
+avoids compiling a scalarized full-order ODE system.
 
 ```@example deim_FitzHugh_Nagumo
 using MethodOfLines
-using SciMLBase: symbolic_discretize
+using SciMLBase: discretize
 N = 15 # equidistant discretization intervals
 dx = (L - 0.0) / N
 dxs = [x => dx]
 order = 2
 discretization = MOLFiniteDifference(dxs, t; approx_order = order)
-ode_sys, tspan = symbolic_discretize(pde_sys, discretization)
-simp_sys = mtkcompile(ode_sys)
-ode_prob = ODEProblem(simp_sys, nothing, tspan;
-    missing_guess_value = ModelingToolkit.MissingGuessValue.Constant(0.0))
+full_prob = discretize(pde_sys, discretization; fallback = false)
 nothing # hide
 ```
 
-The snapshot trajectories are obtained by solving the full-order system.
+The snapshot trajectories are obtained by solving the full-order DAE with its default solver.
 
 ```@example deim_FitzHugh_Nagumo
 using DifferentialEquations
-sol = solve(ode_prob, Tsit5())
+sol = solve(full_prob)
 sol_x = sol[x]
 nₓ = length(sol_x) # number of discretization points in x
 nₜ = length(sol[t]) # number of discretization points in time
@@ -136,18 +136,35 @@ terms. This can be done by simply calling [`deim`](@ref).
 
 ```@example deim_FitzHugh_Nagumo
 using ModelOrderReduction
-snapshot_simpsys = Array(sol.original_sol)
 pod_dim = deim_dim = 5
-deim_sys = deim(simp_sys, snapshot_simpsys, pod_dim)
-deim_prob = ODEProblem(deim_sys, nothing, tspan)
-deim_sol = solve(deim_prob, Tsit5())
+deim_sys = deim(full_prob, sol, pod_dim; deim_dim = deim_dim)
+deim_prob = DAEProblem(deim_sys, nothing, full_prob.tspan; build_initializeprob = false)
+deim_sol = solve(deim_prob)
 nₜ_deim = length(deim_sol[t])
 sol_deim_x = deim_sol[x]
 sol_deim_v = deim_sol[v(x, t)]
 sol_deim_w = deim_sol[w(x, t)]
 nothing # hide
 ```
 
+The ROM retains one array differential equation and one reconstruction array per PDE
+field. For fixed POD/DEIM dimensions, field count, forcing expression count, and local nonlinear stencil size,
+the generated symbolic graph is independent of the spatial grid. Use `DAEProblem`
+for array code generation; scalarizing the returned system loses that property.
+
+This compilation bound does not cover offline structural elimination, POD/DEIM training,
+problem construction, or full-field evaluation. The reconstruction matrix requires
+`O(nk)` storage. Differential states use the POD lift; eliminated algebraic states use
+a least-squares fit to the training snapshots, so reconstructed boundary constraints
+are approximate. The DAE must reduce to an explicit first-order ODE after algebraic
+elimination. General descriptor-DAE projection is not implemented by this method.
+
+For a symbolic array field outside the PDE solution interface, construct its array
+reconstruction function with `SymbolicIndexingInterface.observed(deim_sys, field)`
+and evaluate it with `(reduced_state, parameters, time)`. Convenience `sol[field]`
+indexing currently has an [upstream array-observable issue](https://github.com/SciML/ModelingToolkit.jl/issues/5072).
+The PDE field indexing used in this tutorial uses the explicit reconstruction path.
+
 And plot the result from the POD-DEIM reduced system.
 
 ```@example deim_FitzHugh_Nagumo
@@ -184,7 +201,7 @@ plt_2 = plot(xlabel = L"v(x,t)", ylabel = L"x", zlabel = L"w(x,t)", xlims = (-0.
     ylims = (0.0, L), zlims = (0.0, 0.25), xflip = true, camera = (50, 30),
     titlefont = 10, title = "Comparison of full and reduced systems")
 plot!(plt_2, unconnected(snapshot_v), unconnected(sol_x, nₜ), unconnected(snapshot_w),
-    label = "Full$(length(ModelingToolkit.get_eqs(ode_sys)))")
+    label = "Full$(length(full_prob.u0))")
 plot!(plt_2, unconnected(sol_deim_v), unconnected(sol_deim_x, nₜ_deim),
     unconnected(sol_deim_w), label = "POD$(pod_dim)/DEIM$(deim_dim)")
 ```

src/ModelOrderReduction.jl

@@ -2,9 +2,10 @@ module ModelOrderReduction
 
 using DocStringExtensions: DocStringExtensions, FUNCTIONNAME, SIGNATURES, TYPEDSIGNATURES
 
-using ModelingToolkit: ModelingToolkit, @variables, Differential, Equation, Num, ODESystem,
-    SymbolicUtils, Symbolics, arguments, build_function, complete,
-    expand, substitute, tearing_substitution
+using ModelingToolkit: ModelingToolkit, @parameters, @variables, Differential, Equation, Num, ODESystem,
+    SymbolicUtils, Symbolics, arguments, complete,
+    expand, substitute, mtkcompile
+using SciMLBase: SciMLBase
 using SymbolicIndexingInterface: SymbolicIndexingInterface
 using LinearAlgebra: LinearAlgebra, /, \, mul!, qr, svd