Skip to content

Fold RK46NL into the generic LowStorageRK2N dispatch #3683

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ChrisRackauckas merged 1 commit into from
Aug 1, 2026
Merged

Fold RK46NL into the generic LowStorageRK2N dispatch #3683

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion lib/OrdinaryDiffEqLowStorageRK/Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "OrdinaryDiffEqLowStorageRK"
uuid = "b0944070-b475-4768-8dec-fb6eb410534d"
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
version = "3.3.0"
version = "3.4.0"

[deps]
CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
Expand Down
2 changes: 2 additions & 0 deletions lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
6 changes: 5 additions & 1 deletion lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
96 changes: 23 additions & 73 deletions lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,81 +47,30 @@ function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)

@ChrisRackauckas ChrisRackauckas May 27, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

these all already hit the same code, so it's not a substantive change. Do any other caches hit the same form?

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}
Expand Down Expand Up @@ -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)
Expand Down
90 changes: 0 additions & 90 deletions lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
58 changes: 58 additions & 0 deletions lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading