Skip to content

Migrate ABDF2 eulercache to ESDIRKIMEXCache, drop legacy IE caches#3656

Merged
ChrisRackauckas merged 2 commits into
SciML:masterfrom
singhharsh1708:bdf-migrate-abdf2-to-esdirkimex
May 20, 2026
Merged

Migrate ABDF2 eulercache to ESDIRKIMEXCache, drop legacy IE caches#3656
ChrisRackauckas merged 2 commits into
SciML:masterfrom
singhharsh1708:bdf-migrate-abdf2-to-esdirkimex

Conversation

@singhharsh1708

@singhharsh1708 singhharsh1708 commented May 20, 2026

Copy link
Copy Markdown
Contributor

Summary

Migrates ABDF2’s first-step ImplicitEuler fallback to the unified ESDIRKIMEXCache path and removes the legacy ImplicitEulerCache / ImplicitEulerConstantCache from OrdinaryDiffEqSDIRK.

This was the last remaining non-generic SDIRK cache path.

Verification

  • ABDF2 matches master bit-exactly (naccept, u_end) for both IIP/OOP
  • SDIRK convergence: 82 pass / 2 broken
  • DAE suite: 8/8 pass
  • Existing FBDF/DFBDF allocation failures are unchanged and pre-existing

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

@ChrisRackauckas

Copy link
Copy Markdown
Member

Needs a new minimum on the Project.toml.

@singhharsh1708 singhharsh1708 force-pushed the bdf-migrate-abdf2-to-esdirkimex branch from 78c915d to 414e1ea Compare May 20, 2026 13:54
@singhharsh1708

Copy link
Copy Markdown
Contributor Author

@ChrisRackauckas done..

@ChrisRackauckas ChrisRackauckas merged commit 62e3886 into SciML:master May 20, 2026
93 of 102 checks passed
@singhharsh1708 singhharsh1708 deleted the bdf-migrate-abdf2-to-esdirkimex branch May 20, 2026 15:54
ChrisRackauckas-Claude pushed a commit to ChrisRackauckas-Claude/OrdinaryDiffEq.jl that referenced this pull request May 27, 2026
SciML#3656 retired the ABDF2-specific `ImplicitEulerCache::perform_step!`
and routed the first-step Euler bootstrap through the new generic
`ESDIRKIMEXCache::perform_step!` in
`lib/OrdinaryDiffEqSDIRK/src/generic_imex_perform_step.jl`. That path's
two stage-1 nlsolve-initial-guess branches (the `Predictor.MaxOrder`
current-extrapolant branch and the `Predictor.Linear` `dt·fsalfirst`
branch) were gated on:

    alg isa Union{
        OrdinaryDiffEqNewtonAdaptiveSDIRKAlgorithm,
        OrdinaryDiffEqNewtonNonAdaptiveSDIRKAlgorithm,
        ImplicitEuler, Trapezoid,
    }

That alg-allow-list is misaligned with what the code actually wants to
discriminate: the IE tableau is the only one in this file with an
implicit first stage (`explicit_first_stage = false`) and
`stage1_extrapolation = true` (Trapezoid has `explicit_first_stage =
true` and skips this branch; the SFSDIRK / Cash4 tableaus set
`stage1_extrapolation = false` and also opt out). The alg list was an
indirect way of saying "this stage 1 is an IE step", which is more
precisely expressed as `E === :ie_dd2` on the tableau.

Replacing the alg-Union with `E === :ie_dd2` also covers the
previously-missed ABDF2-bootstrap path. The ABDF2 cache holds a
`cache.eulercache::ESDIRKIMEXCache` built from
`ImplicitEulerESDIRKIMEXTableau` and dispatches its first step through
the same `_perform_step_iip!` / `_perform_step_oop!`. With
`integrator.alg isa ABDF2` (not in the old allow-list), the stage-1
nlsolve initial guess was falling through to `zs[1] = 0` (or
`z1 = zero(u)`) instead of the linear extrapolant. With
`NonlinearSolveAlg(TrustRegion())` (used by the `nlstep_tests.jl`
convergence checks), the loose `iter == 1 && ndz < 1.0e-5` early-exit
in `nlsolve!` then terminated after one Newton step at a residual
that, at small `dt`, dominated the truncation error — collapsing
`ABDF2`'s measured convergence order from ~2 to ~1.534 on
`nlstep_tests.jl:52` and `:108`.

Bisected to SciML#3656 (`62e3886ff`). The pre-SciML#3656 code path was the
explicit `ImplicitEulerCache::perform_step!`, which unconditionally
seeded `nlsolver.z = dt * integrator.fsalfirst`; this commit restores
that behavior in a way scoped to the IE-bootstrap tableau and shared
between IE-as-algorithm and BDF callers.

Verification on Julia 1.11 with this patch:

- `nlstep_tests.jl:52` (autonomous Robertson, `nlstep=true`): errors
  back to `[6.0e-7, 1.5e-7, 3.6e-8, 8.7e-9]`, `𝒪est[:l∞] = 2.039`.
- `nlstep_tests.jl:108` (non-autonomous Robertson, `nlstep=true`):
  same result.
- Full `nlstep_tests.jl` suite on master + this fix: 19 pass / 0 fail /
  6 broken (the `@test_broken` were pre-existing TRBDF2/KenCarp4 issues
  on the 1D-reduced `nlstep` path, orthogonal to this).
- Smoke test of `ImplicitEuler`, `SDIRK2`, `Cash4`, `Hairer4`,
  `Hairer42`, `ImplicitMidpoint`, `Trapezoid`, `ABDF2` on a stiff
  exponential decay (`du = -50u`, `dt = 0.01`): all converge to
  `exp(-50) ≈ 1.93e-22` cleanly; no regression from dropping the
  alg-allow-list.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
ChrisRackauckas-Claude pushed a commit to ChrisRackauckas-Claude/OrdinaryDiffEq.jl that referenced this pull request May 27, 2026
SciML#3656 retired the ABDF2-specific `ImplicitEulerCache::perform_step!`
and routed the first-step Euler bootstrap through the new generic
`ESDIRKIMEXCache::perform_step!` in
`lib/OrdinaryDiffEqSDIRK/src/generic_imex_perform_step.jl`. That path's
stage-1 nlsolve-initial-guess code had two branches:

