Skip to content

Commit eedc675

Browse files
author
Your Name
committed
LowStorageRK: unify 2N family on LowStorageRKTableau
1 parent 8cf2552 commit eedc675

5 files changed

Lines changed: 81 additions & 70 deletions

File tree

‎lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,10 @@ using Reexport
2525
include("arrayfuse.jl")
2626
include("algorithms.jl")
2727
include("alg_utils.jl")
28+
include("low_storage_tableaus.jl")
2829
include("low_storage_rk_caches.jl")
2930
include("low_storage_rk_perform_step.jl")
31+
include("generic_low_storage_perform_step.jl")
3032

3133
import PrecompileTools
3234
import Preferences

‎lib/OrdinaryDiffEqLowStorageRK/src/generic_low_storage_perform_step.jl‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
@muladd function _perform_step_iip!(
2+
integrator, cache, tab::LowStorageRKTableau{:two_n}
3+
)
4+
(; t, dt, u, f, p) = integrator
5+
(; k, tmp, williamson_condition, stage_limiter!, step_limiter!, thread) = cache
6+
(; A2end, B1, B2end, c2end) = tab
7+
8+
f(k, u, p, t)
9+
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
10+
@.. broadcast = false thread = thread tmp = dt * k
11+
@.. broadcast = false thread = thread u = u + B1 * tmp
12+
13+
for i in eachindex(A2end)
14+
if williamson_condition
15+
f(ArrayFuse(tmp, u, (A2end[i], dt, B2end[i])), u, p, t + c2end[i] * dt)
16+
else
17+
@.. broadcast = false thread = thread tmp = A2end[i] * tmp
18+
stage_limiter!(u, integrator, p, t + c2end[i] * dt)
19+
f(k, u, p, t + c2end[i] * dt)
20+
@.. broadcast = false thread = thread tmp = tmp + dt * k
21+
@.. broadcast = false thread = thread u = u + B2end[i] * tmp
22+
end
23+
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
24+
end
25+
stage_limiter!(u, integrator, p, t + dt)
26+
step_limiter!(u, integrator, p, t + dt)
27+
return nothing
28+
end
29+
30+
@muladd function _perform_step_oop!(
31+
integrator, tab::LowStorageRKTableau{:two_n}
32+
)
33+
(; t, dt, u, f, p) = integrator
34+
(; A2end, B1, B2end, c2end) = tab
35+
36+
tmp = dt * integrator.fsalfirst
37+
u = u + B1 * tmp
38+
39+
for i in eachindex(A2end)
40+
k = f(u, p, t + c2end[i] * dt)
41+
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
42+
tmp = A2end[i] * tmp + dt * k
43+
u = u + B2end[i] * tmp
44+
end
45+
46+
integrator.k[1] = integrator.fsalfirst
47+
integrator.fsalfirst = f(u, p, t + dt)
48+
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
49+
integrator.u = u
50+
return nothing
51+
end

‎lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl‎

Lines changed: 9 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -17,13 +17,6 @@ get_fsalfirstlast(cache::LowStorageRKMutableCache, u) = (cache.fsalfirst, cache.
1717
thread::Thread
1818
end
1919

20-
struct LowStorageRK2NConstantCache{N, T, T2} <: OrdinaryDiffEqConstantCache
21-
A2end::NTuple{N, T} # A1 is always zero
22-
B1::T
23-
B2end::NTuple{N, T}
24-
c2end::NTuple{N, T2} # c1 is always zero
25-
end
26-
2720
function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
2821
A2 = convert(T, -1.0)
2922
A3 = convert(T, -1.55798)
@@ -44,7 +37,7 @@ function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
4437
c5 = convert(T2, 0.8)
4538
c2end = (c2, c3, c4, c5)
4639

47-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
40+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
4841
end
4942

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

147-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
140+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
148141
end
149142

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

362-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
355+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
363356
end
364357

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

391-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
384+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
392385
end
393386

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

423-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
416+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
424417
end
425418

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

455-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
448+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
456449
end
457450

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

499-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
492+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
500493
end
501494

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

546-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
539+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
547540
end
548541

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

596-
return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
589+
return LowStorageRKTableau{:two_n}(A2end, B1, B2end, c2end)
597590
end
598591

599592
# 2C low storage methods introduced by Calvo, Franco, Rández (2004)

‎lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl‎

Lines changed: 7 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,38 +1,14 @@
1-
# 2N low storage methods
2-
function initialize!(integrator, cache::LowStorageRK2NConstantCache)
3-
integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
1+
function initialize!(integrator, cache::LowStorageRKTableau{:two_n})
2+
integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
43
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
54
integrator.kshortsize = 1
65
integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
7-
8-
# Avoid undefined entries if k is an array of arrays
96
integrator.fsallast = zero(integrator.fsalfirst)
107
return integrator.k[1] = integrator.fsalfirst
118
end
129

13-
@muladd function perform_step!(
14-
integrator, cache::LowStorageRK2NConstantCache,
15-
repeat_step = false
16-
)
17-
(; t, dt, u, f, p) = integrator
18-
(; A2end, B1, B2end, c2end) = cache
19-
20-
# u1
21-
tmp = dt * integrator.fsalfirst
22-
u = u + B1 * tmp
23-
24-
# other stages
25-
for i in eachindex(A2end)
26-
k = f(u, p, t + c2end[i] * dt)
27-
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
28-
tmp = A2end[i] * tmp + dt * k
29-
u = u + B2end[i] * tmp
30-
end
31-
32-
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
33-
integrator.k[1] = integrator.fsalfirst
34-
integrator.fsalfirst = f(u, p, t + dt) # For interpolation, then FSAL'd
35-
integrator.u = u
10+
function perform_step!(integrator, cache::LowStorageRKTableau{:two_n}, repeat_step = false)
11+
return _perform_step_oop!(integrator, cache)
3612
end
3713

3814
get_fsalfirstlast(cache::LowStorageRK2NCache, u) = (nothing, nothing)
@@ -42,35 +18,12 @@ function initialize!(integrator, cache::LowStorageRK2NCache)
4218
integrator.kshortsize = 1
4319
resize!(integrator.k, integrator.kshortsize)
4420
integrator.k[1] = k
45-
integrator.f(k, integrator.uprev, integrator.p, integrator.t) # FSAL for interpolation
21+
integrator.f(k, integrator.uprev, integrator.p, integrator.t)
4622
return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
4723
end
4824

49-
@muladd function perform_step!(integrator, cache::LowStorageRK2NCache, repeat_step = false)
50-
(; t, dt, u, f, p) = integrator
51-
(; k, tmp, williamson_condition, stage_limiter!, step_limiter!, thread) = cache
52-
(; A2end, B1, B2end, c2end) = cache.tab
53-
54-
# u1
55-
f(k, u, p, t)
56-
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
57-
@.. broadcast = false thread = thread tmp = dt * k
58-
@.. broadcast = false thread = thread u = u + B1 * tmp
59-
# other stages
60-
for i in eachindex(A2end)
61-
if williamson_condition
62-
f(ArrayFuse(tmp, u, (A2end[i], dt, B2end[i])), u, p, t + c2end[i] * dt)
63-
else
64-
@.. broadcast = false thread = thread tmp = A2end[i] * tmp
65-
stage_limiter!(u, integrator, p, t + c2end[i] * dt)
66-
f(k, u, p, t + c2end[i] * dt)
67-
@.. broadcast = false thread = thread tmp = tmp + dt * k
68-
@.. broadcast = false thread = thread u = u + B2end[i] * tmp
69-
end
70-
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
71-
end
72-
stage_limiter!(u, integrator, p, t + dt)
73-
step_limiter!(u, integrator, p, t + dt)
25+
function perform_step!(integrator, cache::LowStorageRK2NCache, repeat_step = false)
26+
return _perform_step_iip!(integrator, cache, cache.tab)
7427
end
7528

7629
# 2C low storage methods

‎lib/OrdinaryDiffEqLowStorageRK/src/low_storage_tableaus.jl‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
struct LowStorageRKTableau{form, N, T, T2} <: OrdinaryDiffEqConstantCache
2+
A2end::NTuple{N, T}
3+
B1::T
4+
B2end::NTuple{N, T}
5+
c2end::NTuple{N, T2}
6+
end
7+
8+
function LowStorageRKTableau{form}(
9+
A2end::NTuple{N, T}, B1::T, B2end::NTuple{N, T}, c2end::NTuple{N, T2}
10+
) where {form, N, T, T2}
11+
return LowStorageRKTableau{form, N, T, T2}(A2end, B1, B2end, c2end)
12+
end

0 commit comments

Comments
 (0)