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

@@ -7,7 +7,7 @@ function alg_cache(alg::Kvaerno3, u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
        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?

ChrisRackauckas · 2025-08-10T11:47:09Z

lib/OrdinaryDiffEqSDIRK/src/unified_sdirk_tableaus.jl

+    A_explicit::Union{SMatrix{S, S, T}, Nothing}
+    b_explicit::Union{SVector{S, T}, Nothing}
+    c_explicit::Union{SVector{S, T2}, Nothing}
+end


the tableaus are missing a step predictor coefficients part?

ChrisRackauckas · 2025-08-10T11:49:46Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+        z₁ = dt * integrator.fsalfirst
+    end
+
+    if s == 4  # 4-stage methods like KenCarp3, Kvaerno3


why this instead of just stages?

ChrisRackauckas · 2025-08-10T11:50:23Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+    end
+
+    if s == 4  # 4-stage methods like KenCarp3, Kvaerno3
+        z₂ = z₁


This is not generally a good prediction? It should be a linear combination from tableau coefficients?

KrishGaur1354 · 2025-08-17T16:01:33Z

@oscardssmith @ChrisRackauckas I would like to apologies for the delay from my side. I have made the necessary changes for the refactoring. Do let me know if there's anything else to be added/removed or if any change is required in this PR.

ChrisRackauckas · 2025-09-03T00:05:31Z

lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl

    γ, c = tab.γ, tab.γ
    nlsolver = build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false))
    Hairer4ConstantCache(nlsolver, tab)
 end

-@cache mutable struct Hairer4Cache{uType, rateType, uNoUnitsType, Tab, N} <:
+@cache mutable struct Hairer4Cache{uType, rateType, uNoUnitsType, Tab, N, StepLimiter} <:


The point of the refactor is that these types don't need to exist. If there's a single SDIRKCache that just has different values for the tableau, then every SDIRK algorithm gives the same cache type, and so they compile all the same. If you have separate cache types, they will all compile different codes. So this doesn't meet the criteria. I recommend looking at the Rosenbrock methods, where a good chunk of them are all the same cache https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl#L748-L832, and all the same perform_step!, and it's just written based on the tableau size. This makes it so that compiling for one Rosenbrock compiles for all of them, effectively removing the per-algorithm compilation.

Thank you for the feedback. I have implemented the refactor following the Rosenbrock pattern you referenced for all algo specific cache types have been removed and replaced with unified SDIRKCache and SDIRKConstantCache types that work across all SDIRK methods.

ChrisRackauckas · 2025-09-08T00:41:36Z

Test failure looks real.

ChrisRackauckas · 2025-09-09T12:01:54Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+    γ = tab.γ
+
+    # Check if this is an IMEX (split function) scheme
+    is_imex = integrator.f isa SplitFunction && tab.has_additive_splitting


tab.has_additive_splitting is not going to be inferred unless it's a value type, so is_imex is not compile-time information, meaning k_explicit is not going to be type-stable. I think the right answer is to make sure has_additive_splitting is type-level information.

ChrisRackauckas · 2025-09-09T12:05:16Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+                θ = c[3] / c[2]
+                α31 = ((1 + (-4θ + 3θ^2)) + (6θ * (1 - θ) / c[2]) * γ)
+                α32 = ((-2θ + 3θ^2)) + (6θ * (1 - θ) / c[2]) * γ
+                z₃ = α31 * z₁ + α32 * z₂


why is this hardcoded here instead of being in the tableau?

ChrisRackauckas · 2025-09-09T12:07:13Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+    generic_sdirk_perform_step!(integrator, cache, repeat_step)
+end
+
+@muladd function generic_additive_sdirk_perform_step!(integrator, cache::SDIRKConstantCache, repeat_step=false)


why is this kept separate? And what is meant by "additive SDRIK"?

ChrisRackauckas · 2025-09-09T12:08:04Z

lib/OrdinaryDiffEqSDIRK/src/generic_sdirk_perform_step.jl

+end
+
+# TRBDF2 special handling for constant cache
+@muladd function perform_step_trbdf2!(integrator, cache::SDIRKConstantCache, repeat_step=false)


What's special about TRBDF2? It should just be a standard ESDIRK

ChrisRackauckas · 2025-09-09T12:08:49Z

lib/OrdinaryDiffEqSDIRK/src/tableau_utils.jl

+    end
+end
+
+# simple 1-stage SDIRK perform step (ImplicitEuler, ImplicitMidpoint)


why would 1 and 2 stages need to be treated differently?

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

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

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

SDIRK tests now pass: JET and Aqua

02fa048

ChrisRackauckas reviewed Aug 10, 2025

View reviewed changes

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 4 commits August 10, 2025 22:16

ensure it passes the test on the buildkite

a8f533f

removed duplicate reexport values

7774123

Merge branch 'master' into master

383e106

made the changes to the toml file

0611319

ChrisRackauckas reviewed Sep 3, 2025

View reviewed changes

KrishGaur1354 and others added 2 commits September 8, 2025 00:01

Updated as required for it

60758a6

Merge branch 'master' into master

17538ba

Test Fixes

e8b11de

ChrisRackauckas reviewed Sep 9, 2025

View reviewed changes

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

Uh oh!

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!

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!

KrishGaur1354 commented Aug 17, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas commented Sep 8, 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!

ChrisRackauckas commented Sep 9, 2025

Uh oh!

Uh oh!

KrishGaur1354 commented Jul 27, 2025 •

edited

Loading