Skip to content

Lowstorage 2c #3684

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

Closed
singhharsh1708 wants to merge 4 commits into from
Closed

Lowstorage 2c #3684

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
eedc675
LowStorageRK: unify 2N family on LowStorageRKTableau
May 26, 2026
0eed72f
LowStorageRK: use @enum LowStorageRKForm instead of Symbol for tablea…
May 27, 2026
05d30e1
LowStorageRK: unify 2C family on LowStorageRKTableau
May 27, 2026
ccc511a
LowStorageRK: use TwoC enum value for 2C tableau dispatch
May 27, 2026
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: 2 additions & 0 deletions lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,10 @@ using Reexport
include("arrayfuse.jl")
include("algorithms.jl")
include("alg_utils.jl")
include("low_storage_tableaus.jl")
include("low_storage_rk_caches.jl")
include("low_storage_rk_perform_step.jl")
include("generic_low_storage_perform_step.jl")

import PrecompileTools
import Preferences
Expand Down
95 changes: 95 additions & 0 deletions lib/OrdinaryDiffEqLowStorageRK/src/generic_low_storage_perform_step.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
@muladd function _perform_step_iip!(
integrator, cache, tab::LowStorageRKTableau{TwoN}
)
(; t, dt, u, f, p) = integrator
(; k, tmp, williamson_condition, stage_limiter!, step_limiter!, thread) = cache
(; A2end, B1, B2end, c2end) = tab

f(k, u, p, t)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
@.. broadcast = false thread = thread tmp = dt * k
@.. broadcast = false thread = thread u = u + B1 * tmp

for i in eachindex(A2end)
if williamson_condition
f(ArrayFuse(tmp, u, (A2end[i], dt, B2end[i])), u, p, t + c2end[i] * dt)
else
@.. broadcast = false thread = thread tmp = A2end[i] * tmp
stage_limiter!(u, integrator, p, t + c2end[i] * dt)
f(k, u, p, t + c2end[i] * dt)
@.. broadcast = false thread = thread tmp = tmp + dt * k
@.. broadcast = false thread = thread u = u + B2end[i] * tmp
end
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
end
stage_limiter!(u, integrator, p, t + dt)
step_limiter!(u, integrator, p, t + dt)
return nothing
end

@muladd function _perform_step_oop!(
integrator, tab::LowStorageRKTableau{TwoN}
)
(; t, dt, u, f, p) = integrator
(; A2end, B1, B2end, c2end) = tab

tmp = dt * integrator.fsalfirst
u = u + B1 * tmp

