Refactoring and implementing unified tableau structure and generic perform_step! for SDIRK solver set#2779
Refactoring and implementing unified tableau structure and generic perform_step! for SDIRK solver set#2779KrishGaur1354 wants to merge 39 commits intoSciML:masterfrom
Conversation
- 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
|
@ChrisRackauckas could you provide any suggestions or any changes I should make? |
|
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? |
|
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. |
| d = T(1 - sqrt(2) / 2) | ||
| ω = T(sqrt(2) / 4) | ||
|
|
||
| A = @SMatrix [0 0 0; |
There was a problem hiding this comment.
are static arrays actually faster here than just normal Vector/Matrix?
There was a problem hiding this comment.
Yes, I did found static arrays being significantly faster than normal Vector/Matrix for SDIRK tableau operations by almost 1.8 times faster.
There was a problem hiding this comment.
That's interesting (and somewhat unfortunate). I would have expected the nonlinear solve cost to dominate.
There was a problem hiding this comment.
Atleast the test that I did, Static Array was faster but the construction was slower due to compilation overhead in this case.
KrishGaur1354
changed the title
| 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)) |
There was a problem hiding this comment.
What's the reason for Symbol rather than alg dispatch here?
There was a problem hiding this comment.
Switched to get_sdirk_tableau(:Kvaerno3, ...) for uniformity during the refactor. I can change it if you want?
| 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, |
There was a problem hiding this comment.
Is there not an abstract type for all of the SDIRK methods to dispatch with? That seems much more sane than a big union
| step_limiter!(u, integrator, p, t + dt) | ||
|
|
||
| if integrator.opts.adaptive | ||
| if tab.has_spice_error && integrator.success_iter > 0 |
There was a problem hiding this comment.
It seems the Trapezoid one is missing? This is just the implicit euler one?
There was a problem hiding this comment.
and it's probably best to just if integrator.alg isa ImplicitEuler here, since it only applies to exactly that truncation error estimate.
|
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. |
|
Sure @ChrisRackauckas , I'll run the stiff ODE benchmarks with TRBDF2, impliciteuler, trapezoid, and KenCarp4 before and after the refactor, and report any regressions. 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. |
|
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
So, @ChrisRackauckas I have addressed some of the things you recommended to look into:
Place where I got stuck: Trapezoid and TRBDF2 fail on first step because SPICE estimators need ImplicitEuler and KenCarp4 work correctly. Once the first-step case is resolved, I can run full benchmarks. |
I thought EEst = 1 would always accept, so it accepts the first step if the nonlinear solver converges? I believe that is how it was intended before. |
|
@oscardssmith should take a deep look |
|
oh did we mix up a |
|
@ KrishGaur1354 I've looked into this and |
…plicit stages. Added regression tests for SDIRK methods with explicit first stages to ensure proper step acceptance.
|
@ChrisRackauckas @oscardssmith Pushed a fix so SDIRK explicit-first stages skip the Newton solve (A[i,i]==0) and added a regression that checks Trapezoid/TRBDF2 accept the first step. Simple solves now succeed |
|
Rebase onto latest mater and tests should be passing. |
|
@ChrisRackauckas @oscardssmith Rebased onto latest master and fixed compatibility issues:
|
|
It looks like a lot of the convergence tests are failing. |
|
Rebase onto master, tests should be fixed there so this should just show up as passing and merge. |
Rebase and cleanup of PR SciML#2779 onto current master. - Implement SDIRKTableau struct consolidating all Butcher coefficients - Add generic perform_step! reducing duplication across 20+ methods - Migrate ImplicitEuler, SDIRK2, Cash4, all ESDIRK/SFSDIRK methods - Preserve TRBDF2 special handling with unified coefficients - Unify KenCarp/Kvaerno tableaus with additive splitting support - Add Cash4Tableau_unified with proper embedding support - Add verbose parameter to alg_cache for upstream compatibility - Add regression tests for SDIRK methods with explicit first stages Note: Performance regression identified - generic perform_step! uses Vector{typeof(u)} for stage storage and dynamic dispatch which causes significant overhead vs the original flat-field specialized methods. This needs optimization before merging. Co-Authored-By: Krish Gaur <csrealracing13d@gmail.com> Co-Authored-By: curd <curdlinrope@gmail.com> Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
4 participants
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 codeduplication. 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:
Part of SciML Small Grants Program for OrdinaryDiffEq.jl solver set refactoring