Skip to content

feat: DAEProblem path that keeps the array (slice-form) discretization - #639

Draft
ChrisRackauckas-Claude wants to merge 4 commits into
SciML:masterfrom
ChrisRackauckas-Claude:array-daeproblem
Draft

feat: DAEProblem path that keeps the array (slice-form) discretization#639
ChrisRackauckas-Claude wants to merge 4 commits into
SciML:masterfrom
ChrisRackauckas-Claude:array-daeproblem

Conversation

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member

What and why

discretize runs mtkcompile, which scalarizes array equations: an ODEProblem needs D(x) = f(x), and isolating the derivative is structural simplification. So the array (slice-form) discretization from #619 is flattened before it ever reaches codegen, and the scaling win is lost at exactly the point it would pay off.

The residuals MethodOfLines already emits, D(u) - f ~ 0, are in implicit-DAE form. This adds a DAEProblem(pdesys, disc) path that skips mtkcompile entirely, so the array equations survive into the generated code:

disc = MOLFiniteDifference([x => n], t; discretization_strategy = ArrayDiscretization())
prob = DAEProblem(pdesys, disc)
sol = solve(prob, DFBDF())
sol[u(t, x)]

Part of #428.

Consistent initialization

A fully implicit DAE needs consistent (u0, du0). initializealg defaults to BrownFullBasicInit(), but only when it is safe to do so: MethodOfLines always emits initialization equations (9 for the heat equation, 18 for the wave equation, 9 for the PDAE), so an isempty guard would be useless. Instead brown_init_offenders inspects what those equations actually constrain, and construction raises BrownFullBasicInitUnsafeError naming the offenders and pointing back at discretize when the algorithm would not honor them. Passing initializealg explicitly overrides both the default and the check. Systems second order in time need the order reduction mtkcompile performs and are rejected outright.

Solution wrapping

The solution is a PDETimeSeriesSolution, the same wrapper discretize produces, so it is indexed and interpolated by the PDESystem's own variables (sol[u(t, x)], sol(t, x)) rather than by discretized variables. This needed problem_type on DAEProblem, which did not exist — ODEProblem, DDEProblem, BVProblem and NonlinearProblem all had one and DAEProblem was the odd one out.

Dependencies — this cannot go green yet

Needed Status
SciML/SciMLBase.jl#1536problem_type field on DAEProblem merged, released as v3.48.0
SciML/ModelingToolkit.jl#4983 — array equations on the implicit-DAE codegen path open
SciML/ModelingToolkit.jl#4984 — forward a system's ProblemTypeCtx into DAEProblem open

SciMLBase compat is bumped to 3.48. The ModelingToolkit compat still needs a bump once #4983 and #4984 land in a release; until then CI here will fail, and I have deliberately not guessed a version number. Please merge those two first.

Verification

New test group test/Array_Discretization/dae_problem.jl (in the existing Array_Discretization group), run against a local environment with all three upstream changes applied:

1D heat, Dirichlet BCs                                |   10     10  2m36.3s
1D advection                                          |    3      3  15.2s
2D diffusion                                          |    3      3  56.4s
Brusselator, periodic BCs                             |    3      3  1m01.6s
PDAE with an algebraic variable                       |    3      3  38.7s
wave equation with a Dt initial condition is rejected |    9      9  3.7s
user-supplied initializealg wins                      |    2      2  13.3s
ScalarizedDiscretization is rejected                  |    2      2  0.1s
safety predicate branches                             |    7      7  0.6s

42 assertions, 0 failures. Every solving case is checked against the discretize (ODEProblem + mtkcompile) path it must reproduce, not just against an analytic solution — dae_vals ≈ ode_vals rtol = 1e-6 — and the 1D heat case additionally asserts the array equations actually survived (any(isarrayeq, get_eqs(prob.f.sys))), which is the whole point of the path.

End-to-end sanity on the heat equation, array form, no mtkcompile:

problem_type(prob) = MOLMetadata
typeof(sol) = PDETimeSeriesSolution
sol[u(t,x)] size = (2, 21)
max|u-exact| = 0.0007565021749123546
interp at (0.1, 0.5) = 0.3734643410283503  exact = 0.37270783885343794
retcode = Success

runic --check and typos are clean on the diff. Comment lines are 6.3% of added lines (36/567).

Not verified

  • The rest of the MethodOfLines suite has not been run against this branch. The environment needed to exercise it is a hand-assembled setup with two unreleased ModelingToolkit PRs deved; a normal Pkg.test() cannot resolve today. Only the Array_Discretization DAE group above was run. Once #4983/#4984 release, the full suite should be run before this merges.
  • No GPU or downstream lanes.

Worth pushing back on

  • test/qa/qa.jl gains an initialization_equations entry in the explicit-imports ignore list, following the existing pattern for names owned by ModelingToolkitBase and re-exported by ModelingToolkit.
  • brown_init_offenders is a heuristic over initialization equations. It is deliberately conservative — it errors rather than silently producing a wrong initialization — but it is the part of this PR most likely to need adjusting as more PDE shapes hit it.

Please ignore until reviewed by @ChrisRackauckas.

ChrisRackauckas and others added 4 commits August 16, 2026 00:59
`discretize` runs `mtkcompile`, which scalarizes the array equations
`ArrayDiscretization` emits: an `ODEProblem` needs `D(x) = f(x)`, and isolating
the derivative is structural simplification. The residuals MethodOfLines already
emits, `D(u) - f ~ 0`, are the implicit-DAE form, so `DAEProblem(pdesys, disc)`
builds a problem without `mtkcompile` and the array equations reach codegen.

`initializealg` defaults to `BrownFullBasicInit()`, the only algorithm that
reproduces the `discretize` result here. Because it takes the differential
variables' values as given and solves for everything else, it is chosen only when
every initialization equation fixes a single differential unknown to a value
involving no other unknown; otherwise construction raises an error naming the
offending equations and pointing at `discretize`. A user-supplied `initializealg`
overrides both. `ScalarizedDiscretization` and systems that are second order in
time are rejected with their own messages.

Requires the ModelingToolkit array-equation DAE fixes; see the PR body.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
PDEBase only emits an initialization equation for a variable whose time-derivative
order in the system exceeds the order of its initial condition, which makes that
variable differential by construction. The algebraic-variable and coupled-unknown
branches of the guard are therefore never exercised by a MethodOfLines-discretized
system; check them on hand-built systems so the guard is known to discriminate.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
With SciMLBase's `problem_type` field on `DAEProblem` and ModelingToolkit
forwarding the system's `ProblemTypeCtx` metadata into it, `wrap_sol` now
reaches the DAE path, so its solutions are indexed and interpolated by the
`PDESystem`'s variables like the `discretize` path.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
The wrapping commit updated the manual but left the docstring claiming the
result is a plain `DAEProblem` indexed by discretized variables. Wrapping needs
`problem_type` on `DAEProblem`, so bump the SciMLBase compat to match.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants