Skip to content

WIP: Initial IMEX tableau abstraction with generic perform_step! (ARS…#3196

Open
singhharsh1708 wants to merge 12 commits intoSciML:masterfrom
singhharsh1708:imex-tableau-prototype
Open

WIP: Initial IMEX tableau abstraction with generic perform_step! (ARS…#3196
singhharsh1708 wants to merge 12 commits intoSciML:masterfrom
singhharsh1708:imex-tableau-prototype

Conversation

@singhharsh1708
Copy link
Copy Markdown

@singhharsh1708 singhharsh1708 commented Mar 19, 2026

Summary

This PR introduces a unified IMEX tableau-based framework and migrates the KenCarp and Kvaerno IMEX SDIRK families to use it.

The goal is to define IMEX methods entirely via tableau coefficients, removing the need for per-method perform_step! implementations and caches.


Included in this PR

  • IMEXTableau abstraction (explicit + implicit Butcher arrays)

  • Generic perform_step! implementation for IMEX methods

  • Migration of:

    • KenCarp3, KenCarp4, KenCarp5, KenCarp47, KenCarp58
    • Kvaerno3, Kvaerno4, Kvaerno5
  • Removal of all per-method caches and perform_step! implementations for these methods (~2000 LOC reduction)

  • ARS343 retained as a reference implementation


Verification

  • Convergence order verified for all migrated methods
  • Works for SplitODEProblem (both IIP and OOP)
  • Adaptive stepping verified
  • All existing SDIRK tests pass
  • No regressions observed in existing test suite

Notes

  • CFNLIRK3 is intentionally left unchanged as it does not fit the IMEX tableau structure (non-adaptive, different abstraction)
  • This completes the migration of IMEX SDIRK methods to a unified tableau-based implementation
  • Future work may extend this approach to additional IMEX families

Checklist

  • Appropriate tests were added
  • No public API breakage
  • Documentation updated where needed
  • Code follows SciML Style Guide and COLPRAC
  • All new documentation uses public API

@singhharsh1708
Copy link
Copy Markdown
Author

singhharsh1708 commented Mar 20, 2026

@ChrisRackauckas This PR introduces a prototype IMEX tableau abstraction and generic perform_step! dispatch.

  • Fixed LTS precompilation issues (removed new algorithms from umbrella exports)
  • Verified SDIRK test suite passes locally
  • Float32 support confirmed
    Remaining CI failures are expected as this is a WIP prototype focused on architecture and design.

@singhharsh1708 singhharsh1708 force-pushed the imex-tableau-prototype branch from ee82a23 to 86409c9 Compare March 20, 2026 19:40
@singhharsh1708
Copy link
Copy Markdown
Author

singhharsh1708 commented Mar 22, 2026

@ChrisRackauckas what about this is this working fine and has correct structure?

@ChrisRackauckas
Copy link
Copy Markdown
Member

It looks like a good start but it needs all of the details of the other SDIRK methods

@singhharsh1708
Copy link
Copy Markdown
Author

singhharsh1708 commented Mar 23, 2026

@ChrisRackauckas Migrated the full KenCarp/Kvaerno family to IMEXTableau and removed the per-method implementations. All tests are passing.

)
end

abstract type OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ} <:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

top level should be in core.

end

function alg_cache(
alg::OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits},
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this doesn't cover all imex

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just sdirk

Comment on lines +230 to +232
if integrator.f isa SplitFunction && !repeat_step && !integrator.last_stepfail
f_impl(zs[1], integrator.uprev, p, integrator.t)
zs[1] .*= dt
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only ESDIRK?

@@ -0,0 +1,1504 @@
struct IMEXTableau{T, T2}
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ESDIRKIMEX?

nlsolver.z = z₂

@.. broadcast = false tmp = uprev + γ * z₁
@.. broadcast = false tmp = uprev
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also while updating these, please replace all uses of @.. broadcast=false with @.. (brodcast=false is the default)

SFSDIRK4, SFSDIRK5, CFNLIRK3, SFSDIRK6,
SFSDIRK7, SFSDIRK8, ESDIRK436L2SA2, ESDIRK437L2SA, ESDIRK547L2SA2, ESDIRK659L2SA
SFSDIRK7, SFSDIRK8, ESDIRK436L2SA2, ESDIRK437L2SA, ESDIRK547L2SA2, ESDIRK659L2SA,
OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm, ARS343
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

abstract types shouldn't export here.

Project.toml Outdated
CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not needed

Project.toml Outdated
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not needed

@ChrisRackauckas
Copy link
Copy Markdown
Member

Looks generally good. Needs a lot more to get the rest of the IMEX methods, but it's probably good to merge a version with just these ESDIRKs supported. Left a bunch of comments.

@singhharsh1708
Copy link
Copy Markdown
Author

singhharsh1708 commented Mar 24, 2026

@ChrisRackauckas Addressed all review comments:

  • Scoped IMEX implementation to ESDIRK methods (renamed types accordingly)
  • Cleaned up abstract type usage and exports
  • Fixed perform_step artifacts
  • Removed unnecessary dependencies
  • Removed all remaining @.. broadcast=false usages in the ESDIRK IMEX code paths.
    All tests passing and formatting clean. Let me know if anything else should be adjusted.
    Removed all remaining @.. broadcast=false usages in the ESDIRK IMEX code paths.

@singhharsh1708 singhharsh1708 force-pushed the imex-tableau-prototype branch 2 times, most recently from 7e4ce55 to 333d24b Compare March 25, 2026 07:20
@singhharsh1708
Copy link
Copy Markdown
Author

@ChrisRackauckas
Final check: all changes have been tested and are working as expected.

  • Migrated KenCarp and Kvaerno methods to the generic ESDIRK IMEX tableau
  • Verified OOP, IIP, and SplitODE behavior
  • All migrated methods pass tests with expected accuracy
  • No regressions observed in non-migrated methods
  • Runic formatting clean
    Additionally, corrected the ARS343 embedded error estimator to satisfy order conditions, improving adaptive stepping behavior.
    Please let me know if anything else should be addressed.

@ChrisRackauckas
Copy link
Copy Markdown
Member

Additionally, corrected the ARS343 embedded error estimator to satisfy order conditions, improving adaptive stepping behavior.

Can you explain that a bit more?

@singhharsh1708
Copy link
Copy Markdown
Author

earlier ARS343 was using a placeholder embedded method for error estimation:
b̃ = b - e₄ (i.e., using only the last stage)
While this keeps Σ b̃ = 0, it does not satisfy the second order condition
Σ b̃ᵢ cᵢ = 0, so it does not form a valid lower-order embedded pair.
As a result, the error estimate was not consistent, and adaptive stepping
could behave poorly at tighter tolerances.
I replaced it with coefficients that satisfy both order conditions:
Σ b̃ᵢ = 0
Σ b̃ᵢ cᵢ = 0
which correspond to a proper order-(p−1) embedded method for this
3rd-order scheme. With this change, the error estimate behaves correctly
and the solver shows expected convergence under tighter tolerances.
The underlying ARS343 method itself is unchanged — only the embedded
error estimator used for adaptivity.

@ChrisRackauckas
Copy link
Copy Markdown
Member

The paper does not have an error estimator, so just don't make it adaptive.

@singhharsh1708
Copy link
Copy Markdown
Author

@ChrisRackauckas got it that makes sense,
I’ve reverted the embedded error estimator change and made ARS343 non-adaptive.

@ChrisRackauckas
Copy link
Copy Markdown
Member

It doesn't look like it's not adaptive. It should set the trait based on whether btilde = nothing.

@singhharsh1708
Copy link
Copy Markdown
Author

@ChrisRackauckas updated ARS343 to set btilde = nothing so adaptivity is disabled via the tableau trait, and guarded the adaptive path accordingly.

@singhharsh1708 singhharsh1708 force-pushed the imex-tableau-prototype branch from 8cc211a to b8bdd3f Compare March 27, 2026 16:27
@singhharsh1708
Copy link
Copy Markdown
Author

@ChrisRackauckas migrated the remaining ESDIRK methods (ESDIRK54I8L2SA, ESDIRK436L2SA2, ESDIRK437L2SA, ESDIRK547L2SA2, ESDIRK659L2SA) to the tableau-based implementation.

@singhharsh1708
Copy link
Copy Markdown
Author

one follow up question :Methods such as ImplicitEuler, Trapezoid, TRBDF2, SDIRK2, SDIRK22, Cash4, Hairer4, and Hairer42 have different structural properties (non-IMEX formulations and different stage structures), and may not directly fit into the current ESDIRK IMEX abstraction. Would it make sense to:
introduce a separate tableau abstraction for pure SDIRK methods, or
generalize the current framework to support both IMEX and non-IMEX SDIRK methods?

@oscardssmith
Copy link
Copy Markdown
Member

IMO either is fine (but we probably want to do either as a separate PR.

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.

3 participants