From 2455dd6ea8b32063ea839d5ab58b209a2c5cd57d Mon Sep 17 00:00:00 2001 From: Codex Date: Sat, 5 Sep 2026 06:55:18 -0400 Subject: [PATCH] Keep POD-DEIM ROMs in array form for v1 Return array dynamics and field reconstruction from both high-level entry points. Evaluate training snapshots and consistent initial derivatives offline, retain the training parameter point, and substitute eliminated algebraic variables before projection. Add numerical and graph-size regressions and a compilation benchmark. Document the offline complexity boundary and fitted algebraic reconstruction. Track the independent array-indexing issue at https://github.com/SciML/ModelingToolkit.jl/issues/5072. Co-Authored-By: Chris Rackauckas Co-Authored-By: Codex Agent-Harness: Codex CLI 0.153.4 Agent-Model: gpt-6-astra Agent-Session: /home/crackauc/.codex/sessions/2026/09/05/rollout-2026-09-05T05-19-42-01a070dd-d97a-7180-9b71-095f2a5f5489.jsonl --- AGENTS.md | 14 + Project.toml | 6 +- README.md | 25 +- benchmark/Project.toml | 11 + benchmark/README.md | 20 + benchmark/array_compilation.jl | 31 ++ docs/Project.toml | 4 +- docs/src/tutorials/deim_FitzHugh-Nagumo.md | 45 ++- src/ModelOrderReduction.jl | 7 +- src/deim.jl | 432 ++++++++++++++++++--- src/utils.jl | 5 +- test/Project.toml | 4 +- test/array_structure.jl | 31 ++ test/deim.jl | 264 ++++++++++++- 14 files changed, 795 insertions(+), 104 deletions(-) create mode 100644 AGENTS.md create mode 100644 benchmark/Project.toml create mode 100644 benchmark/README.md create mode 100644 benchmark/array_compilation.jl create mode 100644 test/array_structure.jl diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 00000000..4917b81e --- /dev/null +++ b/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. diff --git a/Project.toml b/Project.toml index c06c3532..4f3abd70 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelOrderReduction" uuid = "207801d6-6cee-43a9-ad0c-f0c64933efa0" authors = ["Bowen S. Zhu and contributors"] -version = "0.1.11" +version = "1.0.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -9,6 +9,7 @@ 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" diff --git a/README.md b/README.md index d7d4d8a0..8b953396 100644 --- a/README.md +++ b/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. diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 00000000..616a2fc5 --- /dev/null +++ b/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" diff --git a/benchmark/README.md b/benchmark/README.md new file mode 100644 index 00000000..a25737b7 --- /dev/null +++ b/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. diff --git a/benchmark/array_compilation.jl b/benchmark/array_compilation.jl new file mode 100644 index 00000000..42314315 --- /dev/null +++ b/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 diff --git a/docs/Project.toml b/docs/Project.toml index b256827f..148351b3 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/tutorials/deim_FitzHugh-Nagumo.md b/docs/src/tutorials/deim_FitzHugh-Nagumo.md index af11ecb7..fce0bcc7 100644 --- a/docs/src/tutorials/deim_FitzHugh-Nagumo.md +++ b/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,11 +136,10 @@ 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)] @@ -148,6 +147,24 @@ 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)") ``` diff --git a/src/ModelOrderReduction.jl b/src/ModelOrderReduction.jl index b94b1559..b7ef2f80 100644 --- a/src/ModelOrderReduction.jl +++ b/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 diff --git a/src/deim.jl b/src/deim.jl index bbce3651..e12f85fd 100644 --- a/src/deim.jl +++ b/src/deim.jl @@ -156,6 +156,268 @@ function deim( reduced_rhss = Â * ŷ + ĝ + F̂ return reduced_rhss, linear_projection_eqs end + +function _array_state_groups(unknowns) + order = Any[] + rows_by_variable = Dict{Any, Vector{Int}}() + indices_by_variable = Dict{Any, Vector{Tuple}}() + array_variables = Set{Any}() + for (row, unknown) in enumerate(unknowns) + value = Symbolics.unwrap(unknown) + if SymbolicUtils.iscall(value) && SymbolicUtils.operation(value) === getindex + index_arguments = SymbolicUtils.arguments(value) + variable = first(index_arguments) + indices = Tuple(SymbolicUtils.unwrap_const.(index_arguments[2:end])) + push!(array_variables, variable) + else + variable = value + indices = () + end + if !haskey(rows_by_variable, variable) + rows_by_variable[variable] = Int[] + indices_by_variable[variable] = Tuple[] + push!(order, variable) + end + push!(rows_by_variable[variable], row) + push!(indices_by_variable[variable], indices) + end + + return map(order) do variable + rows = rows_by_variable[variable] + is_array = variable in array_variables + shape = is_array ? size(Symbolics.wrap(variable)) : () + is_array && prod(shape) != length(rows) && throw( + DimensionMismatch( + "array variable $variable has shape $shape but $(length(rows)) scalar elements" + ) + ) + if is_array + ordered_rows = Vector{Int}(undef, prod(shape)) + assigned = falses(prod(shape)) + linear_indices = LinearIndices(shape) + for (row, indices) in zip(rows, indices_by_variable[variable]) + all(index -> index isa Integer, indices) || throw( + ArgumentError("array variable $variable has non-integer indices $indices") + ) + linear_index = linear_indices[indices...] + assigned[linear_index] && throw( + ArgumentError("array variable $variable repeats index $indices") + ) + ordered_rows[linear_index] = row + assigned[linear_index] = true + end + all(assigned) || throw( + ArgumentError("array variable $variable does not contain every array element") + ) + rows = ordered_rows + end + return (; variable, rows, shape, is_array) + end +end + +function _evaluate_nonlinear_snapshot!( + output::AbstractMatrix, expressions::AbstractVector, + variables::AbstractVector, snapshot::AbstractMatrix; + parameter_values::AbstractDict = Dict{Any, Any}(), + iv = nothing, snapshot_times = nothing, kwargs... + ) + replacements = Dict{Any, Any}(parameter_values) + for column in axes(snapshot, 2) + isnothing(snapshot_times) || (replacements[iv] = snapshot_times[column]) + for (variable, value) in zip(variables, @view(snapshot[:, column])) + replacements[variable] = value + replacements[Num(variable)] = value + end + for row in axes(output, 1) + value = substitute( + expressions[row], replacements; fold = Val(true), kwargs... + ) + value = SymbolicUtils.unwrap_const(Symbolics.unwrap(value)) + value isa Number || throw( + ArgumentError( + "nonlinear expression $(expressions[row]) did not evaluate to a number; provide numeric parameter defaults and snapshot_times for time-dependent expressions" + ) + ) + output[row, column] = value + end + end + return output +end + +function _array_deim( + source_system::ODESystem, source_snapshot::AbstractMatrix, + pod_dim::Integer, deim_dim::Integer, name::Symbol; + snapshot_times = nothing, training_parameters = Dict{Any, Any}(), kwargs... + ) + source_system = deepcopy(source_system) + source_unknowns = ModelingToolkit.unknowns(source_system) + number_of_source_unknowns = length(source_unknowns) + size(source_snapshot, 1) == number_of_source_unknowns || throw( + DimensionMismatch( + "the snapshot has $(size(source_snapshot, 1)) rows, but the source system has $number_of_source_unknowns unknowns" + ) + ) + isnothing(snapshot_times) || length(snapshot_times) == size(source_snapshot, 2) || + throw(DimensionMismatch("snapshot_times must contain one time per snapshot column")) + state_groups = _array_state_groups(source_unknowns) + source_observed = ModelingToolkit.observed(source_system) + + source_deqs, source_constraints = get_deqs(source_system) + if isempty(source_constraints) && length(source_deqs) == length(source_unknowns) + sys = source_system + else + ModelingToolkit.iscomplete(source_system) && @set! source_system.complete = false + sys = mtkcompile(source_system) + end + iv = ModelingToolkit.independent_variable(sys) + D = Differential(iv) + dynamic_unknowns = ModelingToolkit.unknowns(sys) + model_parameters = filter(ModelingToolkit.parameters(sys)) do parameter + !ModelingToolkit.isinitial(parameter) + end + model_parameter_keys = Set( + Symbolics.unwrap(parameter) for parameter in model_parameters + ) + system_defaults = copy(ModelingToolkit.initial_conditions(sys)) + merge!(system_defaults, training_parameters) + parameter_values = Dict{Any, Any}() + for parameter in model_parameters + key = Symbolics.unwrap(parameter) + haskey(system_defaults, key) || continue + value = system_defaults[key] + parameter_values[key] = value + parameter_values[Symbolics.wrap(key)] = value + end + differential_equations, algebraic_equations = get_deqs(ModelingToolkit.full_equations(sys)) + isempty(algebraic_equations) || throw( + ArgumentError( + "the source system must structurally compile to an explicit first-order ODE" + ) + ) + equations_by_unknown = Dict{Any, Equation}() + for equation in differential_equations + unknown = only(SymbolicUtils.arguments(Symbolics.unwrap(equation.lhs))) + haskey(equations_by_unknown, unknown) && throw( + ArgumentError("the compiled system has duplicate equations for $unknown") + ) + equations_by_unknown[unknown] = equation + end + length(equations_by_unknown) == length(dynamic_unknowns) || throw( + ArgumentError( + "the compiled system does not have one differential equation per unknown" + ) + ) + source_rows = Dict( + Symbolics.unwrap(unknown) => row for (row, unknown) in enumerate(source_unknowns) + ) + dynamic_rows = map(dynamic_unknowns) do unknown + row = get(source_rows, Symbolics.unwrap(unknown), 0) + iszero(row) && throw( + ArgumentError( + "compiled unknown $unknown is not present in the source system" + ) + ) + return row + end + dynamic_snapshot = Matrix{Float64}(source_snapshot[dynamic_rows, :]) + + state_reducer = POD(dynamic_snapshot, pod_dim) + reduce!(state_reducer, TSVD()) + state_basis = state_reducer.rbasis + + reduced_name = gensym(:ŷ) + reduced_state = (@variables $reduced_name(iv)[1:pod_dim])[1] + + right_hand_sides = map(dynamic_unknowns) do unknown + equation = get(equations_by_unknown, Symbolics.unwrap(unknown), nothing) + isnothing(equation) && throw( + ArgumentError("the compiled system has no differential equation for $unknown") + ) + equation.rhs + end + linear_coefficients, constant_part, nonlinear_part = separate_terms( + right_hand_sides, dynamic_unknowns, iv + ) + + nonlinear_snapshot = similar(dynamic_snapshot) + _evaluate_nonlinear_snapshot!( + nonlinear_snapshot, nonlinear_part, dynamic_unknowns, dynamic_snapshot; + parameter_values, iv, snapshot_times, kwargs... + ) + nonlinear_reducer = POD(nonlinear_snapshot, deim_dim) + reduce!(nonlinear_reducer, TSVD()) + nonlinear_basis = nonlinear_reducer.rbasis + + reduced_rhs, _ = deim( + dynamic_unknowns, linear_coefficients, constant_part, nonlinear_part, + reduced_state, state_basis, nonlinear_basis; kwargs... + ) + reduced_equation = D(reduced_state) ~ reduced_rhs + + reduced_snapshot = state_basis' * dynamic_snapshot + reconstruction_basis = Matrix{Float64}(source_snapshot) / reduced_snapshot + reconstruction_basis[dynamic_rows, :] = state_basis + reconstruction_rows = reduce(vcat, (collect(group.rows) for group in state_groups)) + reconstruction_basis = reconstruction_basis[reconstruction_rows, :] + basis_name = gensym(:reconstruction_basis) + basis_parameter = ( + @parameters $basis_name[1:number_of_source_unknowns, 1:pod_dim] + )[1] + basis_parameter = ModelingToolkit.setdefault( + basis_parameter, reconstruction_basis + ) + state_replacements = Dict{Any, Any}() + next_reconstruction_row = 1 + reconstruction_equations = map(state_groups) do group + rows = next_reconstruction_row:(next_reconstruction_row + length(group.rows) - 1) + next_reconstruction_row = last(rows) + 1 + reconstruction = Symbolics.wrap(basis_parameter[rows, :] * reduced_state) + value = group.is_array ? reshape(reconstruction, group.shape) : + only(Symbolics.scalarize(reconstruction)) + state_replacements[group.variable] = value + Symbolics.wrap(group.variable) ~ value + end + source_state_variables = Set(Symbolics.unwrap.(source_unknowns)) + union!(source_state_variables, (group.variable for group in state_groups)) + preserved_observed = Equation[] + for equation in source_observed + Symbolics.unwrap(equation.lhs) in source_state_variables && continue + push!( + preserved_observed, + equation.lhs ~ substitute(equation.rhs, state_replacements; kwargs...) + ) + end + reduced_observed = _sort_observed_equations( + [ + reconstruction_equations; preserved_observed + ] + ) + reduced_parameters = [model_parameters; Symbolics.unwrap(basis_parameter)] + + initial_conditions = system_defaults + filter!(pair -> first(pair) in model_parameter_keys, initial_conditions) + initial_conditions[Symbolics.unwrap(reduced_state)] = reduced_snapshot[:, 1] + initial_derivative = zeros(eltype(reduced_snapshot), pod_dim, 1) + _evaluate_nonlinear_snapshot!( + initial_derivative, Symbolics.scalarize(Symbolics.wrap(reduced_rhs)), + Symbolics.scalarize(reduced_state), + reduced_snapshot[:, 1:1]; parameter_values, iv, snapshot_times, kwargs... + ) + initial_conditions[Symbolics.unwrap(D(reduced_state))] = vec(initial_derivative) + initial_conditions[Symbolics.unwrap(basis_parameter)] = reconstruction_basis + reduced_system = ModelingToolkit.System( + [reduced_equation], iv, Symbolics.unwrap.(Symbolics.scalarize(reduced_state)), + reduced_parameters; name, observed = reduced_observed, initial_conditions + ) + for key in (ModelingToolkit.ProblemTypeCtx, ModelingToolkit.MiscSystemData) + SymbolicUtils.hasmetadata(source_system, key) || continue + reduced_system = SymbolicUtils.setmetadata( + reduced_system, key, SymbolicUtils.getmetadata(source_system, key, nothing) + ) + end + return complete(reduced_system) +end + """ $(FUNCTIONNAME)( sys::ModelingToolkit.ODESystem, @@ -163,6 +425,7 @@ end pod_dim::Integer; deim_dim::Integer = pod_dim, name::Symbol = Symbol(nameof(sys), :_deim), + snapshot_times = nothing, kwargs... ) -> ModelingToolkit.ODESystem @@ -174,28 +437,37 @@ the Discrete Empirical Interpolation Method (DEIM). The LHS of equations in `sys` are all assumed to be 1st order derivatives. Use `ModelingToolkit.ode_order_lowering` to transform higher order ODEs before applying DEIM. -`sys` is assumed to have no internal systems. End users are encouraged to call -`ModelingToolkit.mtkcompile` beforehand. +`sys` is assumed to have no internal systems. MOR performs structural compilation as needed. +Compilation may scalarize the input during offline analysis, but the returned reduced +dynamics are always represented by one symbolic array equation. MethodOfLines v1 +`symbolic_discretize` systems should be passed before structural compilation so their full +array reconstruction metadata remains available. The rows of `snapshot` must follow +`ModelingToolkit.unknowns(sys)`. Nonlinear training uses numeric parameter defaults. +Use the problem/solution overload for time-dependent nonlinear terms and problem-specific +parameter values. Construct a `DAEProblem` from the returned system to retain array +code generation. Eliminated algebraic states use a least-squares reconstruction from +the training snapshots. The POD basis used for DEIM interpolation is obtained from the snapshot matrix of the -nonlinear terms, which is computed by executing the runtime-generated function for -nonlinear expressions. +nonlinear terms. It is evaluated symbolically so the offline reduction never generates a +full-order scalar function. # Arguments -- `sys::ModelingToolkit.ODESystem`: compiled first-order ODE system without internal - subsystems. +- `sys::ModelingToolkit.ODESystem`: first-order system without internal subsystems. Pass + MethodOfLines `symbolic_discretize` output before structural compilation. - `snapshot::AbstractMatrix`: state-by-time snapshot matrix for `sys`. - `pod_dim::Integer`: number of POD state modes to retain. # Keywords - `deim_dim::Integer = pod_dim`: number of DEIM modes for nonlinear terms. - `name::Symbol = Symbol(nameof(sys), :_deim)`: name assigned to the reduced system. -- `kwargs...`: keyword arguments forwarded to ModelingToolkit transformations and - generated nonlinear functions. +- `snapshot_times = nothing`: saved times corresponding to snapshot columns; required for + time-dependent expressions. The problem/solution overload supplies these automatically. +- `kwargs...`: keyword arguments forwarded to symbolic substitutions. # Examples ```julia -reduced_system = deim(compiled_system, snapshots, 4; deim_dim = 6) +reduced_system = deim(source_system, snapshots, 4; deim_dim = 6) ``` # Returns @@ -203,63 +475,115 @@ reduced_system = deim(compiled_system, snapshots, 4; deim_dim = 6) # Throws - `ArgumentError`: if the observed equations contain a dependency cycle. +- `DimensionMismatch`: if a source system and its snapshot have different state + dimensions. """ function deim( sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer; deim_dim::Integer = pod_dim, name::Symbol = Symbol(nameof(sys), :_deim), + snapshot_times = nothing, kwargs... )::ODESystem - sys = deepcopy(sys) - @set! sys.name = name + return _array_deim( + sys, snapshot, pod_dim, deim_dim, name; snapshot_times, kwargs... + ) +end - # handle ODESystem.substitutions - # https://github.com/SciML/ModelingToolkit.jl/issues/1754 - sys = tearing_substitution(sys; kwargs...) +""" + $(FUNCTIONNAME)( + prob::SciMLBase.DAEProblem, + sol, + pod_dim::Integer; + deim_dim::Integer = pod_dim, + name::Union{Nothing, Symbol} = nothing, + kwargs... + ) -> ModelingToolkit.ODESystem - iv = ModelingToolkit.get_iv(sys) # the single independent variable - D = Differential(iv) - dvs = ModelingToolkit.get_unknowns(sys) # dependent variables - - pod_reducer = POD(snapshot, pod_dim) - reduce!(pod_reducer, TSVD()) - V = pod_reducer.rbasis # POD basis - - 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) - - # generate an in-place function from the symbolic expression of the nonlinear functions - F_func! = build_function(F, dvs; expression = Val{false}, kwargs...)[2] - nonlinear_snapshot = similar(snapshot) # snapshot matrix of nonlinear terms - for i in 1:size(snapshot, 2) # iterate through time instances - F_func!(view(nonlinear_snapshot, :, i), view(snapshot, :, i)) - end +Reduce a first-order symbolic `DAEProblem` and one of its saved solutions using POD-DEIM. - deim_reducer = POD(nonlinear_snapshot, deim_dim) - reduce!(deim_reducer, TSVD()) - U = deim_reducer.rbasis # DEIM projection basis +This method is intended for array-form problems such as those produced by MethodOfLines v1. +It supports DAEs that become explicit first-order ODEs after algebraic elimination; +it does not project general descriptor or differential-algebraic residual systems. +The full-order solve retains its array-form DAE compilation path. For the symbolic reduction, +ModelingToolkit tearing eliminates algebraic unknowns and isolates the differential equations +without constructing or compiling a full-order `ODEProblem`. The returned system has one array differential equation and one reconstruction array per +source field. Construct a `DAEProblem` with `build_initializeprob = false` from it +to retain array code generation. Consistent state and derivative initial values are +computed at the first training snapshot. - reduced_rhss, linear_projection_eqs = deim(dvs, A, g, F, ŷ, V, U; kwargs...) +For fixed POD/DEIM dimensions, field count, forcing expression count, and bounded nonlinear stencil size, the +returned symbolic graph is independent of the full-order grid. Offline structural +elimination, snapshot evaluation, basis construction, and reconstruction storage still +scale with grid size. Eliminated algebraic states are reconstructed by a least-squares +fit to the training snapshots, so their constraints need not hold exactly away from +those snapshots. Parameters are trained at `prob`'s values, which become the ROM defaults. - reduced_deqs = D.(ŷ) ~ reduced_rhss - @set! sys.eqs = [Symbolics.scalarize(reduced_deqs); eqs] +`sol` may be the `SciMLBase.PDETimeSeriesSolution` returned by MethodOfLines or its underlying +`SciMLBase.AbstractODESolution`. - old_observed = ModelingToolkit.get_observed(sys) - new_observed = [old_observed; linear_projection_eqs] - @set! sys.observed = _sort_observed_equations(new_observed) +# Arguments +- `prob::SciMLBase.DAEProblem`: symbolic first-order DAE problem whose function stores its + ModelingToolkit system in `prob.f.sys`. +- `sol`: saved solution obtained from `prob`. +- `pod_dim::Integer`: number of POD state modes to retain. - # Numeric initial conditions for the reduced unknowns from the snapshot's first column. - # The snapshot is assumed to start at t = tspan[1], matching the FOM initial state. - new_ics = copy(ModelingToolkit.get_initial_conditions(sys)) - new_ics[Symbolics.unwrap(ŷ)] = V' * snapshot[:, 1] - @set! sys.initial_conditions = new_ics +# Keywords +- `deim_dim::Integer = pod_dim`: number of DEIM modes for nonlinear terms. +- `name::Union{Nothing, Symbol} = nothing`: name assigned to the reduced system. The default + appends `_deim` to the full system's name. +- `kwargs...`: keyword arguments forwarded to ModelingToolkit transformations and symbolic + substitutions. + +# Examples +```julia +full_problem = discretize(pde_system, discretization; fallback = false) +full_solution = solve(full_problem) +reduced_system = deim(full_problem, full_solution, 4) +reduced_problem = DAEProblem( + reduced_system, nothing, full_problem.tspan; build_initializeprob = false +) +``` - return complete(sys) +# Returns +- `ModelingToolkit.ODESystem`: the reduced and completed explicit ODE system. + +# Throws +- `ArgumentError`: if `sol` is not from `prob`, does not start at the beginning of the problem, + or the DAE cannot be torn into an explicit ODE whose unknowns match the saved states. +- `DimensionMismatch`: if the saved solution does not match the DAE system. +""" +function deim( + prob::SciMLBase.DAEProblem, sol::SciMLBase.AbstractODESolution, + pod_dim::Integer; deim_dim::Integer = pod_dim, + name::Union{Nothing, Symbol} = nothing, kwargs... + )::ODESystem + raw_sys = SymbolicIndexingInterface.symbolic_container(prob.f) + raw_sys isa ODESystem || + throw(ArgumentError("the DAE problem must contain a symbolic ModelingToolkit system")) + SymbolicIndexingInterface.symbolic_container(sol.prob.f) === 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.unknowns(raw_sys) + full_snapshot = Array(sol) + size(full_snapshot, 1) == length(raw_variables) || + throw(DimensionMismatch("solution states must match the DAE system unknowns")) + + training_parameters = Dict{Any, Any}() + for parameter in ModelingToolkit.parameters(raw_sys) + ModelingToolkit.isinitial(parameter) && continue + training_parameters[Symbolics.unwrap(parameter)] = prob.ps[parameter] + end + reduced_name = isnothing(name) ? Symbol(nameof(raw_sys), :_deim) : name + return _array_deim( + raw_sys, full_snapshot, pod_dim, deim_dim, reduced_name; + snapshot_times = sol.t, training_parameters, kwargs... + ) +end + +function deim( + prob::SciMLBase.DAEProblem, sol::SciMLBase.PDETimeSeriesSolution, + pod_dim::Integer; kwargs... + )::ODESystem + return deim(prob, sol.original_sol, pod_dim; kwargs...) end diff --git a/src/utils.jl b/src/utils.jl index d1341750..8ad2ab41 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -12,7 +12,10 @@ the left hand side, which is typically the result of `ModelingToolkit.mtkcompile Equations from subsystems are not included. """ function get_deqs(sys::ODESystem)::Tuple{Vector{Equation}, Vector{Equation}} - eqs = ModelingToolkit.get_eqs(sys) + return get_deqs(ModelingToolkit.equations(sys)) +end + +function get_deqs(eqs::AbstractVector{Equation}) deqs = Equation[] others = Equation[] for eq in eqs diff --git a/test/Project.toml b/test/Project.toml index 7da12dc2..c04515fd 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" @@ -9,8 +10,9 @@ SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] +SymbolicIndexingInterface = "0.3" LinearAlgebra = "1" -MethodOfLines = "0.11, 0.12" +MethodOfLines = "1" ModelingToolkit = "11" OrdinaryDiffEq = "7" SafeTestsets = "0.1" diff --git a/test/array_structure.jl b/test/array_structure.jl new file mode 100644 index 00000000..db9cc823 --- /dev/null +++ b/test/array_structure.jl @@ -0,0 +1,31 @@ +using Test, ModelOrderReduction, ModelingToolkit + +@independent_variables t +@variables u(t)[1:4] +D = Differential(t) +snapshot = [sin(0.3 * i * j) + cos(0.2 * i * (j + 1)) for i in 1:4, j in 1:8] + +@testset "array ROM regressions" begin + @testset "structure" begin + @mtkcompile sys = System([D(u) ~ -u - u .^ 3], t) + rom = deim(sys, snapshot, 2) + @test length(equations(rom)) == 1 + @test length(observed(rom)) == 1 + end + @testset "parameter defaults" begin + @parameters α = 2.0 + @mtkcompile sys = System([D(u) ~ -u - α * u .^ 3], t) + rom = deim(sys, snapshot, 2) + default_value = initial_conditions(rom)[Symbolics.unwrap(α)] + @test SymbolicUtils.unwrap_const(Symbolics.unwrap(default_value)) == 2.0 + end + @testset "scalar input" begin + @variables z₁(t) z₂(t) + @mtkcompile sys = System([D(z₁) ~ z₁ - z₁^3, D(z₂) ~ -z₂], t) + rom = deim(sys, [1.0 0.8 0.6; 0.5 0.4 0.3], 1) + @test length(unknowns(rom)) == 1 + @test isempty(guesses(rom)) + @test isempty(initialization_equations(rom)) + end + +end diff --git a/test/deim.jl b/test/deim.jl index 436cb3c6..e7fa884c 100644 --- a/test/deim.jl +++ b/test/deim.jl @@ -1,9 +1,12 @@ using Test, ModelOrderReduction +import SymbolicIndexingInterface as SII using ModelingToolkit, MethodOfLines, OrdinaryDiffEq -using SciMLBase: symbolic_discretize +using LinearAlgebra: norm, svd +using SciMLBase: successful_retcode, symbolic_discretize # construct an ModelingToolkit.ODESystem with non-empty field substitutions -@variables x t v(..) w(..) +@independent_variables x t +@variables v(..) w(..) Dx = Differential(x) Dxx = Dx^2 Dt = Differential(t) @@ -37,29 +40,258 @@ 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) # field substitutions is non-empty -# The BCs fix every state to 0 at t=0, so provide those initial values directly. -# Leaving them as missing guesses makes MethodOfLines' initialization system -# structurally singular under MTK v11, and the least-squares solve diverges. -u0 = [u => 0.0 for u in ModelingToolkit.get_unknowns(simp_sys)] -ode_prob = ODEProblem(simp_sys, u0, tspan) -sol = solve(ode_prob, Rodas5P(), saveat = 1.0) - -snapshot_simpsys = Array(sol.original_sol) +full_prob = discretize(pde_sys, discretization; fallback = false) +@test full_prob isa DAEProblem +sol = solve(full_prob; saveat = 1.0) +@test successful_retcode(sol) + +snapshot = Array(sol.original_sol) pod_dim = 3 -deim_sys = @test_nowarn deim(simp_sys, snapshot_simpsys, pod_dim) +deim_sys = @test_nowarn deim(full_prob, sol, pod_dim) # check the number of dependent variables in the new system @test length(ModelingToolkit.get_unknowns(deim_sys)) == pod_dim -deim_prob = ODEProblem(complete(deim_sys), nothing, tspan) +function is_symbolic_array(expression) + value = ModelingToolkit.Symbolics.unwrap(expression) + return value isa AbstractArray || + ModelingToolkit.SymbolicUtils.symtype(value) <: AbstractArray +end + +reduced_equations = ModelingToolkit.get_eqs(deim_sys) +@test length(reduced_equations) == 1 +@test is_symbolic_array(only(reduced_equations).lhs) +reconstruction_equations = ModelingToolkit.get_observed(deim_sys) +@test length(reconstruction_equations) == length(dvs) +@test all(reconstruction_equations) do equation + is_symbolic_array(equation.lhs) && is_symbolic_array(equation.rhs) +end +reconstruction_parameters = filter(ModelingToolkit.get_ps(deim_sys)) do parameter + is_symbolic_array(parameter) && + size(ModelingToolkit.Symbolics.wrap(parameter)) == (size(snapshot, 1), pod_dim) +end +@test length(reconstruction_parameters) == 1 -deim_sol = solve(deim_prob, Rodas5P(), saveat = 1.0) +deim_prob = DAEProblem( + deim_sys, nothing, full_prob.tspan; + build_initializeprob = false, +) + +deim_sol = solve(deim_prob; saveat = 1.0) +@test successful_retcode(deim_sol) nₓ = length(sol[x]) nₜ = length(sol[t]) -# test solution retrieva +# test solution retrieval @test size(deim_sol[v(x, t)]) == (nₓ, nₜ) @test size(deim_sol[w(x, t)]) == (nₓ, nₜ) +full_fields = vcat(sol[v(x, t)], sol[w(x, t)]) +reduced_fields = vcat(deim_sol[v(x, t)], deim_sol[w(x, t)]) +@test norm(reduced_fields - full_fields) / norm(full_fields) < 0.6 + +@parameters a = 2.0 +@variables z(t)[1:4] energy(t) +z_scalars = ModelingToolkit.Symbolics.value.( + ModelingToolkit.Symbolics.scalarize(z) +) +@named parameterized_array_system = System( + [Dt(z) ~ -z - a * (z .^ 3)], t, z_scalars, [a]; + observed = [energy ~ sum(z)], +) +parameterized_snapshot = [ + sin(0.3 * i * j) + cos(0.2 * i * (j + 1)) for i in 1:4, j in 1:8 +] +@test_throws DimensionMismatch deim( + parameterized_array_system, parameterized_snapshot, 2; snapshot_times = [0.0] +) +parameterized_reduced = deim(parameterized_array_system, parameterized_snapshot, 2) +parameterized_prob = DAEProblem( + parameterized_reduced, nothing, (0.0, 1.0); + build_initializeprob = false, +) +@test parameterized_prob.ps[a] == 2.0 +lift_parameter = only( + filter(ModelingToolkit.parameters(parameterized_reduced)) do p + is_symbolic_array(p) && size(ModelingToolkit.Symbolics.wrap(p)) == (4, 2) + end +) +lift = parameterized_prob.ps[lift_parameter] +nonlinear_basis = svd(-2.0 .* parameterized_snapshot .^ 3).U[:, 1:2] +first_index = argmax(abs.(nonlinear_basis[:, 1])) +second_residual = nonlinear_basis[:, 2] - nonlinear_basis[:, 1] * + (nonlinear_basis[first_index, 2] / nonlinear_basis[first_index, 1]) +indices = [first_index, argmax(abs.(second_residual))] +for state in ([0.1, -0.2], [-0.3, 0.4]) + full_state = lift * state + expected = -state + lift' * nonlinear_basis * + (nonlinear_basis[indices, :] \ (-2.0 .* full_state[indices] .^ 3)) + residual = zeros(2) + parameterized_prob.f(residual, expected, state, parameterized_prob.p, 0.0) + @test norm(residual) < 1.0e-10 +end + +parameterized_sol = solve(parameterized_prob; saveat = 0.2) +@test successful_retcode(parameterized_sol) +@test length(parameterized_sol[energy]) == 6 +completed_array_reduced = deim( + complete(deepcopy(parameterized_array_system)), parameterized_snapshot, 2 +) +@test length(ModelingToolkit.get_eqs(completed_array_reduced)) == 1 +@test is_symbolic_array(only(ModelingToolkit.get_eqs(completed_array_reduced)).lhs) +@test is_symbolic_array(only(ModelingToolkit.get_eqs(completed_array_reduced)).rhs) + +# `mtkcompile` scalarizes this array equation and reverses its registered element order. +# DEIM must still emit one correctly ordered array equation. +@variables compiled_z(t)[1:4] +@mtkcompile scalarized_array_system = System( + [Dt(compiled_z) ~ -compiled_z - compiled_z .^ 3], t; + name = :scalarized_array_system, +) +@test length(ModelingToolkit.get_eqs(scalarized_array_system)) == 4 +@test all(ModelingToolkit.get_eqs(scalarized_array_system)) do equation + !is_symbolic_array(equation.lhs) && !is_symbolic_array(equation.rhs) +end +scalarized_snapshot = [ + (0.2 + 0.1 * i) * exp(-0.1 * j) + (1.0 - 0.05 * i) * sin(0.2 * j) + for i in 1:4, j in 1:8 +] +array_codegen_reduced = deim(scalarized_array_system, scalarized_snapshot, 2) +array_codegen_equations = ModelingToolkit.get_eqs(array_codegen_reduced) +@test length(array_codegen_equations) == 1 +@test is_symbolic_array(only(array_codegen_equations).lhs) +@test is_symbolic_array(only(array_codegen_equations).rhs) +array_codegen_observed = ModelingToolkit.get_observed(array_codegen_reduced) +@test length(array_codegen_observed) == 1 +@test is_symbolic_array(only(array_codegen_observed).lhs) +@test is_symbolic_array(only(array_codegen_observed).rhs) +array_codegen_prob = DAEProblem( + array_codegen_reduced, nothing, (0.0, 1.0); + build_initializeprob = false, +) +array_codegen_sol = solve(array_codegen_prob; saveat = 0.2) +@test successful_retcode(array_codegen_sol) +@test size(Array(array_codegen_sol)) == (2, 6) +array_observed_function = SII.observed(array_codegen_reduced, compiled_z) +initial_reconstruction = array_observed_function( + array_codegen_prob.u0, array_codegen_prob.p, first(array_codegen_prob.tspan) +) +@test initial_reconstruction ≈ + reverse(scalarized_snapshot[:, 1]) + +# Multidimensional arrays may be scalarized into a noncontiguous permutation rather than a +# simple reversal. Reconstruction must restore their canonical Cartesian index order. +@variables compiled_matrix(t)[1:2, 1:3] +@mtkcompile scalarized_matrix_system = System( + [Dt(compiled_matrix) ~ -compiled_matrix - compiled_matrix .^ 3], t; + name = :scalarized_matrix_system, +) +matrix_unknowns = ModelingToolkit.get_unknowns(scalarized_matrix_system) +matrix_rows_by_index = Dict{Tuple, Int}() +for (row, unknown) in enumerate(matrix_unknowns) + index_arguments = ModelingToolkit.SymbolicUtils.arguments( + ModelingToolkit.Symbolics.unwrap(unknown) + ) + index = Tuple( + ModelingToolkit.SymbolicUtils.unwrap_const.(index_arguments[2:end]) + ) + matrix_rows_by_index[index] = row +end +matrix_canonical_rows = vec( + [ + matrix_rows_by_index[Tuple(index)] for index in CartesianIndices((2, 3)) + ] +) +@test length(unique(diff(matrix_canonical_rows))) > 1 +matrix_snapshot = [ + (0.3 + 0.07 * i) * exp(-0.1 * j) + (0.8 - 0.03 * i) * sin(0.2 * j) + for i in 1:6, j in 1:8 +] +matrix_reduced = deim(scalarized_matrix_system, matrix_snapshot, 2) +matrix_equations = ModelingToolkit.get_eqs(matrix_reduced) +@test length(matrix_equations) == 1 +@test is_symbolic_array(only(matrix_equations).lhs) +@test is_symbolic_array(only(matrix_equations).rhs) +matrix_prob = DAEProblem( + matrix_reduced, nothing, (0.0, 1.0); + build_initializeprob = false, +) +matrix_sol = solve(matrix_prob; saveat = 0.2) +@test successful_retcode(matrix_sol) +matrix_observed_function = SII.observed(matrix_reduced, compiled_matrix) +initial_matrix_reconstruction = matrix_observed_function( + matrix_prob.u0, matrix_prob.p, first(matrix_prob.tspan) +) +@test initial_matrix_reconstruction ≈ + reshape(matrix_snapshot[matrix_canonical_rows, 1], 2, 3) + +function tree_size(expression) + value = ModelingToolkit.Symbolics.unwrap(expression) + ModelingToolkit.SymbolicUtils.iscall(value) || return 1 + return 1 + sum( + tree_size, ModelingToolkit.SymbolicUtils.arguments(value); init = 0 + ) +end + +function reduced_system_size(intervals) + local_dx = L / intervals + local_discretization = MOLFiniteDifference( + [x => local_dx], t; approx_order = order + ) + local_sys, _ = symbolic_discretize(pde_sys, local_discretization) + rows = length(ModelingToolkit.get_unknowns(local_sys)) + local_snapshot = [ + sin(0.13 * i * j) + cos(0.07 * i * (j + 1)) for i in 1:rows, j in 1:8 + ] + local_reduced = deim(local_sys, local_snapshot, 2; snapshot_times = collect(0.0:0.1:0.7)) + local_equations = ModelingToolkit.get_eqs(local_reduced) + local_observed = ModelingToolkit.get_observed(local_reduced) + return ( + equations = length(local_equations), + equation_tree = sum(local_equations; init = 0) do equation + tree_size(equation.lhs) + tree_size(equation.rhs) + end, + observed = length(local_observed), + observed_tree = sum(local_observed; init = 0) do equation + tree_size(equation.lhs) + tree_size(equation.rhs) + end, + array_parameters = count( + is_symbolic_array, ModelingToolkit.get_ps(local_reduced) + ), + ) +end + +# Fixed POD/DEIM dimensions produce a grid-independent generated symbolic graph. +@test reduced_system_size(5) == reduced_system_size(20) + +@testset "training parameters and nonautonomous nonlinearities" begin + @parameters α = 2.0 + @variables q(t)[1:4] + q0 = [0.4, 0.3, 0.2, 0.1] + qstates = ModelingToolkit.Symbolics.unwrap.(ModelingToolkit.Symbolics.scalarize(q)) + @named training_system = System( + [Dt(q) ~ -q - α * (1 + sin(t)) * q .^ 3], t, qstates, [α]; + initial_conditions = [ + q => q0, Dt(q) => -q0 - 3 .* q0 .^ 3, + ] + ) + training_problem = DAEProblem( + complete(training_system), [α => 3.0], (0.0, 0.5); build_initializeprob = false + ) + training_solution = solve(training_problem; saveat = 0.05, abstol = 1.0e-10, reltol = 1.0e-10) + @test successful_retcode(training_solution) + trained = deim(training_problem, training_solution, 2) + trained_problem = DAEProblem(trained, nothing, training_problem.tspan; build_initializeprob = false) + @test trained_problem.ps[α] == 3.0 + trained_solution = solve(trained_problem; saveat = 0.05, abstol = 1.0e-10, reltol = 1.0e-10) + @test successful_retcode(trained_solution) + truth = Array(training_solution) + reconstruct = SII.observed(trained, q) + reconstruction = reduce( + hcat, [ + reconstruct(state, trained_problem.p, time) + for (state, time) in zip(trained_solution.u, trained_solution.t) + ] + ) + @test norm(reconstruction - truth) / norm(truth) < 0.01 +end