- a `Predictor.MaxOrder` branch (use the previous step's interpolant)
- a `Predictor.Linear` branch (use `dt·fsalfirst`)

both gated on an alg-allow-list:

    alg isa Union{
        OrdinaryDiffEqNewtonAdaptiveSDIRKAlgorithm,
        OrdinaryDiffEqNewtonNonAdaptiveSDIRKAlgorithm,
        ImplicitEuler, Trapezoid,
    }

That alg-allow-list was misaligned with what the code actually wants to
discriminate: in this `else` branch (implicit first stage), the IE
tableau is the only configuration that reaches this code with
`stage1_extrapolation = true` — Trapezoid and the other Newton-SDIRKs in
the Union all have `explicit_first_stage = true` and skip into the `if
tab.explicit_first_stage` arm above. The alg-Union was an indirect way
of saying "this stage 1 is an IE step", which is more precisely
expressed as `E === :ie_dd2` on the tableau.

This commit replaces the alg-Union with `E === :ie_dd2` and also drops
the `Predictor.MaxOrder` branch — it was dead code given that the only
algorithm reaching it is `ImplicitEuler`, whose interpolant *is* linear,
so the MaxOrder predictor produces essentially the same guess as the
default Linear predictor. Net result: -27 / +5 (smaller diff than my
first attempt, which kept both branches).

The replacement covers the previously-missed ABDF2-bootstrap path. The
ABDF2 cache holds a `cache.eulercache::ESDIRKIMEXCache` built from
`ImplicitEulerESDIRKIMEXTableau` and dispatches its first step through
the same `_perform_step_iip!` / `_perform_step_oop!`. With
`integrator.alg isa ABDF2` (not in the old allow-list), the stage-1
nlsolve initial guess was falling through to `zs[1] = 0` (or
`z1 = zero(u)`) instead of the linear extrapolant. With
`NonlinearSolveAlg(TrustRegion())` (used by the `nlstep_tests.jl`
convergence checks), the loose `iter == 1 && ndz < 1.0e-5` early-exit
in `nlsolve!` then terminated after one Newton step at a residual
that, at small `dt`, dominated the truncation error — collapsing
`ABDF2`'s measured convergence order from ~2 to ~1.534 on
`nlstep_tests.jl:52` and `:108`.

Bisected to SciML#3656 (`62e3886ff`). The pre-SciML#3656 code path was the
explicit `ImplicitEulerCache::perform_step!`, which unconditionally
seeded `nlsolver.z = dt * integrator.fsalfirst`; this commit restores
that behavior in a way scoped to the IE-bootstrap tableau and shared
between IE-as-algorithm and BDF callers.

Verification on Julia 1.11 with this patch:

- `nlstep_tests.jl:52` (autonomous Robertson, `nlstep=true`): errors
  back to `[6.0e-7, 1.5e-7, 3.6e-8, 8.7e-9]`, `𝒪est[:l∞] = 2.039`.
- `nlstep_tests.jl:108` (non-autonomous Robertson, `nlstep=true`):
  same result.
- Full `nlstep_tests.jl` suite on master + this fix: 19 pass / 0 fail /
  6 broken (the `@test_broken` were pre-existing TRBDF2/KenCarp4 issues
  on the 1D-reduced `nlstep` path, orthogonal to this).
- Smoke test of `ImplicitEuler`, `SDIRK2`, `Cash4`, `Hairer4`,
  `Hairer42`, `ImplicitMidpoint`, `Trapezoid`, `ABDF2` on a stiff
  exponential decay (`du = -50u`, `dt = 0.01`): all converge to
  `exp(-50) ≈ 1.93e-22` cleanly; no regression from dropping the
  alg-allow-list or the MaxOrder branch.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
ChrisRackauckas-Claude pushed a commit to ChrisRackauckas-Claude/OrdinaryDiffEq.jl that referenced this pull request May 27, 2026
SciML#3656 retired the ABDF2-specific `ImplicitEulerCache::perform_step!`
and routed the first-step Euler bootstrap through the new generic
`ESDIRKIMEXCache::perform_step!` in
`lib/OrdinaryDiffEqSDIRK/src/generic_imex_perform_step.jl`. That path's
stage-1 nlsolve-initial-guess code had two branches:

- a `Predictor.MaxOrder` branch (use the previous step's interpolant)
- a `Predictor.Linear` branch (use `dt·fsalfirst`)

both gated on an alg-allow-list:

    alg isa Union{
        OrdinaryDiffEqNewtonAdaptiveSDIRKAlgorithm,
        OrdinaryDiffEqNewtonNonAdaptiveSDIRKAlgorithm,
        ImplicitEuler, Trapezoid,
    }

That alg-allow-list was misaligned with what the code actually wants to
discriminate: in this `else` branch (implicit first stage), the IE
tableau is the only configuration that reaches this code with
`stage1_extrapolation = true` — Trapezoid and the other Newton-SDIRKs in
the Union all have `explicit_first_stage = true` and skip into the `if
tab.explicit_first_stage` arm above. The alg-Union was an indirect way
of saying "this stage 1 is an IE step", which is more precisely
expressed as `E === :ie_dd2` on the tableau.

This commit replaces the alg-Union with `E === :ie_dd2` and also drops
the `Predictor.MaxOrder` branch — it was dead code given that the only
algorithm reaching it is `ImplicitEuler`, whose interpolant *is* linear,
so the MaxOrder predictor produces essentially the same guess as the
default Linear predictor. Net result: -27 / +5 (smaller diff than my
first attempt, which kept both branches).

The replacement covers the previously-missed ABDF2-bootstrap path. The
ABDF2 cache holds a `cache.eulercache::ESDIRKIMEXCache` built from
`ImplicitEulerESDIRKIMEXTableau` and dispatches its first step through
the same `_perform_step_iip!` / `_perform_step_oop!`. With
`integrator.alg isa ABDF2` (not in the old allow-list), the stage-1
nlsolve initial guess was falling through to `zs[1] = 0` (or
`z1 = zero(u)`) instead of the linear extrapolant. With
`NonlinearSolveAlg(TrustRegion())` (used by the `nlstep_tests.jl`
convergence checks), the loose `iter == 1 && ndz < 1.0e-5` early-exit
in `nlsolve!` then terminated after one Newton step at a residual
that, at small `dt`, dominated the truncation error — collapsing
`ABDF2`'s measured convergence order from ~2 to ~1.534 on
`nlstep_tests.jl:52` and `:108`.

Bisected to SciML#3656 (`62e3886ff`). The pre-SciML#3656 code path was the
explicit `ImplicitEulerCache::perform_step!`, which unconditionally
seeded `nlsolver.z = dt * integrator.fsalfirst`; this commit restores
that behavior in a way scoped to the IE-bootstrap tableau and shared
between IE-as-algorithm and BDF callers.

Verification on Julia 1.11 with this patch:

- `nlstep_tests.jl:52` (autonomous Robertson, `nlstep=true`): errors
  back to `[6.0e-7, 1.5e-7, 3.6e-8, 8.7e-9]`, `𝒪est[:l∞] = 2.039`.
- `nlstep_tests.jl:108` (non-autonomous Robertson, `nlstep=true`):
  same result.
- Full `nlstep_tests.jl` suite on master + this fix: 19 pass / 0 fail /
  6 broken (the `@test_broken` were pre-existing TRBDF2/KenCarp4 issues
  on the 1D-reduced `nlstep` path, orthogonal to this).
- Smoke test of `ImplicitEuler`, `SDIRK2`, `Cash4`, `Hairer4`,
  `Hairer42`, `ImplicitMidpoint`, `Trapezoid`, `ABDF2` on a stiff
  exponential decay (`du = -50u`, `dt = 0.01`): all converge to
  `exp(-50) ≈ 1.93e-22` cleanly; no regression from dropping the
  alg-allow-list or the MaxOrder branch.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
ChrisRackauckas added a commit that referenced this pull request May 28, 2026
#3694)

#3656 retired the ABDF2-specific `ImplicitEulerCache::perform_step!`
and routed the first-step Euler bootstrap through the new generic
`ESDIRKIMEXCache::perform_step!` in
`lib/OrdinaryDiffEqSDIRK/src/generic_imex_perform_step.jl`. That path's
stage-1 nlsolve-initial-guess code had two branches:

- a `Predictor.MaxOrder` branch (use the previous step's interpolant)
- a `Predictor.Linear` branch (use `dt·fsalfirst`)

both gated on an alg-allow-list:

    alg isa Union{
        OrdinaryDiffEqNewtonAdaptiveSDIRKAlgorithm,
        OrdinaryDiffEqNewtonNonAdaptiveSDIRKAlgorithm,
        ImplicitEuler, Trapezoid,
    }

That alg-allow-list was misaligned with what the code actually wants to
discriminate: in this `else` branch (implicit first stage), the IE
tableau is the only configuration that reaches this code with
`stage1_extrapolation = true` — Trapezoid and the other Newton-SDIRKs in
the Union all have `explicit_first_stage = true` and skip into the `if
tab.explicit_first_stage` arm above. The alg-Union was an indirect way
of saying "this stage 1 is an IE step", which is more precisely
expressed as `E === :ie_dd2` on the tableau.

This commit replaces the alg-Union with `E === :ie_dd2` and also drops
the `Predictor.MaxOrder` branch — it was dead code given that the only
algorithm reaching it is `ImplicitEuler`, whose interpolant *is* linear,
so the MaxOrder predictor produces essentially the same guess as the
default Linear predictor. Net result: -27 / +5 (smaller diff than my
first attempt, which kept both branches).

The replacement covers the previously-missed ABDF2-bootstrap path. The
ABDF2 cache holds a `cache.eulercache::ESDIRKIMEXCache` built from
`ImplicitEulerESDIRKIMEXTableau` and dispatches its first step through
the same `_perform_step_iip!` / `_perform_step_oop!`. With
`integrator.alg isa ABDF2` (not in the old allow-list), the stage-1
nlsolve initial guess was falling through to `zs[1] = 0` (or
`z1 = zero(u)`) instead of the linear extrapolant. With
`NonlinearSolveAlg(TrustRegion())` (used by the `nlstep_tests.jl`
convergence checks), the loose `iter == 1 && ndz < 1.0e-5` early-exit
in `nlsolve!` then terminated after one Newton step at a residual
that, at small `dt`, dominated the truncation error — collapsing
`ABDF2`'s measured convergence order from ~2 to ~1.534 on
`nlstep_tests.jl:52` and `:108`.

Bisected to #3656 (`62e3886ff`). The pre-#3656 code path was the
explicit `ImplicitEulerCache::perform_step!`, which unconditionally
seeded `nlsolver.z = dt * integrator.fsalfirst`; this commit restores
that behavior in a way scoped to the IE-bootstrap tableau and shared
between IE-as-algorithm and BDF callers.

Verification on Julia 1.11 with this patch:

- `nlstep_tests.jl:52` (autonomous Robertson, `nlstep=true`): errors
  back to `[6.0e-7, 1.5e-7, 3.6e-8, 8.7e-9]`, `𝒪est[:l∞] = 2.039`.
- `nlstep_tests.jl:108` (non-autonomous Robertson, `nlstep=true`):
  same result.
- Full `nlstep_tests.jl` suite on master + this fix: 19 pass / 0 fail /
  6 broken (the `@test_broken` were pre-existing TRBDF2/KenCarp4 issues
  on the 1D-reduced `nlstep` path, orthogonal to this).
- Smoke test of `ImplicitEuler`, `SDIRK2`, `Cash4`, `Hairer4`,
  `Hairer42`, `ImplicitMidpoint`, `Trapezoid`, `ABDF2` on a stiff
  exponential decay (`du = -50u`, `dt = 0.01`): all converge to
  `exp(-50) ≈ 1.93e-22` cleanly; no regression from dropping the
  alg-allow-list or the MaxOrder branch.

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.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