From 09177724eb905188b65f1ee33ed4472890ed2c93 Mon Sep 17 00:00:00 2001 From: Harsh Singh Date: Mon, 30 Mar 2026 18:07:52 +0530 Subject: [PATCH 1/6] refactor(Rosenbrock): migrate Rosenbrock23/32 to generic RodasTableau framework Unify Rosenbrock23 and Rosenbrock32 with the existing generic Rosenbrock/Rodas infrastructure by expressing their coefficients as RodasTableau entries. This removes ~900 lines of duplicated cache structs, perform_step!, interpolation, and addsteps code. Co-Authored-By: Claude Opus 4.6 --- .../src/interp_func.jl | 15 - .../src/rosenbrock_caches.jl | 279 +---------- .../src/rosenbrock_interpolants.jl | 136 ------ .../src/rosenbrock_perform_step.jl | 454 +----------------- .../src/rosenbrock_tableaus.jl | 92 +++- .../src/stiff_addsteps.jl | 52 -- 6 files changed, 89 insertions(+), 939 deletions(-) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl b/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl index efc66a63a1e..314f2792aa9 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl @@ -1,18 +1,3 @@ -function SciMLBase.interp_summary( - ::Type{cacheType}, - dense::Bool - ) where { - cacheType <: - Union{ - Rosenbrock23ConstantCache, - Rosenbrock32ConstantCache, - Rosenbrock23Cache, - Rosenbrock32Cache, - }, - } - return dense ? "specialized 2nd order \"free\" stiffness-aware interpolation" : - "1st order linear" -end function SciMLBase.interp_summary( ::Type{cacheType}, dense::Bool diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl index 5985e4ad254..4a37aab74c2 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl @@ -115,277 +115,6 @@ struct RosenbrockCombinedConstantCache{TF, UF, Tab, JType, WType, F, AD, JRType} jac_reuse::JRType end -@cache mutable struct Rosenbrock23Cache{ - uType, rateType, uNoUnitsType, JType, WType, - TabType, TFType, UFType, F, JCType, GCType, - RTolType, A, AV, StepLimiter, StageLimiter, JRType, - } <: RosenbrockMutableCache - u::uType - uprev::uType - k₁::rateType - k₂::rateType - k₃::rateType - du1::rateType - du2::rateType - f₁::rateType - fsalfirst::rateType - fsallast::rateType - dT::rateType - J::JType - W::WType - tmp::rateType - atmp::uNoUnitsType - weight::uNoUnitsType - tab::TabType - tf::TFType - uf::UFType - linsolve_tmp::rateType - linsolve::F - jac_config::JCType - grad_config::GCType - reltol::RTolType - alg::A - algebraic_vars::AV - step_limiter!::StepLimiter - stage_limiter!::StageLimiter - jac_reuse::JRType -end - -@cache mutable struct Rosenbrock32Cache{ - uType, rateType, uNoUnitsType, JType, WType, - TabType, TFType, UFType, F, JCType, GCType, - RTolType, A, AV, StepLimiter, StageLimiter, JRType, - } <: RosenbrockMutableCache - u::uType - uprev::uType - k₁::rateType - k₂::rateType - k₃::rateType - du1::rateType - du2::rateType - f₁::rateType - fsalfirst::rateType - fsallast::rateType - dT::rateType - J::JType - W::WType - tmp::rateType - atmp::uNoUnitsType - weight::uNoUnitsType - tab::TabType - tf::TFType - uf::UFType - linsolve_tmp::rateType - linsolve::F - jac_config::JCType - grad_config::GCType - reltol::RTolType - alg::A - algebraic_vars::AV - step_limiter!::StepLimiter - stage_limiter!::StageLimiter - jac_reuse::JRType -end - -function alg_cache( - alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}, verbose - ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - k₁ = zero(rate_prototype) - k₂ = zero(rate_prototype) - k₃ = zero(rate_prototype) - du1 = zero(rate_prototype) - du2 = zero(rate_prototype) - # f₀ = zero(u) fsalfirst - f₁ = zero(rate_prototype) - fsalfirst = zero(rate_prototype) - fsallast = zero(rate_prototype) - dT = zero(rate_prototype) - tmp = zero(rate_prototype) - atmp = similar(u, uEltypeNoUnits) - recursivefill!(atmp, false) - weight = similar(u, uEltypeNoUnits) - recursivefill!(weight, false) - tab = Rosenbrock23Tableau(constvalue(uBottomEltypeNoUnits)) - tf = TimeGradientWrapper(f, uprev, p) - uf = UJacobianWrapper(f, t, p) - linsolve_tmp = zero(rate_prototype) - - grad_config = build_grad_config(alg, f, tf, du1, t) - jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2) - - J, W = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true)) - - linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp)) - Pl, - Pr = wrapprecs( - alg.precs( - W, nothing, u, p, t, nothing, nothing, nothing, - nothing - )..., weight, tmp - ) - linsolve = init( - linprob, alg.linsolve, alias = LinearAliasSpecifier(alias_A = true, alias_b = true), - Pl = Pl, Pr = Pr, - abstol = reltol, reltol = reltol, - assumptions = LinearSolve.OperatorAssumptions(true), - verbose = verbose.linear_verbosity - ) - - algebraic_vars = f.mass_matrix === I ? nothing : - [all(iszero, x) for x in eachcol(f.mass_matrix)] - - return Rosenbrock23Cache( - u, uprev, k₁, k₂, k₃, du1, du2, f₁, - fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf, - linsolve_tmp, - linsolve, jac_config, grad_config, reltol, alg, algebraic_vars, alg.step_limiter!, - alg.stage_limiter!, _make_jac_reuse_state(zero(dt), alg.max_jac_age) - ) -end - -function alg_cache( - alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}, verbose - ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - k₁ = zero(rate_prototype) - k₂ = zero(rate_prototype) - k₃ = zero(rate_prototype) - du1 = zero(rate_prototype) - du2 = zero(rate_prototype) - # f₀ = zero(u) fsalfirst - f₁ = zero(rate_prototype) - fsalfirst = zero(rate_prototype) - fsallast = zero(rate_prototype) - dT = zero(rate_prototype) - tmp = zero(rate_prototype) - atmp = similar(u, uEltypeNoUnits) - recursivefill!(atmp, false) - weight = similar(u, uEltypeNoUnits) - recursivefill!(weight, false) - tab = Rosenbrock32Tableau(constvalue(uBottomEltypeNoUnits)) - - tf = TimeGradientWrapper(f, uprev, p) - uf = UJacobianWrapper(f, t, p) - linsolve_tmp = zero(rate_prototype) - - grad_config = build_grad_config(alg, f, tf, du1, t) - jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2) - - J, W = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true)) - - linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp)) - - Pl, - Pr = wrapprecs( - alg.precs( - W, nothing, u, p, t, nothing, nothing, nothing, - nothing - )..., weight, tmp - ) - linsolve = init( - linprob, alg.linsolve, alias = LinearAliasSpecifier(alias_A = true, alias_b = true), - Pl = Pl, Pr = Pr, - abstol = reltol, reltol = reltol, - assumptions = LinearSolve.OperatorAssumptions(true), - verbose = verbose.linear_verbosity - ) - - algebraic_vars = f.mass_matrix === I ? nothing : - [all(iszero, x) for x in eachcol(f.mass_matrix)] - - return Rosenbrock32Cache( - u, uprev, k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, - tmp, atmp, weight, tab, tf, uf, linsolve_tmp, linsolve, jac_config, - grad_config, reltol, alg, algebraic_vars, alg.step_limiter!, alg.stage_limiter!, - _make_jac_reuse_state(zero(dt), alg.max_jac_age) - ) -end - -struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F, AD, JRType} <: - RosenbrockConstantCache - c₃₂::T - d::T - tf::TF - uf::UF - J::JType - W::WType - linsolve::F - autodiff::AD - jac_reuse::JRType -end - -function Rosenbrock23ConstantCache( - ::Type{T}, tf, uf, J, W, linsolve, autodiff, max_jac_age::Int = 20 - ) where {T} - tab = Rosenbrock23Tableau(T) - return Rosenbrock23ConstantCache( - tab.c₃₂, tab.d, tf, uf, J, W, linsolve, autodiff, - _make_jac_reuse_state(zero(T), max_jac_age) - ) -end - -function alg_cache( - alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}, verbose - ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - tf = TimeDerivativeWrapper(f, u, p) - uf = UDerivativeWrapper(f, t, p) - J, W = build_J_W(alg, u, uprev, p, t, dt, f, nothing, uEltypeNoUnits, Val(false)) - linprob = nothing #LinearProblem(W,copy(u); u0=copy(u)) - linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true) - return Rosenbrock23ConstantCache( - constvalue(uBottomEltypeNoUnits), tf, uf, J, W, linsolve, - alg_autodiff(alg), alg.max_jac_age - ) -end - -struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F, AD, JRType} <: - RosenbrockConstantCache - c₃₂::T - d::T - tf::TF - uf::UF - J::JType - W::WType - linsolve::F - autodiff::AD - jac_reuse::JRType -end - -function Rosenbrock32ConstantCache( - ::Type{T}, tf, uf, J, W, linsolve, autodiff, max_jac_age::Int = 20 - ) where {T} - tab = Rosenbrock32Tableau(T) - return Rosenbrock32ConstantCache( - tab.c₃₂, tab.d, tf, uf, J, W, linsolve, autodiff, - _make_jac_reuse_state(zero(T), max_jac_age) - ) -end - -function alg_cache( - alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}, verbose - ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - tf = TimeDerivativeWrapper(f, u, p) - uf = UDerivativeWrapper(f, t, p) - J, W = build_J_W(alg, u, uprev, p, t, dt, f, nothing, uEltypeNoUnits, Val(false)) - linprob = nothing #LinearProblem(W,copy(u); u0=copy(u)) - linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true) - return Rosenbrock32ConstantCache( - constvalue(uBottomEltypeNoUnits), tf, uf, J, W, linsolve, - alg_autodiff(alg), alg.max_jac_age - ) -end - ### Rodas4+ methods and consolidated Rosenbrock methods (using RodasTableau) # Helper accessors for step_limiter!/stage_limiter! — algorithms that have these fields @@ -438,6 +167,8 @@ tabtype(::GRK4T) = GRK4TRodasTableau tabtype(::GRK4A) = GRK4ARodasTableau tabtype(::Ros4LStab) = Ros4LStabRodasTableau tabtype(::RosenbrockW6S4OS) = RosenbrockW6S4OSRodasTableau +tabtype(::Rosenbrock23) = Rosenbrock23RodasTableau +tabtype(::Rosenbrock32) = Rosenbrock32RodasTableau # Union of all algorithms using RodasTableau-based RosenbrockCache const RodasTableauAlgorithms = Union{ @@ -449,6 +180,7 @@ const RodasTableauAlgorithms = Union{ ROS34PRw, ROS3PRL, ROS3PRL2, ROK4a, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A, Ros4LStab, RosenbrockW6S4OS, + Rosenbrock23, Rosenbrock32, } function alg_cache( @@ -562,10 +294,7 @@ function alg_cache( end function get_fsalfirstlast( - cache::Union{ - Rosenbrock23Cache, Rosenbrock32Cache, - RosenbrockCache, - }, + cache::RosenbrockCache, u ) return (cache.fsalfirst, cache.fsallast) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_interpolants.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_interpolants.jl index 55578f44c1c..52a6e70930b 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_interpolants.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_interpolants.jl @@ -1,7 +1,5 @@ ### Fallbacks to capture ROSENBROCKS_WITH_INTERPOLATIONS = Union{ - Rosenbrock23ConstantCache, Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, RosenbrockCombinedConstantCache, RosenbrockCache, HybridExplicitImplicitConstantCache, HybridExplicitImplicitCache, @@ -25,140 +23,6 @@ end #### -""" -From MATLAB ODE Suite by Shampine -""" -@def rosenbrock2332unpack begin - if cache isa OrdinaryDiffEqMutableCache - d = cache.tab.d - else - d = cache.d - end -end - -@def rosenbrock2332pre0 begin - @rosenbrock2332unpack - c1 = Θ * (1 - Θ) / (1 - 2d) - c2 = Θ * (Θ - 2d) / (1 - 2d) -end - -@muladd function _ode_interpolant( - Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, - Rosenbrock32ConstantCache, - }, idxs::Nothing, - T::Type{Val{0}}, differential_vars - ) - @rosenbrock2332pre0 - @inbounds y₀ + dt * (c1 * k[1] + c2 * k[2]) -end - -@muladd function _ode_interpolant( - Θ, dt, y₀, y₁, k, - cache::Union{Rosenbrock23Cache, Rosenbrock32Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars - ) - @rosenbrock2332pre0 - @inbounds @.. y₀ + dt * (c1 * k[1] + c2 * k[2]) -end - -@muladd function _ode_interpolant( - Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{0}}, differential_vars - ) - @rosenbrock2332pre0 - @.. y₀[idxs] + dt * (c1 * k[1][idxs] + c2 * k[2][idxs]) -end - -@muladd function _ode_interpolant!( - out, Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, - Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs::Nothing, T::Type{Val{0}}, differential_vars - ) - @rosenbrock2332pre0 - @inbounds @.. out = y₀ + dt * (c1 * k[1] + c2 * k[2]) - out -end - -@muladd function _ode_interpolant!( - out, Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, - Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{0}}, differential_vars - ) - @rosenbrock2332pre0 - @views @.. out = y₀[idxs] + dt * (c1 * k[1][idxs] + c2 * k[2][idxs]) - out -end - -# First Derivative of the dense output -@def rosenbrock2332pre1 begin - @rosenbrock2332unpack - c1diff = (1 - 2 * Θ) / (1 - 2 * d) - c2diff = (2 * Θ - 2 * d) / (1 - 2 * d) -end - -@muladd function _ode_interpolant( - Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs::Nothing, T::Type{Val{1}}, differential_vars - ) - @rosenbrock2332pre1 - @.. c1diff * k[1] + c2diff * k[2] -end - -@muladd function _ode_interpolant( - Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{1}}, differential_vars - ) - @rosenbrock2332pre1 - @.. c1diff * k[1][idxs] + c2diff * k[2][idxs] -end - -@muladd function _ode_interpolant!( - out, Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, - Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs::Nothing, T::Type{Val{1}}, differential_vars - ) - @rosenbrock2332pre1 - @.. out = c1diff * k[1] + c2diff * k[2] - out -end - -@muladd function _ode_interpolant!( - out, Θ, dt, y₀, y₁, k, - cache::Union{ - Rosenbrock23ConstantCache, - Rosenbrock23Cache, - Rosenbrock32ConstantCache, Rosenbrock32Cache, - }, idxs, T::Type{Val{1}}, differential_vars - ) - @rosenbrock2332pre1 - @views @.. out = c1diff * k[1][idxs] + c2diff * k[2][idxs] - out -end - -""" -From MATLAB ODE Suite by Shampine -""" - @muladd function _ode_interpolant( Θ, dt, y₀, diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl index 368f57169fc..bce48e7a770 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl @@ -1,442 +1,3 @@ -function initialize!( - integrator, cache::Union{ - Rosenbrock23Cache, - Rosenbrock32Cache, - } - ) - integrator.kshortsize = 2 - (; k₁, k₂, fsalfirst, fsallast) = cache - resize!(integrator.k, integrator.kshortsize) - integrator.k[1] = k₁ - integrator.k[2] = k₂ - integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) - return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) -end - -function initialize!( - integrator, - cache::Union{ - Rosenbrock23ConstantCache, - Rosenbrock32ConstantCache, - } - ) - integrator.kshortsize = 2 - integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) - integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - # Avoid undefined entries if k is an array of arrays - integrator.fsallast = zero(integrator.fsalfirst) - integrator.k[1] = zero(integrator.fsalfirst) - return integrator.k[2] = zero(integrator.fsalfirst) -end - -@muladd function perform_step!(integrator, cache::Rosenbrock23Cache, repeat_step = false) - (; t, dt, uprev, u, f, p, opts) = integrator - (; - k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, - linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter!, - ) = cache - (; c₃₂, d) = cache.tab - - # Assignments - sizeu = size(u) - mass_matrix = integrator.f.mass_matrix - - # Precalculations - dtγ = dt * d - neginvdtγ = -inv(dtγ) - dto2 = dt / 2 - dto6 = dt / 6 - - if repeat_step - f(integrator.fsalfirst, uprev, p, t) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - end - - new_W = calc_rosenbrock_differentiation!(integrator, cache, dtγ, dtγ, repeat_step) - - calculate_residuals!( - weight, fill!(weight, one(eltype(u))), uprev, uprev, - integrator.opts.abstol, integrator.opts.reltol, - integrator.opts.internalnorm, t - ) - - if repeat_step || !new_W - linres = dolinsolve( - integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp), - du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight, - solverdata = (; gamma = dtγ) - ) - else - linres = dolinsolve( - integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp), - du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight, - solverdata = (; gamma = dtγ) - ) - end - - vecu = _vec(linres.u) - veck₁ = _vec(k₁) - - @.. veck₁ = vecu * neginvdtγ - integrator.stats.nsolve += 1 - - @.. u = uprev + dto2 * k₁ - stage_limiter!(u, integrator, p, t + dto2) - f(f₁, u, p, t + dto2) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - copyto!(tmp, k₁) - else - mul!(_vec(tmp), mass_matrix, _vec(k₁)) - end - - @.. linsolve_tmp = f₁ - tmp - - linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp)) - vecu = _vec(linres.u) - veck₂ = _vec(k₂) - - @.. veck₂ = vecu * neginvdtγ + veck₁ - integrator.stats.nsolve += 1 - - @.. u = uprev + dt * k₂ - stage_limiter!(u, integrator, p, t + dt) - step_limiter!(u, integrator, p, t + dt) - - if integrator.opts.adaptive - f(fsallast, u, p, t + dt) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - @.. broadcast = false linsolve_tmp = fsallast - c₃₂ * (k₂ - f₁) - - 2(k₁ - fsalfirst) + dt * dT - else - @.. broadcast = false du2 = c₃₂ * k₂ + 2k₁ - mul!(_vec(du1), mass_matrix, _vec(du2)) - @.. broadcast = false linsolve_tmp = fsallast - du1 + c₃₂ * f₁ + 2fsalfirst + - dt * dT - end - - linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp)) - vecu = _vec(linres.u) - veck3 = _vec(k₃) - @.. veck3 = vecu * neginvdtγ - - integrator.stats.nsolve += 1 - - if mass_matrix === I - @.. broadcast = false tmp = dto6 * (k₁ - 2 * k₂ + k₃) - else - veck₁ = _vec(k₁) - veck₂ = _vec(k₂) - veck₃ = _vec(k₃) - vectmp = _vec(tmp) - @.. broadcast = false vectmp = ifelse( - cache.algebraic_vars, - false, dto6 * (veck₁ - 2 * veck₂ + veck₃) - ) - end - calculate_residuals!( - atmp, tmp, uprev, u, integrator.opts.abstol, - integrator.opts.reltol, integrator.opts.internalnorm, t - ) - integrator.EEst = integrator.opts.internalnorm(atmp, t) - - if mass_matrix !== I - algvar = reshape(cache.algebraic_vars, size(u)) - invatol = inv(integrator.opts.abstol) - @.. atmp = ifelse(algvar, fsallast, false) * invatol - integrator.EEst += integrator.opts.internalnorm(atmp, t) - end - end - cache.linsolve = linres.cache -end - -@muladd function perform_step!(integrator, cache::Rosenbrock32Cache, repeat_step = false) - (; t, dt, uprev, u, f, p, opts) = integrator - (; - k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W, tmp, uf, tf, - linsolve_tmp, jac_config, atmp, weight, stage_limiter!, step_limiter!, - ) = cache - (; c₃₂, d) = cache.tab - - # Assignments - sizeu = size(u) - mass_matrix = integrator.f.mass_matrix - - # Precalculations - dtγ = dt * d - neginvdtγ = -inv(dtγ) - dto2 = dt / 2 - dto6 = dt / 6 - - if repeat_step - f(integrator.fsalfirst, uprev, p, t) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - end - - new_W = calc_rosenbrock_differentiation!(integrator, cache, dtγ, dtγ, repeat_step) - - calculate_residuals!( - weight, fill!(weight, one(eltype(u))), uprev, uprev, - integrator.opts.abstol, integrator.opts.reltol, - integrator.opts.internalnorm, t - ) - - if repeat_step || !new_W - linres = dolinsolve( - integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp), - du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight, - solverdata = (; gamma = dtγ) - ) - else - linres = dolinsolve( - integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp), - du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight, - solverdata = (; gamma = dtγ) - ) - end - - vecu = _vec(linres.u) - veck₁ = _vec(k₁) - - @.. veck₁ = vecu * neginvdtγ - integrator.stats.nsolve += 1 - - @.. broadcast = false u = uprev + dto2 * k₁ - stage_limiter!(u, integrator, p, t + dto2) - f(f₁, u, p, t + dto2) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - tmp .= k₁ - else - mul!(_vec(tmp), mass_matrix, _vec(k₁)) - end - - @.. broadcast = false linsolve_tmp = f₁ - tmp - - linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp)) - vecu = _vec(linres.u) - veck₂ = _vec(k₂) - - @.. veck₂ = vecu * neginvdtγ + veck₁ - integrator.stats.nsolve += 1 - - @.. tmp = uprev + dt * k₂ - stage_limiter!(u, integrator, p, t + dt) - f(fsallast, tmp, p, t + dt) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - @.. broadcast = false linsolve_tmp = fsallast - c₃₂ * (k₂ - f₁) - 2(k₁ - fsalfirst) + - dt * dT - else - @.. broadcast = false du2 = c₃₂ * k₂ + 2k₁ - mul!(_vec(du1), mass_matrix, _vec(du2)) - @.. broadcast = false linsolve_tmp = fsallast - du1 + c₃₂ * f₁ + 2fsalfirst + dt * dT - end - - linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp)) - vecu = _vec(linres.u) - veck3 = _vec(k₃) - - @.. veck3 = vecu * neginvdtγ - integrator.stats.nsolve += 1 - - @.. broadcast = false u = uprev + dto6 * (k₁ + 4k₂ + k₃) - - step_limiter!(u, integrator, p, t + dt) - - if integrator.opts.adaptive - @.. broadcast = false tmp = dto6 * (k₁ - 2 * k₂ + k₃) - calculate_residuals!( - atmp, tmp, uprev, u, integrator.opts.abstol, - integrator.opts.reltol, integrator.opts.internalnorm, t - ) - integrator.EEst = integrator.opts.internalnorm(atmp, t) - - if mass_matrix !== I - invatol = inv(integrator.opts.abstol) - @.. atmp = ifelse(cache.algebraic_vars, fsallast, false) * invatol - integrator.EEst += integrator.opts.internalnorm(atmp, t) - end - end - cache.linsolve = linres.cache -end - -@muladd function perform_step!( - integrator, cache::Rosenbrock23ConstantCache, - repeat_step = false - ) - (; t, dt, uprev, u, f, p) = integrator - (; c₃₂, d, tf, uf) = cache - - # Precalculations - dtγ = dt * d - neginvdtγ = -inv(dtγ) - dto2 = dt / 2 - dto6 = dt / 6 - - if repeat_step - integrator.fsalfirst = f(uprev, p, t) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - end - - mass_matrix = integrator.f.mass_matrix - - # Time derivative and W matrix (with Jacobian reuse for W-methods) - dT, W = calc_rosenbrock_differentiation(integrator, cache, dtγ, repeat_step) - if !issuccess_W(W) - integrator.EEst = 2 - return nothing - end - - k₁ = _reshape(W \ _vec((integrator.fsalfirst + dtγ * dT)), axes(uprev)) * neginvdtγ - integrator.stats.nsolve += 1 - tmp = @.. uprev + dto2 * k₁ - f₁ = f(tmp, p, t + dto2) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - k₂ = _reshape(W \ _vec(f₁ - k₁), axes(uprev)) - else - k₂ = _reshape(W \ _vec(f₁ - mass_matrix * k₁), axes(uprev)) - end - k₂ = @.. k₂ * neginvdtγ + k₁ - integrator.stats.nsolve += 1 - u = uprev + dt * k₂ - - if integrator.opts.adaptive - integrator.fsallast = f(u, p, t + dt) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - linsolve_tmp = @.. ( - integrator.fsallast - c₃₂ * (k₂ - f₁) - - 2 * (k₁ - integrator.fsalfirst) + dt * dT - ) - else - linsolve_tmp = mass_matrix * (@.. c₃₂ * k₂ + 2 * k₁) - linsolve_tmp = @.. ( - integrator.fsallast - linsolve_tmp + - c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT - ) - end - k₃ = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) * neginvdtγ - integrator.stats.nsolve += 1 - - if u isa Number - utilde = dto6 * f.mass_matrix[1, 1] * (k₁ - 2 * k₂ + k₃) - else - utilde = f.mass_matrix * (@.. dto6 * (k₁ - 2 * k₂ + k₃)) - end - atmp = calculate_residuals( - utilde, uprev, u, integrator.opts.abstol, - integrator.opts.reltol, integrator.opts.internalnorm, t - ) - integrator.EEst = integrator.opts.internalnorm(atmp, t) - - if mass_matrix !== I - invatol = inv(integrator.opts.abstol) - atmp = @. ifelse(integrator.differential_vars, false, integrator.fsallast) * - invatol - integrator.EEst += integrator.opts.internalnorm(atmp, t) - end - end - integrator.k[1] = k₁ - integrator.k[2] = k₂ - integrator.u = u - return nothing -end - -@muladd function perform_step!( - integrator, cache::Rosenbrock32ConstantCache, - repeat_step = false - ) - (; t, dt, uprev, u, f, p) = integrator - (; c₃₂, d, tf, uf) = cache - - # Precalculations - dtγ = dt * d - neginvdtγ = -inv(dtγ) - dto2 = dt / 2 - dto6 = dt / 6 - - mass_matrix = integrator.f.mass_matrix - - if repeat_step - integrator.fsalfirst = f(uprev, p, t) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - end - - # Time derivative and W matrix (with Jacobian reuse for W-methods) - dT, W = calc_rosenbrock_differentiation(integrator, cache, dtγ, repeat_step) - if !issuccess_W(W) - integrator.EEst = 2 - return nothing - end - - k₁ = _reshape(W \ -_vec((integrator.fsalfirst + dtγ * dT)), axes(uprev)) / dtγ - integrator.stats.nsolve += 1 - tmp = @.. uprev + dto2 * k₁ - f₁ = f(tmp, p, t + dto2) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - k₂ = _reshape(W \ _vec(f₁ - k₁), axes(uprev)) - else - linsolve_tmp = f₁ - mass_matrix * k₁ - k₂ = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) - end - k₂ = @.. k₂ * neginvdtγ + k₁ - - integrator.stats.nsolve += 1 - tmp = @.. uprev + dt * k₂ - integrator.fsallast = f(tmp, p, t + dt) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - if mass_matrix === I - linsolve_tmp = @.. ( - integrator.fsallast - c₃₂ * (k₂ - f₁) - - 2(k₁ - integrator.fsalfirst) + dt * dT - ) - else - linsolve_tmp = mass_matrix * (@.. c₃₂ * k₂ + 2 * k₁) - linsolve_tmp = @.. ( - integrator.fsallast - linsolve_tmp + - c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT - ) - end - k₃ = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) * neginvdtγ - integrator.stats.nsolve += 1 - u = @.. uprev + dto6 * (k₁ + 4k₂ + k₃) - - if integrator.opts.adaptive - utilde = @.. dto6 * (k₁ - 2k₂ + k₃) - atmp = calculate_residuals( - utilde, uprev, u, integrator.opts.abstol, - integrator.opts.reltol, integrator.opts.internalnorm, t - ) - integrator.EEst = integrator.opts.internalnorm(atmp, t) - - if mass_matrix !== I - invatol = inv(integrator.opts.abstol) - atmp = ifelse(integrator.differential_vars, false, integrator.fsallast) .* - invatol - integrator.EEst += integrator.opts.internalnorm(atmp, t) - end - end - - integrator.k[1] = k₁ - integrator.k[2] = k₂ - integrator.u = u - return nothing -end - #### Rodas4 type method — unified perform_step for all Rosenbrock/Rodas methods function initialize!(integrator, cache::RosenbrockCombinedConstantCache) @@ -569,10 +130,12 @@ end integrator.EEst = max(EEst, integrator.EEst) end else - # No H matrix: set k[1]=f(uprev), k[2]=f(u_new) for Hermite interpolation - integrator.k[1] = fsalfirst_cache - integrator.k[2] = f(u, p, t + dt) + # No H matrix: compute Hermite interpolation coefficients + # k[1] = dt*f₀ - (u - uprev), k[2] = 2*(u - uprev) - dt*(f₀ + f₁) + f1 = f(u, p, t + dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + integrator.k[1] = @.. dt * fsalfirst_cache - (u - uprev) + integrator.k[2] = @.. 2 * (u - uprev) - dt * (fsalfirst_cache + f1) end end @@ -731,11 +294,12 @@ end integrator.EEst = max(EEst, integrator.EEst) end else - # No H matrix: set k[1]=fsalfirst, k[2]=f(u_new) for interpolation - integrator.k[1] .= cache.fsalfirst + # No H matrix: compute Hermite interpolation coefficients + # k[1] = dt*f₀ - (u - uprev), k[2] = 2*(u - uprev) - dt*(f₀ + f₁) f(du, u, p, t + dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - integrator.k[2] .= du + @.. integrator.k[1] = dt * cache.fsalfirst - (u - uprev) + @.. integrator.k[2] = 2 * (u - uprev) - dt * (cache.fsalfirst + du) end end cache.linsolve = linres.cache diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl index 9489201e8fb..983de487590 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl @@ -1,23 +1,83 @@ -struct Rosenbrock23Tableau{T} - c₃₂::T - d::T -end +################################################################################ +# Rosenbrock23 / Rosenbrock32 (Shampine, 3-stage) +# +# Original Shampine formulation uses scaled stages k̃ᵢ with J*k coupling. +# Converted to Rodas form (C/dt coupling, no explicit J*k) via the identity +# J = W + M/(γdt) +# which absorbs the Jacobian coupling into the LHS operator. +# +# Stage variable relationship: k̃ᵢ_old = Σⱼ L[i,j]*ks[j]/(dt*γ) +# where L = [1 0 0; 1 1 0; 2 c₃₂ 1]. +################################################################################ -function Rosenbrock23Tableau(T) - c₃₂ = convert(T, 6 + sqrt(2)) - d = convert(T, 1 / (2 + sqrt(2))) - return Rosenbrock23Tableau(c₃₂, d) -end +function Rosenbrock23RodasTableau(T, T2) + gamma = convert(T2, 1 / (2 + sqrt(2))) + igamma = inv(gamma) + c32_old = convert(T, 6 + sqrt(2)) + + A = zeros(T, 3, 3) + A[2, 1] = convert(T, igamma / 2) # 1/(2γ) + A[3, 1] = convert(T, igamma) # 1/γ + A[3, 2] = convert(T, igamma) # 1/γ + + C = zeros(T, 3, 3) + C[2, 1] = convert(T, -igamma) # -1/γ + C[3, 1] = convert(T, -2 * igamma) # -2/γ + C[3, 2] = convert(T, -c32_old * igamma) # -c₃₂/γ + + c = T2[zero(T2), convert(T2, 1 // 2), one(T2)] + d = T[convert(T, gamma), zero(T), convert(T, 1 - 2 * gamma)] + + # Rosenbrock23: order-2 solution uses stages 1,2 only + b = T[convert(T, igamma), convert(T, igamma), zero(T)] + + # Error estimate: dt/6*(k̃₁ - 2k̃₂ + k̃₃) in new stage variables + btilde = T[ + convert(T, igamma / 6), + convert(T, (c32_old - 2) * igamma / 6), + convert(T, igamma / 6), + ] + + H = zeros(T, 0, 3) -struct Rosenbrock32Tableau{T} - c₃₂::T - d::T + return RodasTableau(A, C, gamma, c, d, H, b, btilde) end -function Rosenbrock32Tableau(T) - c₃₂ = convert(T, 6 + sqrt(2)) - d = convert(T, 1 / (2 + sqrt(2))) - return Rosenbrock32Tableau(c₃₂, d) +function Rosenbrock32RodasTableau(T, T2) + gamma = convert(T2, 1 / (2 + sqrt(2))) + igamma = inv(gamma) + c32_old = convert(T, 6 + sqrt(2)) + + A = zeros(T, 3, 3) + A[2, 1] = convert(T, igamma / 2) # 1/(2γ) + A[3, 1] = convert(T, igamma) # 1/γ + A[3, 2] = convert(T, igamma) # 1/γ + + C = zeros(T, 3, 3) + C[2, 1] = convert(T, -igamma) # -1/γ + C[3, 1] = convert(T, -2 * igamma) # -2/γ + C[3, 2] = convert(T, -c32_old * igamma) # -c₃₂/γ + + c = T2[zero(T2), convert(T2, 1 // 2), one(T2)] + d = T[convert(T, gamma), zero(T), convert(T, 1 - 2 * gamma)] + + # Rosenbrock32: order-3 solution dt/6*(k̃₁ + 4k̃₂ + k̃₃) in new stages + b = T[ + convert(T, 7 * igamma / 6), + convert(T, (4 + c32_old) * igamma / 6), + convert(T, igamma / 6), + ] + + # Error estimate (same as Rosenbrock23) + btilde = T[ + convert(T, igamma / 6), + convert(T, (c32_old - 2) * igamma / 6), + convert(T, igamma / 6), + ] + + H = zeros(T, 0, 3) + + return RodasTableau(A, C, gamma, c, d, H, b, btilde) end ################################################################################ diff --git a/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl b/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl index 07ed3fa1c41..a1db3975210 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl @@ -1,55 +1,3 @@ -function _ode_addsteps!( - k, t, uprev, u, dt, f, p, - cache::Union{ - Rosenbrock23ConstantCache, - Rosenbrock32ConstantCache, - }, - always_calc_begin = false, allow_calc_end = true, - force_calc_end = false - ) - if length(k) < 2 || always_calc_begin - (; tf, uf, d) = cache - dtγ = dt * d - neginvdtγ = -inv(dtγ) - dto2 = dt / 2 - tf.u = uprev - - autodiff_alg = cache.autodiff - - if autodiff_alg isa AutoFiniteDiff - autodiff_alg = SciMLBase.@set autodiff_alg.dir = sign(dt) - end - - dT = DI.derivative(tf, autodiff_alg, t) - - mass_matrix = f.mass_matrix - if uprev isa Number - J = DI.derivative(uf, autodiff_alg, uprev) - W = neginvdtγ .+ J - else - J = DI.jacobian(uf, autodiff_alg, uprev) - if mass_matrix isa UniformScaling - W = neginvdtγ * mass_matrix + J - else - W = @.. neginvdtγ * mass_matrix .+ J - end - end - f₀ = f(uprev, p, t) - k₁ = _reshape(W \ _vec((f₀ + dtγ * dT)), axes(uprev)) * neginvdtγ - tmp = @.. uprev + dto2 * k₁ - f₁ = f(tmp, p, t + dto2) - if mass_matrix === I - k₂ = _reshape(W \ _vec(f₁ - k₁), axes(uprev)) - else - k₂ = _reshape(W \ _vec(f₁ - mass_matrix * k₁), axes(uprev)) - end - k₂ = @.. k₂ * neginvdtγ + k₁ - copyat_or_push!(k, 1, k₁) - copyat_or_push!(k, 2, k₂) - end - return nothing -end - function _ode_addsteps!( k, t, uprev, u, dt, f, p, cache::RosenbrockCombinedConstantCache, always_calc_begin = false, allow_calc_end = true, From f458907303c5255aae25135b084d808151ed7767 Mon Sep 17 00:00:00 2001 From: Harsh Singh Date: Wed, 15 Apr 2026 18:15:24 +0530 Subject: [PATCH 2/6] fix: strip_cache + resize test field names for RosenbrockCache --- .../src/interp_func.jl | 28 +++++++++++++++++++ test/integrators/resize_tests.jl | 18 ++++-------- 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl b/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl index 314f2792aa9..eed96fd916d 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/interp_func.jl @@ -27,3 +27,31 @@ function SciMLBase.interp_summary( "specialized 4th (Rodas6P = 5th) order \"free\" stiffness-aware interpolation" : "1st order linear" end + +# strip_cache for RosenbrockCache: the generic OrdinaryDiffEqCore version passes +# all-Nothing args, but several fields (dense, dtC, dtd, ks, interp_order, jac_reuse) have +# concrete types that don't accept Nothing. Provide a custom override that clears +# only the interpolation-related fields to zero-length arrays / nothing. +function OrdinaryDiffEqCore.strip_cache(cache::RosenbrockCache) + return RosenbrockCache( + cache.u, cache.uprev, + similar(cache.dense, 0), # dense::Vector{rateType} — must not be Nothing + cache.du, cache.du1, cache.du2, + similar(cache.dtC, 0, 0), # dtC::Matrix{tabType} — must not be Nothing + similar(cache.dtd, 0), # dtd::Vector{tabType} — must not be Nothing + similar(cache.ks, 0), # ks::Vector{rateType} — must not be Nothing + cache.fsalfirst, cache.fsallast, cache.dT, + cache.J, cache.W, + cache.tmp, cache.atmp, cache.weight, + cache.tab, + nothing, # tf + nothing, # uf + cache.linsolve_tmp, cache.linsolve, + nothing, # jac_config + nothing, # grad_config + cache.reltol, cache.alg, + cache.step_limiter!, cache.stage_limiter!, + cache.interp_order, # interp_order::Int — must not be Nothing + cache.jac_reuse # jac_reuse — preserve state across strip + ) +end diff --git a/test/integrators/resize_tests.jl b/test/integrators/resize_tests.jl index 68483eec8e3..eb0b888e5ac 100644 --- a/test/integrators/resize_tests.jl +++ b/test/integrators/resize_tests.jl @@ -76,12 +76,10 @@ i = init(prob, Rosenbrock23()) resize!(i, 5) @test length(i.cache.u) == 5 @test length(i.cache.uprev) == 5 -@test length(i.cache.k₁) == 5 -@test length(i.cache.k₂) == 5 -@test length(i.cache.k₃) == 5 +@test length(i.cache.ks) == 3 +@test all(length.(i.cache.ks) .== 5) @test length(i.cache.du1) == 5 @test length(i.cache.du2) == 5 -@test length(i.cache.f₁) == 5 @test length(i.cache.fsalfirst) == 5 @test length(i.cache.fsallast) == 5 @test length(i.cache.dT) == 5 @@ -104,12 +102,10 @@ i = init(prob, Rosenbrock23(autodiff = AutoForwardDiff(), linsolve = KrylovJL_GM resize!(i, 5) @test length(i.cache.u) == 5 @test length(i.cache.uprev) == 5 -@test length(i.cache.k₁) == 5 -@test length(i.cache.k₂) == 5 -@test length(i.cache.k₃) == 5 +@test length(i.cache.ks) == 3 +@test all(length.(i.cache.ks) .== 5) @test length(i.cache.du1) == 5 @test length(i.cache.du2) == 5 -@test length(i.cache.f₁) == 5 @test length(i.cache.fsalfirst) == 5 @test length(i.cache.fsallast) == 5 @test length(i.cache.dT) == 5 @@ -123,12 +119,10 @@ i = init(prob, Rosenbrock23(; autodiff = AutoFiniteDiff())) resize!(i, 5) @test length(i.cache.u) == 5 @test length(i.cache.uprev) == 5 -@test length(i.cache.k₁) == 5 -@test length(i.cache.k₂) == 5 -@test length(i.cache.k₃) == 5 +@test length(i.cache.ks) == 3 +@test all(length.(i.cache.ks) .== 5) @test length(i.cache.du1) == 5 @test length(i.cache.du2) == 5 -@test length(i.cache.f₁) == 5 @test length(i.cache.fsalfirst) == 5 @test length(i.cache.fsallast) == 5 @test length(i.cache.dT) == 5 From f513f8031174df0f48dd86740dfb23694c2f0d2b Mon Sep 17 00:00:00 2001 From: Harsh Singh Date: Thu, 16 Apr 2026 00:46:14 +0530 Subject: [PATCH 3/6] fix: OOP addsteps guard for H=0 methods (Rosenbrock23/32) --- lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl b/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl index a1db3975210..760832619ca 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/stiff_addsteps.jl @@ -3,7 +3,9 @@ function _ode_addsteps!( always_calc_begin = false, allow_calc_end = true, force_calc_end = false ) - if length(k) < size(cache.tab.H, 1) || always_calc_begin + H_rows = size(cache.tab.H, 1) + kneeded = H_rows == 0 ? 2 : H_rows + if length(k) < kneeded || always_calc_begin (; tf, uf) = cache (; A, C, gamma, c, d, H) = cache.tab From f71e15e549e7a27c9cd86ee1384b71ca02871b86 Mon Sep 17 00:00:00 2001 From: Harsh Singh Date: Thu, 16 Apr 2026 01:18:58 +0530 Subject: [PATCH 4/6] fix: restore Shampine stiff interpolation via 2-row H matrix --- .../src/rosenbrock_tableaus.jl | 14 ++++++++++++-- test/integrators/resize_tests.jl | 3 +++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl index 983de487590..fed89efccff 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl @@ -38,7 +38,15 @@ function Rosenbrock23RodasTableau(T, T2) convert(T, igamma / 6), ] - H = zeros(T, 0, 3) + # Shampine 2nd-order "free" stiff interpolation encoded as a 2-row H matrix + # (interp_order=2) so the standard Rodas 2-vector formula applies: + # p(Θ) = (1-Θ)y₀ + Θ(y₁ + (1-Θ)(K[1] + Θ·K[2])) + # K[1] = -ks[2]/(γ(1-2γ)), K[2] = 0 + # Algebraically equals the MATLAB ODE Suite c1·k̃₁ + c2·k̃₂ dense output + # with c1=Θ(1-Θ)/(1-2γ), c2=Θ(Θ-2γ)/(1-2γ); K[2]=0 because the Θ² term + # from c2 is already absorbed into K[1] via the y₁ endpoint. + H = T[zero(T) convert(T, -1/(gamma*(1 - 2*gamma))) zero(T) + zero(T) zero(T) zero(T)] return RodasTableau(A, C, gamma, c, d, H, b, btilde) end @@ -75,7 +83,9 @@ function Rosenbrock32RodasTableau(T, T2) convert(T, igamma / 6), ] - H = zeros(T, 0, 3) + # Same Shampine interpolation as Rosenbrock23 (identical γ, same stencil) + H = T[zero(T) convert(T, -1/(gamma*(1 - 2*gamma))) zero(T) + zero(T) zero(T) zero(T)] return RodasTableau(A, C, gamma, c, d, H, b, btilde) end diff --git a/test/integrators/resize_tests.jl b/test/integrators/resize_tests.jl index eb0b888e5ac..acc3d8c48f8 100644 --- a/test/integrators/resize_tests.jl +++ b/test/integrators/resize_tests.jl @@ -78,6 +78,7 @@ resize!(i, 5) @test length(i.cache.uprev) == 5 @test length(i.cache.ks) == 3 @test all(length.(i.cache.ks) .== 5) +@test all(length.(i.cache.dense) .== 5) @test length(i.cache.du1) == 5 @test length(i.cache.du2) == 5 @test length(i.cache.fsalfirst) == 5 @@ -104,6 +105,7 @@ resize!(i, 5) @test length(i.cache.uprev) == 5 @test length(i.cache.ks) == 3 @test all(length.(i.cache.ks) .== 5) +@test all(length.(i.cache.dense) .== 5) @test length(i.cache.du1) == 5 @test length(i.cache.du2) == 5 @test length(i.cache.fsalfirst) == 5 @@ -121,6 +123,7 @@ resize!(i, 5) @test length(i.cache.uprev) == 5 @test length(i.cache.ks) == 3 @test all(length.(i.cache.ks) .== 5) +@test all(length.(i.cache.dense) .== 5) @test length(i.cache.du1) == 5 @test length(i.cache.du2) == 5 @test length(i.cache.fsalfirst) == 5 From d2e864b42b924948854fd5d3a31bfd0d6755be61 Mon Sep 17 00:00:00 2001 From: Harsh Singh Date: Fri, 17 Apr 2026 01:36:16 +0530 Subject: [PATCH 5/6] comment fixes --- .../src/rosenbrock_perform_step.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl index bce48e7a770..4a2556a0412 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl @@ -130,12 +130,12 @@ end integrator.EEst = max(EEst, integrator.EEst) end else - # No H matrix: compute Hermite interpolation coefficients - # k[1] = dt*f₀ - (u - uprev), k[2] = 2*(u - uprev) - dt*(f₀ + f₁) + # No H matrix: store raw derivatives f₀, f₁ for standard Hermite interpolation. + # _ode_interpolant(interp_order=-1) expects k[1]=f₀, k[2]=f₁. f1 = f(u, p, t + dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - integrator.k[1] = @.. dt * fsalfirst_cache - (u - uprev) - integrator.k[2] = @.. 2 * (u - uprev) - dt * (fsalfirst_cache + f1) + integrator.k[1] = fsalfirst_cache # f₀ = f(uprev, p, t) + integrator.k[2] = f1 # f₁ = f(u, p, t+dt) end end @@ -294,12 +294,12 @@ end integrator.EEst = max(EEst, integrator.EEst) end else - # No H matrix: compute Hermite interpolation coefficients - # k[1] = dt*f₀ - (u - uprev), k[2] = 2*(u - uprev) - dt*(f₀ + f₁) + # No H matrix: store raw derivatives f₀, f₁ for standard Hermite interpolation. + # _ode_interpolant(interp_order=-1) expects k[1]=f₀, k[2]=f₁. f(du, u, p, t + dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. integrator.k[1] = dt * cache.fsalfirst - (u - uprev) - @.. integrator.k[2] = 2 * (u - uprev) - dt * (cache.fsalfirst + du) + @.. integrator.k[1] = cache.fsalfirst # f₀ = f(uprev, p, t) + @.. integrator.k[2] = du # f₁ = f(u, p, t+dt) end end cache.linsolve = linres.cache From 1e70c8b7e3412b627f2089e495fee6e80da61264 Mon Sep 17 00:00:00 2001 From: Harsh Singh Date: Fri, 17 Apr 2026 03:03:24 +0530 Subject: [PATCH 6/6] fix: runic formatting for H matrix literals in Rosenbrock23/32 tableaus --- .../src/rosenbrock_tableaus.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl index fed89efccff..6ea68ad3dec 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl @@ -45,8 +45,10 @@ function Rosenbrock23RodasTableau(T, T2) # Algebraically equals the MATLAB ODE Suite c1·k̃₁ + c2·k̃₂ dense output # with c1=Θ(1-Θ)/(1-2γ), c2=Θ(Θ-2γ)/(1-2γ); K[2]=0 because the Θ² term # from c2 is already absorbed into K[1] via the y₁ endpoint. - H = T[zero(T) convert(T, -1/(gamma*(1 - 2*gamma))) zero(T) - zero(T) zero(T) zero(T)] + H = T[ + zero(T) convert(T, -1 / (gamma * (1 - 2 * gamma))) zero(T) + zero(T) zero(T) zero(T) + ] return RodasTableau(A, C, gamma, c, d, H, b, btilde) end @@ -84,8 +86,10 @@ function Rosenbrock32RodasTableau(T, T2) ] # Same Shampine interpolation as Rosenbrock23 (identical γ, same stencil) - H = T[zero(T) convert(T, -1/(gamma*(1 - 2*gamma))) zero(T) - zero(T) zero(T) zero(T)] + H = T[ + zero(T) convert(T, -1 / (gamma * (1 - 2 * gamma))) zero(T) + zero(T) zero(T) zero(T) + ] return RodasTableau(A, C, gamma, c, d, H, b, btilde) end