Skip to content

feat!: make array-form DAE discretization the v1 default - #650

Merged
ChrisRackauckas merged 13 commits into
SciML:masterfrom
ChrisRackauckas-Claude:v1
Aug 21, 2026
Merged

feat!: make array-form DAE discretization the v1 default#650
ChrisRackauckas merged 13 commits into
SciML:masterfrom
ChrisRackauckas-Claude:v1

Conversation

@ChrisRackauckas-Claude

@ChrisRackauckas-Claude ChrisRackauckas-Claude commented Aug 20, 2026

Copy link
Copy Markdown
Member

Ignore this draft until it has been reviewed by @ChrisRackauckas.

What changed and why

This is the combined MethodOfLines v1 change. Time-dependent PDE systems now discretize to an implicit DAEProblem by default, so array-form equations reach code generation without mtkcompile scalarizing them. Systems that cannot use the first-order DAE path fall back to the compiled ODEProblem path; fallback = false makes that fallback an error. Stationary systems continue to return NonlinearProblems, and every time-dependent solve path retains PDETimeSeriesSolution wrapping.

Discretization strategy selection is removed entirely. MethodOfLines always attempts the array representation and automatically emits pointwise equations for individual patterns that have no array form. The use_ODAE switch is also removed. Callers that need an ODEProblem now use the documented symbolic_discretize + mtkcompile + ODEProblem route explicitly. DAE examples call solve(prob) so OrdinaryDiffEq chooses its default DAE algorithm instead of MethodOfLines documentation hard-coding DFBDF.

The package version moves to 1.0.0. Documentation and tests were migrated according to the problem they need: default DAE solves for implicit examples, and explicit symbolic-to-ODE construction for SSPRK/WENO/staggered examples.

This supersedes the earlier array-discretization PR and subsumes the default flip requested by the older scalar-default PR. The older branch was not replayed because its stale diff would revert newer array-discretization work; its behavioral intent is implemented here by removing the selector entirely.

The branch is rebased through the periodic nonuniform WENO work linked below. That implementation and its seam-coordinate, accuracy, minimal-grid, and unsupported-pattern coverage are preserved under the unconditional array API.

Dependency boundary

This requires the registered DiffEqBase = "7.18.1" fix linked below. ModelingToolkitBase = "1.67" is the tested minimum: 1.66 carries DAEProblem.problem_type through construction but does not yet assemble array residuals, while 1.67 does. The SciMLBase floor supplies DAEProblem.problem_type support.

The newly direct runtime dependency ModelingToolkitBase is MIT licensed. No new solver dependency is introduced; the temporary direct OrdinaryDiffEqBDF test/docs dependency was removed once examples changed to solve(prob).

Verification

Focused runtime behavior on Julia 1.12.6:

problem_type=DAEProblem
algorithm_type=DFBDF  # selected by solve(prob) in the resolved local environment
retcode=Success
max_error=0.0007547373427254089
public_array_strategy=false
public_scalar_strategy=false
legacy_option_names_present=false
ode_problem_type=ODEProblem
ode_retcode=Success

The full rewritten array group passed:

Array discretization             | 3110 / 3110 passed | 27m39.3s
DAEProblem on the array form      |   49 /   49 passed |  1m39.9s
Testing MethodOfLines tests passed

After DiffEqBase 7.18.1 was registered, 68/69 CI checks passed. The remaining WENO
benchmark failed because its explicit SSPRK33 solve still consumed the new
DAEProblem returned by discretize:

Problem type: SciMLBase.DAEProblem
Solver type: OrdinaryDiffEqSSPRK.SSPRK33
Problem types compatible with the chosen solver: SciMLBase.ODEProblem

The benchmark now constructs its RHS/solve problem through symbolic_discretize,
mtkcompile, and ODEProblem, while its dedicated discretization benchmark continues
to exercise discretize. The identical failing leaf passed locally after the change:

run(SUITE["solve"]["burgers"]["stretched"]["N=128"]; samples=1, evals=1)
Trial(25.353 ms)

The resulting CI run passed all 69 checks, including benchmarks, documentation,
downgrade, QA, the complete test matrix, and PDESystemLibrary downstream.

The exact downgrade environment discriminated the ModelingToolkitBase floor with the
same unchanged GROUP=MOL_Interface2 command:

ModelingToolkitBase 1.65 | 3 passed, 1 failed, 2 errored, 8 broken
ModelingToolkitBase 1.66 | 4 passed, 1 failed, 0 errored, 8 broken
ModelingToolkitBase 1.67 | 5 passed, 0 failed, 0 errored, 8 broken | 9m56.0s
Testing MethodOfLines tests passed

Final-tree QA passed with the warmed absolute depot:

QA | 12 passed, 6 existing broken, 18 total | 13m26.0s
Testing MethodOfLines tests passed

The first QA attempt had one Aqua infrastructure failure: its cold, relative-depot persistent-task subprocess exited without creating done.log; the other 11 checks passed. No test or timeout setting changed for the successful retry.

Static checks:

julia +release -m Runic --check .                                            # exit 0
git diff --name-only --diff-filter=ACMRT -z upstream/master...HEAD |
  xargs -0 typos --force-exclude                                              # exit 0
git diff --check upstream/master...HEAD                                       # exit 0

The added source/test comment share is 5.11% (49 comment lines among 958 added lines). The explicit-ODE helper is defined once in test/shared/ode_discretize.jl; no duplicated test bodies or new public definitions were added.

With a local checkout of the required DiffEqBase fix, the unchanged strict docs build completed every example, cross-reference, missing-docstring check, and document check. It then failed external link checking when GitHub returned HTTP 429 for existing source links. Repeated unchanged attempts reached the same external limit; link checking was not disabled or weakened. The initial CI run predated the DiffEqBase 7.18.1 registration and failed at dependency resolution. After registration, a fresh detached checkout resolved DiffEqBase 7.18.1 from General; commit 00788307 retriggered the affected workflows without changing the source tree and commit 972d2586 fixes the sole resulting benchmark failure. All 69 checks on 972d2586 pass.

GROUP=Everything julia +release --project=. -e 'using Pkg; Pkg.test()' ran the
runtime groups on parent commit 545c3231. Representative results were Array
3111/3111, DAE array 49/49, and MOL_Interface2 5 passed / 8 existing broken. It
then failed, before DiffEqBase 7.18.1 was registered, when SciMLTesting created an
isolated QA environment that could only resolve DiffEqBase 7.17.1. The current source
difference is removal of the two legacy-keyword guards; the current tree is covered by
the focused, exact-downgrade, and QA runs above.

The PDESystemLibrary downstream job reaches dependency resolution and fails because its
current MethodOfLines = "0.12" compat excludes 1.0.0. A focused downstream compat and
documented symbolic-to-ODE migration passed 427/427 examples against MethodOfLines 0.12.
Against this v1 branch, the unchanged pre-v1 call shape errored 427/427 times on the
intentional DAEProblem/ODE-solver incompatibility; the migrated v1 run passed 427/427 in
12m19.1s. All nine checks on the supporting downstream draft pass, including Core on
Julia LTS, current, and prerelease, QA, downgrade, Runic, and spelling. The draft is
linked below.

Not run locally: Julia prerelease and GPU jobs.

Reviewer decisions

  • The default return-type and public API removals are intentionally breaking and paired with the 1.0.0 bump.
  • There is no public or internal whole-system strategy selector. Pointwise generation remains only as the automatic per-equation fallback used by the array implementation.
  • There is no use_ODAE compatibility switch. Explicit ODE integrators require the documented symbolic_discretize route.
  • DAE examples rely on solve(prob) rather than promising a specific OrdinaryDiffEq default algorithm.
  • The compiled ODE fallback can be disabled with fallback = false when implicit DAE construction is required.

Links

ChrisRackauckas and others added 9 commits August 20, 2026 19:47
`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>
The array-equation implicit-DAE codegen landed in ModelingToolkitBase 1.67.0,
released as ModelingToolkit 11.39.1. The `DAEProblem` path here does not work
against earlier versions, which reject array equations outright.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
`discretize` now builds a `DAEProblem` for time-dependent systems. The residuals
MethodOfLines emits are already implicit-DAE form, so no `mtkcompile` is needed
and the array equations survive into the generated code; isolating the
derivative for an `ODEProblem` is what scalarizes them. Systems that cannot be
posed as a first-order DAE fall back to `mtkcompile` plus an `ODEProblem`,
which `fallback = false` turns into an error.

`ScalarizedDiscretization` is removed. `ArrayDiscretization` is the default and
only strategy, with the scalar form as its fallback for patterns that have no
slice representation.

`symbolic_discretize` is documented as the entry point for building other
problem types.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
…ymbolic_discretize

Tutorials solving with an implicit method now solve the `DAEProblem` that
`discretize` returns, with `DFBDF()`. The hyperbolic showcases use explicit
solvers, which need an `ODEProblem`, so those build one from
`symbolic_discretize` and `mtkcompile` and say why.

Tests exercise the pointwise form through the internal
`PointwiseDiscretization`, which is what `ArrayDiscretization` falls back to,
so the array-vs-pointwise equivalence checks are preserved.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Complete the solver-path migration across the functional test groups and documentation, keep compiled ODE coverage where explicit solvers or exact interface parity require it, and use the owning package's public initialization-equation API.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Remove the public strategy selector and the use_ODAE compatibility switch. Array-form equations are now unconditional with automatic per-equation pointwise fallback, while callers that need an ODEProblem use symbolic_discretize and mtkcompile explicitly. DAE examples and tests defer to the default solver.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Remove the remaining compatibility guards and their tests so the retired discretization_strategy and use_ODAE names are absent entirely. Require the DiffEqBase release containing the DAE initialization despecialization fix.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Require the first ModelingToolkitBase release that assembles array residuals for implicit DAE code generation.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
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