From 4a6bfc8f8003514d77f0c6f4afdf744e51dd54d3 Mon Sep 17 00:00:00 2001 From: singhharsh1708 Date: Fri, 31 Jul 2026 02:26:40 +0530 Subject: [PATCH] Fold RK46NL into the generic LowStorageRK2N dispatch RK46NL was the last plain-2N method still carrying a bespoke constant cache struct, mutable cache struct, and four hand-unrolled initialize!/perform_step! methods. Its recurrence is the same Williamson form the shared LowStorageRK2N kernels already run, so the tableau moves into a RK46NLConstantCache function returning a LowStorageRK2NConstantCache with unchanged coefficients, RK46NL joins the LowStorageRK2NAlgorithm union, and the bespoke caches and kernels are deleted. The struct gains the williamson_condition keyword (default false) that the shared alg_cache reads, and isfsal/uses_uprev are declared false to match the rest of the family. The generic in-place kernel re-evaluates f at the start of each step instead of using FSAL storage, keeping the per-step f-call count at six. The new field is fourth, so positional construction now takes four arguments, the same break the nine sibling 2N methods already took when they gained the keyword. Fixed-dt end states are bit-identical to the removed implementation for in-place problems and for out-of-place problems whose state is an array. Out-of-place problems with a scalar state differ in the last bits: 3 ulp after 128 steps of u' = 1.01u, 1 ulp after 256 steps of a nonlinear scalar problem, roughly four thousand times below the method's own discretization error there. Both kernels expand to the same muladd(dt, k, A*tmp), but the straight-line code folds the other product into the fused multiply-add, and only scalar states are affected because muladd over arrays never fuses. RK46NL gains the williamson_condition testset the nine siblings already have. Version goes to 3.4.0 for the new public keyword. --- lib/OrdinaryDiffEqLowStorageRK/Project.toml | 2 +- .../src/alg_utils.jl | 2 + .../src/algorithms.jl | 6 +- .../src/low_storage_rk_caches.jl | 96 +++++-------------- .../src/low_storage_rk_perform_step.jl | 90 ----------------- .../test/ode_low_storage_rk_tests.jl | 58 +++++++++++ 6 files changed, 89 insertions(+), 165 deletions(-) diff --git a/lib/OrdinaryDiffEqLowStorageRK/Project.toml b/lib/OrdinaryDiffEqLowStorageRK/Project.toml index 5e86e1345ea..4c0497a761d 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/Project.toml +++ b/lib/OrdinaryDiffEqLowStorageRK/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqLowStorageRK" uuid = "b0944070-b475-4768-8dec-fb6eb410534d" authors = ["ParamThakkar123 "] -version = "3.3.0" +version = "3.4.0" [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl index 0495abd9280..e437c6c346f 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl @@ -60,6 +60,7 @@ alg_order(alg::SHLDDRK_2N) = 4 alg_order(alg::SHLDDRK52) = 2 isfsal(alg::ORK256) = false +isfsal(alg::RK46NL) = false isfsal(alg::CarpenterKennedy2N54) = false isfsal(alg::DGLDDRK84_F) = false isfsal(alg::DGLDDRK73_C) = false @@ -73,6 +74,7 @@ isfsal(alg::NDBLSRK124) = false isfsal(alg::NDBLSRK144) = false uses_uprev(alg::ORK256, adaptive::Bool) = false +uses_uprev(alg::RK46NL, adaptive::Bool) = false uses_uprev(alg::SHLDDRK64, adaptive::Bool) = false uses_uprev(alg::CarpenterKennedy2N54, adaptive::Bool) = false uses_uprev(alg::NDBLSRK124, adaptive::Bool) = false diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl b/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl index 48a84d82741..a2946c6635f 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl @@ -195,12 +195,16 @@ end @doc explicit_rk_docstring( "6-stage, fourth order low-stage, low-dissipation, low-dispersion scheme. Fixed timestep only.", "RK46NL", - references = "Julien Berland, Christophe Bogey, Christophe Bailly. Low-Dissipation and Low-Dispersion Fourth-Order Runge-Kutta Algorithm. Computers & Fluids, 35(10), pp 1459-1463, 2006. doi: https://doi.org/10.1016/j.compfluid.2005.04.003" + references = "Julien Berland, Christophe Bogey, Christophe Bailly. Low-Dissipation and Low-Dispersion Fourth-Order Runge-Kutta Algorithm. Computers & Fluids, 35(10), pp 1459-1463, 2006. doi: https://doi.org/10.1016/j.compfluid.2005.04.003", + extra_keyword_description = """- `williamson_condition`: allows for an optimization that allows fusing broadcast expressions with the function call `f`. However, it only works for `Array` types. + """, + extra_keyword_default = "williamson_condition = false" ) Base.@kwdef struct RK46NL{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm stage_limiter!::StageLimiter = trivial_limiter! step_limiter!::StepLimiter = trivial_limiter! thread::Thread = Serial() + williamson_condition::Bool = false end @doc explicit_rk_docstring( diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl index 6a436b6ce67..0c71fb3aeef 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl @@ -47,81 +47,30 @@ function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2} return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end) end -function alg_cache( - alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}, verbose - ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - return RK46NLConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) -end +function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T, T2} + A2 = convert(T, -0.737101392796) + A3 = convert(T, -1.634740794343) + A4 = convert(T, -0.74473900378) + A5 = convert(T, -1.469897351522) + A6 = convert(T, -2.813971388035) + A2end = (A2, A3, A4, A5, A6) -struct RK46NLConstantCache{T, T2} <: OrdinaryDiffEqConstantCache - α2::T - α3::T - α4::T - α5::T - α6::T - β1::T - β2::T - β3::T - β4::T - β5::T - β6::T - c2::T2 - c3::T2 - c4::T2 - c5::T2 - c6::T2 - - function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T, T2} - α2 = T(-0.737101392796) - α3 = T(-1.634740794343) - α4 = T(-0.74473900378) - α5 = T(-1.469897351522) - α6 = T(-2.813971388035) - β1 = T(0.032918605146) - β2 = T(0.8232569982) - β3 = T(0.3815309489) - β4 = T(0.200092213184) - β5 = T(1.718581042715) - β6 = T(0.27) - c2 = T2(0.032918605146) - c3 = T2(0.249351723343) - c4 = T2(0.466911705055) - c5 = T2(0.582030414044) - c6 = T2(0.847252983783) - return new{T, T2}(α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) - end -end + B1 = convert(T, 0.032918605146) + B2 = convert(T, 0.8232569982) + B3 = convert(T, 0.3815309489) + B4 = convert(T, 0.200092213184) + B5 = convert(T, 1.718581042715) + B6 = convert(T, 0.27) + B2end = (B2, B3, B4, B5, B6) -function alg_cache( - alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}, verbose - ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - tmp = zero(u) - k = zero(rate_prototype) - fsalfirst = zero(rate_prototype) - tab = RK46NLConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) - return RK46NLCache( - u, uprev, k, tmp, fsalfirst, tab, alg.stage_limiter!, alg.step_limiter!, - alg.thread - ) -end + c2 = convert(T2, 0.032918605146) + c3 = convert(T2, 0.249351723343) + c4 = convert(T2, 0.466911705055) + c5 = convert(T2, 0.582030414044) + c6 = convert(T2, 0.847252983783) + c2end = (c2, c3, c4, c5, c6) -@cache struct RK46NLCache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <: - LowStorageRKMutableCache - u::uType - uprev::uType - k::rateType - tmp::uType - fsalfirst::rateType - tab::TabType - stage_limiter!::StageLimiter - step_limiter!::StepLimiter - thread::Thread + return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end) end function CarpenterKennedy2N54ConstantCache(::Type{T}, ::Type{T2}) where {T, T2} @@ -4357,11 +4306,12 @@ end # --- LowStorageRK2N --- const LowStorageRK2NAlgorithm = Union{ - ORK256, CarpenterKennedy2N54, SHLDDRK64, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F, + ORK256, RK46NL, CarpenterKennedy2N54, SHLDDRK64, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F, NDBLSRK124, NDBLSRK134, NDBLSRK144, } _lowstorage_2n_tableau(::ORK256, ::Type{T}, ::Type{T2}) where {T, T2} = ORK256ConstantCache(T, T2) +_lowstorage_2n_tableau(::RK46NL, ::Type{T}, ::Type{T2}) where {T, T2} = RK46NLConstantCache(T, T2) _lowstorage_2n_tableau(::CarpenterKennedy2N54, ::Type{T}, ::Type{T2}) where {T, T2} = CarpenterKennedy2N54ConstantCache(T, T2) _lowstorage_2n_tableau(::SHLDDRK64, ::Type{T}, ::Type{T2}) where {T, T2} = SHLDDRK64ConstantCache(T, T2) _lowstorage_2n_tableau(::DGLDDRK73_C, ::Type{T}, ::Type{T2}) where {T, T2} = DGLDDRK73_CConstantCache(T, T2) diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl index 24ef9f39fa4..21cf4e024bc 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl @@ -903,96 +903,6 @@ end OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) end -function initialize!(integrator, cache::RK46NLCache) - (; k, fsalfirst) = cache - - integrator.kshortsize = 1 - resize!(integrator.k, integrator.kshortsize) - integrator.k[1] = integrator.fsalfirst - integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # FSAL for interpolation - return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) -end - -@muladd function perform_step!(integrator, cache::RK46NLCache, repeat_step = false) - (; t, dt, uprev, u, f, p) = integrator - (; k, fsalfirst, tmp, thread) = cache - stage_limiter! = integrator.opts.stage_limiter! - (; α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) = cache.tab - - # u1 - @.. broadcast = false thread = thread tmp = dt * fsalfirst - @.. broadcast = false thread = thread u = uprev + β1 * tmp - stage_limiter!(u, integrator, p, t + c2 * dt) - # u2 - f(k, u, p, t + c2 * dt) - @.. broadcast = false thread = thread tmp = α2 * tmp + dt * k - @.. broadcast = false thread = thread u = u + β2 * tmp - stage_limiter!(u, integrator, p, t + c3 * dt) - # u3 - f(k, u, p, t + c3 * dt) - @.. broadcast = false thread = thread tmp = α3 * tmp + dt * k - @.. broadcast = false thread = thread u = u + β3 * tmp - stage_limiter!(u, integrator, p, t + c4 * dt) - # u4 - f(k, u, p, t + c4 * dt) - @.. broadcast = false thread = thread tmp = α4 * tmp + dt * k - @.. broadcast = false thread = thread u = u + β4 * tmp - stage_limiter!(u, integrator, p, t + c5 * dt) - # u5 = u - f(k, u, p, t + c5 * dt) - @.. broadcast = false thread = thread tmp = α5 * tmp + dt * k - @.. broadcast = false thread = thread u = u + β5 * tmp - stage_limiter!(u, integrator, p, t + c6 * dt) - - f(k, u, p, t + c6 * dt) - @.. broadcast = false thread = thread tmp = α6 * tmp + dt * k - @.. broadcast = false thread = thread u = u + β6 * tmp - stage_limiter!(u, integrator, p, t + dt) - - f(k, u, p, t + dt) - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 6) -end - -function initialize!(integrator, cache::RK46NLConstantCache) - integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - integrator.kshortsize = 1 - integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) - - # Avoid undefined entries if k is an array of arrays - integrator.fsallast = zero(integrator.fsalfirst) - return integrator.k[1] = integrator.fsalfirst -end - -@muladd function perform_step!(integrator, cache::RK46NLConstantCache, repeat_step = false) - (; t, dt, uprev, u, f, p) = integrator - (; α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) = cache - - # u1 - tmp = dt * integrator.fsalfirst - u = uprev + β1 * tmp - # u2 - tmp = α2 * tmp + dt * f(u, p, t + c2 * dt) - u = u + β2 * tmp - # u3 - tmp = α3 * tmp + dt * f(u, p, t + c3 * dt) - u = u + β3 * tmp - # u4 - tmp = α4 * tmp + dt * f(u, p, t + c4 * dt) - u = u + β4 * tmp - # u5 = u - tmp = α5 * tmp + dt * f(u, p, t + c5 * dt) - u = u + β5 * tmp - # u6 - tmp = α6 * tmp + dt * f(u, p, t + c6 * dt) - u = u + β6 * tmp - - integrator.fsallast = f(u, p, t + dt) # For interpolation, then FSAL'd - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 6) - integrator.k[1] = integrator.fsalfirst - integrator.u = u -end - function initialize!(integrator, cache::SHLDDRK52ConstantCache) integrator.kshortsize = 2 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) diff --git a/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl b/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl index 27724b5889e..df998c58daa 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl @@ -143,6 +143,64 @@ end @test sol_old.u[end] ≈ sol_new.u[end] end +@testset "RK46NL" begin + alg = RK46NL(; williamson_condition = true) + alg2 = RK46NL(; williamson_condition = false) + dts = 1 ./ 2 .^ (7:-1:3) + for prob in test_problems_only_time + sim = test_convergence(dts, prob, alg) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol + sim = test_convergence(dts, prob, alg2) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol + end + for prob in test_problems_linear + sim = test_convergence(dts, prob, alg) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol + sim = test_convergence(dts, prob, alg2) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol + end + for prob in test_problems_nonlinear + sim = test_convergence(dts, prob, alg) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol + sim = test_convergence(dts, prob, alg2) + @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol + end + integ = init( + prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false, + save_everystep = false + ) + @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3 + integ = init( + prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false, + save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true) + ) + @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2 + integ = init( + prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false, + save_everystep = false + ) + @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4 + integ = init( + prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false, + save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true) + ) + @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3 + # test whether aliasing u0 is bad + new_prob_ode_nonlinear_inplace = ODEProblem( + prob_ode_nonlinear_inplace.f, [1.0], + (0.0, 0.5) + ) + sol_old = solve( + prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false, + save_start = false + ) + sol_new = solve( + new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false, + save_start = false, alias = ODEAliasSpecifier(alias_u0 = true) + ) + @test sol_old.u[end] ≈ sol_new.u[end] +end + @testset "CarpenterKennedy2N54" begin alg = CarpenterKennedy2N54(; williamson_condition = true) alg2 = CarpenterKennedy2N54(; williamson_condition = false)