Refactoring and implementing unified tableau structure and generic perform_step! for SDIRK solver set #2779

KrishGaur1354 · 2025-07-27T18:40:27Z

This PR refactors the SDIRK solver implementation to use a unified tableau structure (SDIRKTableau)
and a single generic perform_step! function, consolidating 95% of SDIRK methods and reducing code
duplication. All methods now share consistent tableau coefficients while maintaining
full backward compatibility and performance. Special handling preserved for TRBDF2's unique structure.
KenCarp/Kvaerno methods retain separate implementations for additive splitting complexity.

Remaining Work:

KenCarp/Kvaerno Family (18 cache types):
- Kvaerno3/4/5 (ConstantCache + Cache variants)
- KenCarp3/4/5/47/58 (ConstantCache + Cache variants)
- CFNLIRK3 (ConstantCache + Cache variants)
- Kvaerno3/4/5 (ConstantCache + Cache variants)
- KenCarp3/4/5/47/58 (ConstantCache + Cache variants)
- CFNLIRK3 (ConstantCache + Cache variants)
Hairer Methods (2 cache types):
- Hairer4ConstantCache, Hairer4Cache
Special Array Variant:
- TRBDF2Cache{<:Array} (specialized version)
Duplication Removal of SDIRK related code
Re-run tests to check if working proper

Part of SciML Small Grants Program for OrdinaryDiffEq.jl solver set refactoring

- Implement SDIRKTableau struct consolidating all Butcher coefficients - Add generic_sdirk_perform_step! reducing duplication across 20+ methods - Migrate ImplicitEuler, SDIRK2/22, Cash4, all ESDIRK/SFSDIRK methods - Preserve TRBDF2 special handling with unified coefficients - Unify KenCarp/Kvaerno tableaus (separate perform_step! for additive splitting) - Reduce codebase by ~3000 lines while maintaining full compatibility - Pass all tests (JET, Aqua, integration) Part of SciML Small Grants: OrdinaryDiffEq.jl tableau refactoring Claimed by: Krish Gaur (July 4 - August 4, 2025) Ref: https://github.com/SciML/sciml.ai/blob/master/small_grants.md#refactor-ordinarydiffeqjl-solver-sets-to-reuse-perform_step-implementations-via-tableaus-100solver-set

…d Added missing tableau definitions for Hairer4, CFNLIRK3, KenCarp47, and KenCarp58

KrishGaur1354 · 2025-07-31T10:08:23Z

@ChrisRackauckas could you provide any suggestions or any changes I should make?

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

oscardssmith · 2025-07-31T12:11:33Z

The primary thing this PR seems to be missing is a bunch of deletions. If this PR unifies the SDIRK tableau structure and perform steps, why isn't it deleting the old tableaus and perform step methods?

ChrisRackauckas · 2025-07-31T12:18:24Z

What's the unique structure? TRBDF2 and KenCarps should be able to join just fine. Any step predictor here is just a linear combination of previous values, which these match, so you just need an alpha array of step predictor coefficients.

KrishGaur1354 · 2025-08-01T17:23:48Z

What's the unique structure? TRBDF2 and KenCarps should be able to join just fine. Any step predictor here is just a linear combination of previous values, which these match, so you just need an alpha array of step predictor coefficients.

Sure, I'll just ensure that the TRBDF2 and KenCarp methods would now be unified under consistent structure, which will utilize shared tableau and step predictor coefficients for it.

oscardssmith · 2025-08-04T14:41:02Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+    d = T(1 - sqrt(2) / 2) 
+    ω = T(sqrt(2) / 4)
+
+    A = @SMatrix [0  0  0;


are static arrays actually faster here than just normal Vector/Matrix?

Yes, I did found static arrays being significantly faster than normal Vector/Matrix for SDIRK tableau operations by almost 1.8 times faster.

That's interesting (and somewhat unfortunate). I would have expected the nonlinear solve cost to dominate.

Atleast the test that I did, Static Array was faster but the construction was slower due to compilation overhead in this case.

…ompat; bump version to 1.5.0

ChrisRackauckas · 2025-08-10T11:45:42Z

lib/OrdinaryDiffEqSDIRK/src/kencarp_kvaerno_caches.jl

        uprev, uprev2, f, t, dt, reltol, p, calck,
        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    tab = Kvaerno3Tableau(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+    tab = get_sdirk_tableau(:Kvaerno3, constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))


What's the reason for Symbol rather than alg dispatch here?

Switched to get_sdirk_tableau(:Kvaerno3, ...) for uniformity during the refactor. I can change it if you want?

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

ChrisRackauckas · 2025-09-09T12:10:47Z

Yes this is right on track. Looking fairly close. Left a bunch of comments.

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

ChrisRackauckas · 2025-09-29T06:24:59Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+        # For IMEX schemes, include explicit error contributions
+        # (This would require additional error estimation tableau components)


This is the same as before?

ChrisRackauckas · 2025-09-29T06:25:42Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+# Dispatch for unified caches - all SDIRK algorithms use these same cache types
+@muladd function perform_step!(integrator, cache::SDIRKCache, repeat_step=false)
+    generic_sdirk_perform_step!(integrator, cache, repeat_step)
+end
+
+@muladd function perform_step!(integrator, cache::SDIRKConstantCache, repeat_step=false)
+    generic_sdirk_perform_step!(integrator, cache, repeat_step)


Can just move the code into here instead of having a different function call?

ChrisRackauckas · 2025-09-29T06:26:39Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+    is_A_stable::Bool
+    is_L_stable::Bool
+    predictor_type::Symbol
+    A_explicit::Union{SMatrix{S, S, T}, Nothing}


It's not type-stable to grab these: should just parameterize this?

ChrisRackauckas · 2025-09-29T06:27:22Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+                 b_embed=b_embed, embedded_order=3,
+                 is_fsal=false, is_stiffly_accurate=true,
+                 is_A_stable=true, is_L_stable=true,
+                 predictor_type=:trbdf2_special)


Did TRBDF2 not have alpha_preds?

ChrisRackauckas · 2025-09-29T06:28:09Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+    c = @SVector [T2(1.0)]
+    γ = T(1.0)
+
+    SDIRKTableau(A, b, c, γ, 1;


How is the special error estimators of ImplicitEuler and Trapezoidal handled?

ChrisRackauckas · 2025-09-29T06:29:04Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+                 is_A_stable=true, is_L_stable=false,
+                 predictor_type=:kencarp_additive,
+                 has_additive_splitting=true,
+                 A_explicit=A_explicit, b_explicit=b_explicit, c_explicit=c_explicit)


This didn't have alpha preds?

ChrisRackauckas · 2025-09-29T06:30:27Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+                                                               A_explicit, b_explicit, c_explicit, α_pred)
+end
+
+function TRBDF2Tableau_unified(T=Float64, T2=Float64)


Suggested change

function TRBDF2Tableau_unified(T=Float64, T2=Float64)

function TRBDF2Tableau_unified(::Type{T}=Float64, ::Type{T2}=Float64) where {T,T2}

Need to force specialization on these in order to get inference.

ChrisRackauckas · 2025-09-29T06:30:54Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+    end
+end
+
+get_sdirk_tableau(alg, T=Float64, T2=Float64) = get_sdirk_tableau(nameof(typeof(alg)), T, T2)


Suggested change

get_sdirk_tableau(alg, T=Float64, T2=Float64) = get_sdirk_tableau(nameof(typeof(alg)), T, T2)

get_sdirk_tableau(alg, ::Type{T}=Float64, ::Type{T2}=Float64) where {T,T2} = get_sdirk_tableau(nameof(typeof(alg)), T, T2)

ChrisRackauckas · 2025-09-29T06:31:13Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+                 predictor_type=:hermite)
+end
+
+function get_sdirk_tableau(alg::Symbol, T=Float64, T2=Float64)


Suggested change

function get_sdirk_tableau(alg::Symbol, T=Float64, T2=Float64)

function get_sdirk_tableau(alg::Symbol, ::Type{T}=Float64, ::Type{T2}=Float64) where {T,T2}

ChrisRackauckas · 2025-09-29T06:33:09Z

Okay this is generally looking right. Left some comments. Probably just needs one more small round? Most are just little inference things. I think the only thing that's not trivial inference stuff is:

ImplicitEuler and Trapezoidal had the SPICE error estimators, how are those handled? I don't see that in there?
There are some methods that had built-in step predictors like TRBDF2 and KenCarps that I don't see changed to using the alphas? Curious if that was intentional.

ChrisRackauckas · 2025-11-06T10:48:26Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+                 predictor_type=:hermite)
+end
+
+function get_sdirk_tableau(alg::Symbol, ::Type{T}=Float64, ::Type{T2}=Float64) where {T, T2}