for i in eachindex(A2end)
k = f(u, p, t + c2end[i] * dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
tmp = A2end[i] * tmp + dt * k
u = u + B2end[i] * tmp
end

integrator.k[1] = integrator.fsalfirst
integrator.fsalfirst = f(u, p, t + dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
integrator.u = u
return nothing
end

@muladd function _perform_step_iip!(
integrator, cache, tab::LowStorageRKTableau{TwoC}
)
(; t, dt, u, f, p) = integrator
(; k, fsalfirst, tmp, stage_limiter!, step_limiter!, thread) = cache
(; A2end, B1, B2end, c2end) = tab

@.. broadcast = false thread = thread k = integrator.fsalfirst
@.. broadcast = false thread = thread u = u + B1 * dt * k

for i in eachindex(A2end)
@.. broadcast = false thread = thread tmp = u + A2end[i] * dt * k
f(k, tmp, p, t + c2end[i] * dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
@.. broadcast = false thread = thread u = u + B2end[i] * dt * k
end
step_limiter!(u, integrator, p, t + dt)
f(k, u, p, t + dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
return nothing
end

@muladd function _perform_step_oop!(
integrator, tab::LowStorageRKTableau{TwoC}
)
(; t, dt, u, f, p) = integrator
(; A2end, B1, B2end, c2end) = tab

k = integrator.fsalfirst = integrator.f(u, p, t)
integrator.k[1] = integrator.fsalfirst
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
u = u + B1 * dt * k

for i in eachindex(A2end)
tmp = u + A2end[i] * dt * k
k = integrator.f(tmp, p, t + c2end[i] * dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
u = u + B2end[i] * dt * k
end

integrator.u = u
return nothing
end
36 changes: 11 additions & 25 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 @@ -17,13 +17,6 @@ get_fsalfirstlast(cache::LowStorageRKMutableCache, u) = (cache.fsalfirst, cache.
thread::Thread
end

struct LowStorageRK2NConstantCache{N, T, T2} <: OrdinaryDiffEqConstantCache
A2end::NTuple{N, T} # A1 is always zero
B1::T
B2end::NTuple{N, T}
c2end::NTuple{N, T2} # c1 is always zero
end

function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
A2 = convert(T, -1.0)
A3 = convert(T, -1.55798)
Expand All @@ -44,7 +37,7 @@ function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c5 = convert(T2, 0.8)
c2end = (c2, c3, c4, c5)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

function alg_cache(
Expand Down Expand Up @@ -144,7 +137,7 @@ function CarpenterKennedy2N54ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c5 = convert(T2, 2802321613138 // 2924317926251)
c2end = (c2, c3, c4, c5)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

@cache mutable struct SHLDDRK_2NCache{
Expand Down Expand Up @@ -359,7 +352,7 @@ function SHLDDRK64ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c6 = convert(T2, 0.9271047)
c2end = (c2, c3, c4, c5, c6)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

function DGLDDRK73_CConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
Expand Down Expand Up @@ -388,7 +381,7 @@ function DGLDDRK73_CConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c7 = convert(T2, 0.998046608462379)
c2end = (c2, c3, c4, c5, c6, c7)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

function DGLDDRK84_CConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
Expand Down Expand Up @@ -420,7 +413,7 @@ function DGLDDRK84_CConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c8 = convert(T2, 0.919902896453866)
c2end = (c2, c3, c4, c5, c6, c7, c8)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

function DGLDDRK84_FConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
Expand Down Expand Up @@ -452,7 +445,7 @@ function DGLDDRK84_FConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c8 = convert(T2, 0.9484087623348481)
c2end = (c2, c3, c4, c5, c6, c7, c8)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

function NDBLSRK124ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
Expand Down Expand Up @@ -496,7 +489,7 @@ function NDBLSRK124ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c12 = convert(T2, 0.9032588871651854)
c2end = (c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

function NDBLSRK134ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
Expand Down Expand Up @@ -543,7 +536,7 @@ function NDBLSRK134ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c13 = convert(T2, 0.9126827615920843)
c2end = (c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

function NDBLSRK144ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
Expand Down Expand Up @@ -593,7 +586,7 @@ function NDBLSRK144ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c14 = convert(T2, 0.8734213127600976)
c2end = (c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14)

return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoN}(A2end, B1, B2end, c2end)
end

# 2C low storage methods introduced by Calvo, Franco, Rández (2004)
Expand All @@ -612,13 +605,6 @@ end
thread::Thread
end

struct LowStorageRK2CConstantCache{N, T, T2} <: OrdinaryDiffEqConstantCache
A2end::NTuple{N, T} # A1 is always zero
B1::T
B2end::NTuple{N, T}
c2end::NTuple{N, T2} # c1 is always zero
end

function CFRLDDRK64ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
A2 = convert(T, 0.17985400977138)
A3 = convert(T, 0.14081893152111)
Expand All @@ -642,7 +628,7 @@ function CFRLDDRK64ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c6 = convert(T2, 0.83050587987157)
c2end = (c2, c3, c4, c5, c6)

return LowStorageRK2CConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoC}(A2end, B1, B2end, c2end)
end

function TSLDDRK74ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
Expand Down Expand Up @@ -671,7 +657,7 @@ function TSLDDRK74ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
c7 = convert(T2, 0.91124223849547205)
c2end = (c2, c3, c4, c5, c6, c7)

return LowStorageRK2CConstantCache(A2end, B1, B2end, c2end)
return LowStorageRKTableau{TwoC}(A2end, B1, B2end, c2end)
end

# 3S low storage methods introduced by Ketcheson
Expand Down
117 changes: 14 additions & 103 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
@@ -1,38 +1,14 @@
# 2N low storage methods
function initialize!(integrator, cache::LowStorageRK2NConstantCache)
integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
function initialize!(integrator, cache::LowStorageRKTableau{TwoN})
integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
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::LowStorageRK2NConstantCache,
repeat_step = false
)
(; t, dt, u, f, p) = integrator
(; A2end, B1, B2end, c2end) = cache

# u1
tmp = dt * integrator.fsalfirst
u = u + B1 * tmp

# other stages
for i in eachindex(A2end)
k = f(u, p, t + c2end[i] * dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
tmp = A2end[i] * tmp + dt * k
u = u + B2end[i] * tmp
end

OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
integrator.k[1] = integrator.fsalfirst
integrator.fsalfirst = f(u, p, t + dt) # For interpolation, then FSAL'd
integrator.u = u
function perform_step!(integrator, cache::LowStorageRKTableau{TwoN}, repeat_step = false)
return _perform_step_oop!(integrator, cache)
end

get_fsalfirstlast(cache::LowStorageRK2NCache, u) = (nothing, nothing)
Expand All @@ -42,102 +18,37 @@ function initialize!(integrator, cache::LowStorageRK2NCache)
integrator.kshortsize = 1
resize!(integrator.k, integrator.kshortsize)
integrator.k[1] = k
integrator.f(k, integrator.uprev, integrator.p, integrator.t) # FSAL for interpolation
integrator.f(k, integrator.uprev, integrator.p, integrator.t)
return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
end

@muladd function perform_step!(integrator, cache::LowStorageRK2NCache, repeat_step = false)
(; t, dt, u, f, p) = integrator
(; k, tmp, williamson_condition, stage_limiter!, step_limiter!, thread) = cache
(; A2end, B1, B2end, c2end) = cache.tab

# u1
f(k, u, p, t)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
@.. broadcast = false thread = thread tmp = dt * k
@.. broadcast = false thread = thread u = u + B1 * tmp
# other stages
for i in eachindex(A2end)
if williamson_condition
f(ArrayFuse(tmp, u, (A2end[i], dt, B2end[i])), u, p, t + c2end[i] * dt)
else
@.. broadcast = false thread = thread tmp = A2end[i] * tmp
stage_limiter!(u, integrator, p, t + c2end[i] * dt)
f(k, u, p, t + c2end[i] * dt)
@.. broadcast = false thread = thread tmp = tmp + dt * k
@.. broadcast = false thread = thread u = u + B2end[i] * tmp
end
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
end
stage_limiter!(u, integrator, p, t + dt)
step_limiter!(u, integrator, p, t + dt)
function perform_step!(integrator, cache::LowStorageRK2NCache, repeat_step = false)
return _perform_step_iip!(integrator, cache, cache.tab)
end

# 2C low storage methods
function initialize!(integrator, cache::LowStorageRK2CConstantCache)
integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
function initialize!(integrator, cache::LowStorageRKTableau{TwoC})
integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
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::LowStorageRK2CConstantCache,
repeat_step = false
)
(; t, dt, u, f, p) = integrator
(; A2end, B1, B2end, c2end) = cache

# u1
k = integrator.fsalfirst = f(u, p, t)
integrator.k[1] = integrator.fsalfirst
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
u = u + B1 * dt * k

# other stages
for i in eachindex(A2end)
tmp = u + A2end[i] * dt * k
k = f(tmp, p, t + c2end[i] * dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
u = u + B2end[i] * dt * k
end

integrator.u = u
function perform_step!(integrator, cache::LowStorageRKTableau{TwoC}, repeat_step = false)
return _perform_step_oop!(integrator, cache)
end

function initialize!(integrator, cache::LowStorageRK2CCache)
(; 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
integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
end

@muladd function perform_step!(integrator, cache::LowStorageRK2CCache, repeat_step = false)
(; t, dt, u, f, p) = integrator
(; k, fsalfirst, tmp, stage_limiter!, step_limiter!, thread) = cache
(; A2end, B1, B2end, c2end) = cache.tab

# u1
@.. broadcast = false thread = thread k = integrator.fsalfirst
@.. broadcast = false thread = thread u = u + B1 * dt * k

# other stages
for i in eachindex(A2end)
@.. broadcast = false thread = thread tmp = u + A2end[i] * dt * k
f(k, tmp, p, t + c2end[i] * dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
@.. broadcast = false thread = thread u = u + B2end[i] * dt * k
end
step_limiter!(u, integrator, p, t + dt)
f(k, u, p, t + dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
function perform_step!(integrator, cache::LowStorageRK2CCache, repeat_step = false)
return _perform_step_iip!(integrator, cache, cache.tab)
end

# 3S low storage methods
Expand Down
Loading
Loading