This is a bit wonky. It should instead just be like SDIRKTableau(KenCarp4(), ...) i.e. dispatch on the algorithm rather than different functions and this would all go away. But I can just refactor that at some later point. I wouldn't worry about this for now but it is a very weird way to do dispatch.

ChrisRackauckas · 2025-11-06T10:51:01Z

lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl

-    SFSDIRK5Cache(u, uprev, fsalfirst, z₁, z₂, z₃, z₄, z₅, z₆, atmp, nlsolver, tab)
-end
+# Unified alg_cache functions for mutable cache
+function alg_cache(alg::Union{ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22, SSPSDIRK2,


Is there not an abstract type for all of the SDIRK methods to dispatch with? That seems much more sane than a big union

ChrisRackauckas · 2025-11-06T10:52:13Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+    step_limiter!(u, integrator, p, t + dt)
+
+    if integrator.opts.adaptive
+        if tab.has_spice_error && integrator.success_iter > 0


It seems the Trapezoid one is missing? This is just the implicit euler one?

and it's probably best to just if integrator.alg isa ImplicitEuler here, since it only applies to exactly that truncation error estimate.

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

ChrisRackauckas · 2025-11-06T10:57:43Z

I think this looks good?

@oscardssmith take a look as well.

And @KrishGaur1354 run some of the stiff ODE benchmarks before and after like https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/Hires/ and https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/Pollution/ to see if there are any regressions, TRBDF2, ImplicitEuler, Trapezoid, and KenCarp4 in particular.

There are a few refactor things, specifically how the tableaus dispatch and how the cache dispatch is written that I noted, but those are relatively minor and I'd be willing to just Claude that after merging if there's no regressions and if @oscardssmith also approves.

KrishGaur1354 · 2025-11-07T12:07:46Z

Sure @ChrisRackauckas , I'll run the stiff ODE benchmarks with TRBDF2, impliciteuler, trapezoid, and KenCarp4 before and after the refactor, and report any regressions.
If everything looks stable, I’ll update here so we can move forward with the merge.

Also, apologies for the long gaps between updates, was a bit occupied, but I'm fully back on it now.

Let me know if you'd like me to include any additional benchmark cases.

ChrisRackauckas · 2025-11-09T11:30:56Z

No worries, take the time you need. Near the end now!

…proper isa checks - Implemented α_pred predictor coefficients for TRBDF2, KenCarp3, and Kvaerno3 - Fixed SPICE error estimator broadcast bug by removing incorrect internalnorm calls - Added ImplicitEuler cache compatibility for BDF module integration

…fEq.jl

KrishGaur1354 · 2025-11-12T06:25:01Z

No worries, take the time you need. Near the end now!

So, @ChrisRackauckas I have addressed some of the things you recommended to look into:

ImplicitEuler and Trapezoid SPICE error estimators which is now implemented with isa checks
TRBDF2, KenCarp3, Kvaerno3 which now use α_pred predictors
Package compiles successfully

Place where I got stuck: Trapezoid and TRBDF2 fail on first step because SPICE estimators need success_iter > 0. With no history and no embedded method, they set EEst = 1 causing immediate failure. So, how should these handle the first step before SPICE history exists?

ImplicitEuler and KenCarp4 work correctly. Once the first-step case is resolved, I can run full benchmarks.

KrishGaur1354 marked this pull request as ready for review July 28, 2025 12:55

KrishGaur1354 and others added 2 commits July 29, 2025 12:30

Merge branch 'SciML:master' into master

bbe9847

Implemented generic additive splitting for KenCarp/Kvaerno methods an…

9d1201c

…d Added missing tableau definitions for Hairer4, CFNLIRK3, KenCarp47, and KenCarp58

oscardssmith reviewed Jul 31, 2025

View reviewed changes

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl Outdated Show resolved Hide resolved

Merge branch 'SciML:master' into master

8288658

KrishGaur1354 and others added 3 commits August 4, 2025 12:48

Merge branch 'SciML:master' into master

ba76fb1

Made the requested changes

1c64aad

comments added

d69d9e6

oscardssmith reviewed Aug 4, 2025

View reviewed changes

Fix: Add missing perform_step! methods for SDIRK algorithms

e2aec3b

KrishGaur1354 requested a review from oscardssmith August 6, 2025 16:58

Fixed: SSPSDIRK2Cache and TrapezoidCache Error should be fixed for now

38c7af4

KrishGaur1354 changed the title ~~[WIP]: Refactoring and implementing unified tableau structure and generic perform_step! for SDIRK solver set~~ Refactoring and implementing unified tableau structure and generic perform_step! for SDIRK solver set Aug 6, 2025

KrishGaur1354 and others added 5 commits August 8, 2025 23:24

SDIRK Project.toml: resolve merge conflicts, add test/extras, align c…

11cc123

…ompat; bump version to 1.5.0

Resolve all merge conflicts in OrdinaryDiffEqSDIRK/Project.toml

498f3ac

.

885c77e

Merge branch 'master' into master

8cb50fd

Merge branch 'SciML:master' into master

12e7c16

ChrisRackauckas reviewed Aug 10, 2025

View reviewed changes

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl Show resolved Hide resolved

SDIRK tests now pass: JET and Aqua

02fa048

ChrisRackauckas reviewed Aug 10, 2025

View reviewed changes

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Aug 10, 2025

View reviewed changes

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl Outdated Show resolved Hide resolved

KrishGaur1354 and others added 2 commits August 10, 2025 17:54

Replaced the hardcoded element and extended SDIRKTableau

fcb60ef

Merge branch 'master' into master

8b8b00d

KrishGaur1354 and others added 3 commits September 17, 2025 19:12

Merge branch 'SciML:master' into master

f03df64

Made suggested changes

8152d18

Fixed SDIRKConstantCache naming conflict

eedab04

ChrisRackauckas reviewed Sep 29, 2025

View reviewed changes

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl Show resolved Hide resolved

ChrisRackauckas reviewed Sep 29, 2025

View reviewed changes

oscardssmith mentioned this pull request Sep 29, 2025

WIP: Implementation of SDIRK22 #2876

Draft

5 tasks

KrishGaur1354 and others added 2 commits November 1, 2025 00:19

Merge branch 'SciML:master' into master

63abf0f

Made the requested changes

9fcc496

ChrisRackauckas reviewed Nov 6, 2025

View reviewed changes

Delete lib/OrdinaryDiffEqSDIRK/src/sdirk_perform_step.jl

9f162fa

ChrisRackauckas reviewed Nov 6, 2025

View reviewed changes

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl Show resolved Hide resolved

oscardssmith self-assigned this Nov 7, 2025

curdlinrope added 2 commits November 12, 2025 11:43

Merge branch 'master' of https://github.com/KrishGaur1354/OrdinaryDif…

df753ce

…fEq.jl

		# For IMEX schemes, include explicit error contributions
		# (This would require additional error estimation tableau components)

	function TRBDF2Tableau_unified(T=Float64, T2=Float64)
	function TRBDF2Tableau_unified(::Type{T}=Float64, ::Type{T2}=Float64) where {T,T2}

	get_sdirk_tableau(alg, T=Float64, T2=Float64) = get_sdirk_tableau(nameof(typeof(alg)), T, T2)
	get_sdirk_tableau(alg, ::Type{T}=Float64, ::Type{T2}=Float64) where {T,T2} = get_sdirk_tableau(nameof(typeof(alg)), T, T2)

	function get_sdirk_tableau(alg::Symbol, T=Float64, T2=Float64)
	function get_sdirk_tableau(alg::Symbol, ::Type{T}=Float64, ::Type{T2}=Float64) where {T,T2}

Uh oh!

Refactoring and implementing unified tableau structure and generic perform_step! for SDIRK solver set #2779

Are you sure you want to change the base?

Refactoring and implementing unified tableau structure and generic perform_step! for SDIRK solver set #2779

Conversation

KrishGaur1354 commented Jul 27, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

KrishGaur1354 commented Jul 31, 2025

Uh oh!

Uh oh!

oscardssmith commented Jul 31, 2025

Uh oh!

ChrisRackauckas commented Jul 31, 2025

Uh oh!

KrishGaur1354 commented Aug 1, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

ChrisRackauckas commented Sep 9, 2025

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas commented Sep 29, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

ChrisRackauckas commented Nov 6, 2025

Uh oh!

KrishGaur1354 commented Nov 7, 2025

Uh oh!

ChrisRackauckas commented Nov 9, 2025

Uh oh!

KrishGaur1354 commented Nov 12, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

KrishGaur1354 commented Jul 27, 2025 •

edited

Loading