Skip to content

Commit 27da511

Browse files
Harsh SinghHarsh Singh
Harsh Singh
authored and
Harsh Singh
committed
Fix IMEX prototype: remove ARS343 re-export (LTS safe)
1 parent 91d7933 commit 27da511

File tree

4 files changed

+215
-16
lines changed

4 files changed

+215
-16
lines changed

‎lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22,
4646
SFSDIRK5, CFNLIRK3, SFSDIRK6, SFSDIRK7, SFSDIRK8, Kvaerno5, KenCarp4, KenCarp5,
4747
SFSDIRK4, SFSDIRK5, CFNLIRK3, SFSDIRK6,
4848
SFSDIRK7, SFSDIRK8, ESDIRK436L2SA2, ESDIRK437L2SA, ESDIRK547L2SA2, ESDIRK659L2SA,
49-
ARS343
49+
OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm, ARS343
5050

5151
import PrecompileTools
5252
import Preferences

‎lib/OrdinaryDiffEqSDIRK/src/algorithms.jl‎

Lines changed: 29 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1594,8 +1594,36 @@ function ESDIRK659L2SA(;
15941594
)
15951595
end
15961596

1597+
abstract type OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ} <:
1598+
OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} end
1599+
1600+
@doc SDIRK_docstring(
1601+
"3rd order L-stable IMEX ARK method. Uses a generic tableau-driven implementation that supports both split and non-split forms.",
1602+
"ARS343";
1603+
references = "@article{ascher1997implicit,
1604+
title={Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations},
1605+
author={Ascher, Uri M and Ruuth, Steven J and Spiteri, Raymond J},
1606+
journal={Applied Numerical Mathematics},
1607+
volume={25},
1608+
number={2-3},
1609+
pages={151--167},
1610+
year={1997},
1611+
publisher={Elsevier}}",
1612+
extra_keyword_description = """
1613+
- `smooth_est`: TBD
1614+
- `extrapolant`: TBD
1615+
- `controller`: TBD
1616+
- `step_limiter!`: function of the form `limiter!(u, integrator, p, t)`
1617+
""",
1618+
extra_keyword_default = """
1619+
smooth_est = true,
1620+
extrapolant = :linear,
1621+
controller = :PI,
1622+
step_limiter! = trivial_limiter!,
1623+
"""
1624+
)
15971625
struct ARS343{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
1598-
OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
1626+
OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
15991627
linsolve::F
16001628
nlsolve::F2
16011629
precs::P

‎lib/OrdinaryDiffEqSDIRK/src/generic_imex_perform_step.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -25,12 +25,12 @@ function full_cache(c::IMEXCache)
2525
end
2626

2727
function alg_cache(
28-
alg::ARS343, u, rate_prototype, ::Type{uEltypeNoUnits},
28+
alg::OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits},
2929
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
3030
uprev, uprev2, f, t, dt, reltol, p, calck,
3131
::Val{false}, verbose
3232
) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
33-
tab = ARS343Tableau(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
33+
tab = IMEXTableau(alg, constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
3434
γ = tab.Ai[2, 2]
3535
c = tab.c[2]
3636
nlsolver = build_nlsolver(
@@ -41,12 +41,12 @@ function alg_cache(
4141
end
4242

4343
function alg_cache(
44-
alg::ARS343, u, rate_prototype, ::Type{uEltypeNoUnits},
44+
alg::OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits},
4545
::Type{uBottomEltypeNoUnits},
4646
::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
4747
::Val{true}, verbose
4848
) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
49-
tab = ARS343Tableau(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
49+
tab = IMEXTableau(alg, constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
5050
γ = tab.Ai[2, 2]
5151
c = tab.c[2]
5252
nlsolver = build_nlsolver(

‎lib/OrdinaryDiffEqSDIRK/src/imex_tableaus.jl‎

Lines changed: 181 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,113 @@ struct IMEXTableau{T, T2}
1111
s::Int
1212
end
1313

14+
# Dispatch: each algorithm type maps to its tableau constructor
15+
IMEXTableau(::ARS343, T, T2) = ARS343Tableau(T, T2)
16+
17+
#
18+
# KenCarp3 IMEX Tableau
19+
#
20+
21+
function KenCarp3IMEXTableau(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
22+
γ = convert(T, 0.435866521508459)
23+
24+
a31 = convert(T, 0.2576482460664272)
25+
a32 = -convert(T, 0.09351476757488625)
26+
a41 = convert(T, 0.18764102434672383)
27+
a42 = -convert(T, 0.595297473576955)
28+
a43 = convert(T, 0.9717899277217721)
29+
30+
btilde1 = convert(T, 0.027099261876665316)
31+
btilde2 = convert(T, 0.11013520969201586)
32+
btilde3 = -convert(T, 0.10306492520138458)
33+
btilde4 = -convert(T, 0.0341695463672966)
34+
35+
c3 = convert(T2, 0.6)
36+
c2 = 2γ
37+
θ = c3 / c2
38+
α31 = ((1 + (-4θ + 3θ^2)) + (6θ * (1 - θ) / c2) * γ)
39+
α32 = ((-2θ + 3θ^2) + (6θ * (1 - θ) / c2) * γ)
40+
θ = 1 / c2
41+
α41 = ((1 + (-4θ + 3θ^2)) + (6θ * (1 - θ) / c2) * γ)
42+
α42 = ((-2θ + 3θ^2) + (6θ * (1 - θ) / c2) * γ)
43+
44+
ea21 = convert(T, 0.871733043016918)
45+
ea31 = convert(T, 0.5275890119763004)
46+
ea32 = convert(T, 0.0724109880236996)
47+
ea41 = convert(T, 0.3990960076760701)
48+
ea42 = -convert(T, 0.4375576546135194)
49+
ea43 = convert(T, 1.0384616469374492)
50+
eb1 = convert(T, 0.18764102434672383)
51+
eb2 = -convert(T, 0.595297473576955)
52+
eb3 = convert(T, 0.9717899277217721)
53+
eb4 = convert(T, 0.435866521508459)
54+
ebtilde1 = convert(T, 0.027099261876665316)
55+
ebtilde2 = convert(T, 0.11013520969201586)
56+
ebtilde3 = -convert(T, 0.10306492520138458)
57+
ebtilde4 = -convert(T, 0.0341695463672966)
58+
59+
s = 4
60+
Ai = zeros(T, s, s)
61+
Ai[2, 1] = γ
62+
Ai[2, 2] = γ
63+
Ai[3, 1] = a31
64+
Ai[3, 2] = a32
65+
Ai[3, 3] = γ
66+
Ai[4, 1] = a41
67+
Ai[4, 2] = a42
68+
Ai[4, 3] = a43
69+
Ai[4, 4] = γ
70+
71+
bi_vec = zeros(T, s)
72+
bi_vec[1] = a41
73+
bi_vec[2] = a42
74+
bi_vec[3] = a43
75+
bi_vec[4] = γ
76+
77+
Ae = zeros(T, s, s)
78+
Ae[2, 1] = ea21
79+
Ae[3, 1] = ea31
80+
Ae[3, 2] = ea32
81+
Ae[4, 1] = ea41
82+
Ae[4, 2] = ea42
83+
Ae[4, 3] = ea43
84+
85+
be_vec = zeros(T, s)
86+
be_vec[1] = eb1
87+
be_vec[2] = eb2
88+
be_vec[3] = eb3
89+
be_vec[4] = eb4
90+
91+
c_vec = zeros(T2, s)
92+
c_vec[1] = zero(T2)
93+
c_vec[2] = convert(T2, 2γ)
94+
c_vec[3] = c3
95+
c_vec[4] = one(T2)
96+
97+
btilde_vec = zeros(T, s)
98+
btilde_vec[1] = btilde1
99+
btilde_vec[2] = btilde2
100+
btilde_vec[3] = btilde3
101+
btilde_vec[4] = btilde4
102+
103+
ebtilde_vec = zeros(T, s)
104+
ebtilde_vec[1] = ebtilde1
105+
ebtilde_vec[2] = ebtilde2
106+
ebtilde_vec[3] = ebtilde3
107+
ebtilde_vec[4] = ebtilde4
108+
109+
α_mat = zeros(T2, s, s)
110+
α_mat[3, 1] = α31
111+
α_mat[3, 2] = α32
112+
α_mat[4, 1] = α41
113+
α_mat[4, 2] = α42
114+
115+
return IMEXTableau(
116+
Ai, bi_vec, Ae, be_vec, c_vec,
117+
btilde_vec, ebtilde_vec, α_mat, 3, s
118+
)
119+
end
120+
14121
function KenCarp3IMEXTableau(T, T2)
15122
γ = convert(T, 1767732205903 // 4055673282236)
16123

@@ -143,19 +250,23 @@ function KenCarp3IMEXTableau(T, T2)
143250
)
144251
end
145252

146-
function ARS343Tableau(T, T2)
147-
γ = convert(T, 4358665215084590 // 10000000000000000)
253+
#
254+
# ARS343 Tableau
255+
#
256+
257+
function ARS343Tableau(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
258+
γ = convert(T, 0.435866521508459)
148259

149260
s = 4
150261

151-
c2 = γ
152-
c3 = (one(T2) + convert(T2, γ)) / 2
262+
c2 = convert(T2, 0.435866521508459)
263+
c3 = convert(T2, 0.7179332607542295)
153264
c4 = one(T2)
154265

155-
a32_i = (one(T) - γ) / 2
266+
a32_i = convert(T, 0.2820667392457705)
156267

157-
b3_i = (one(T) / 2 - 2γ + γ^2) / ((one(T) - γ) / 2)
158-
b2_i = one(T) - γ - b3_i
268+
b3_i = convert(T, -0.644363170684469)
269+
b2_i = convert(T, 1.20849664917601)
159270

160271
Ai = zeros(T, s, s)
161272
Ai[2, 2] = γ
@@ -167,7 +278,67 @@ function ARS343Tableau(T, T2)
167278

168279
bi_vec = T[zero(T), b2_i, b3_i, γ]
169280

170-
ae21 = γ
281+
ae21 = convert(T, 0.435866521508459)
282+
ae31 = convert(T, 0.321278886)
283+
ae32 = convert(T, 0.3966543748)
284+
ae41 = -convert(T, 0.105858296)
285+
ae42 = convert(T, 0.5529291479)
286+
ae43 = convert(T, 0.5529291479)
287+
288+
Ae = zeros(T, s, s)
289+
Ae[2, 1] = ae21
290+
Ae[3, 1] = ae31
291+
Ae[3, 2] = ae32
292+
Ae[4, 1] = ae41
293+
Ae[4, 2] = ae42
294+
Ae[4, 3] = ae43
295+
296+
be_vec = T[zero(T), b2_i, b3_i, γ]
297+
298+
c_vec = T2[zero(T2), c2, c3, c4]
299+
300+
btilde_vec = T[
301+
zero(T), convert(T, 1.20849664917601), convert(T, -0.644363170684469),
302+
convert(T, -0.564133478491541),
303+
]
304+
ebtilde_vec = T[
305+
zero(T), convert(T, 1.20849664917601), convert(T, -0.644363170684469),
306+
convert(T, -0.564133478491541),
307+
]
308+
309+
α_mat = zeros(T2, s, s)
310+
311+
return IMEXTableau(
312+
Ai, bi_vec, Ae, be_vec, c_vec,
313+
btilde_vec, ebtilde_vec, α_mat, 3, s
314+
)
315+
end
316+
317+
function ARS343Tableau(T, T2)
318+
γ = convert(T, 4358665215084590 // 10000000000000000)
319+
320+
s = 4
321+
322+
c2 = convert(T2, γ)
323+
c3 = (one(T2) + convert(T2, γ)) / 2
324+
c4 = one(T2)
325+
326+
a32_i = (one(T) - γ) / 2
327+
328+
b3_i = (one(T) / 2 - 2γ + γ^2) / ((one(T) - γ) / 2)
329+
b2_i = one(T) - γ - b3_i
330+
331+
Ai = zeros(T, s, s)
332+
Ai[2, 2] = convert(T, γ)
333+
Ai[3, 2] = convert(T, a32_i)
334+
Ai[3, 3] = convert(T, γ)
335+
Ai[4, 2] = convert(T, b2_i)
336+
Ai[4, 3] = convert(T, b3_i)
337+
Ai[4, 4] = convert(T, γ)
338+
339+
bi_vec = T[zero(T), convert(T, b2_i), convert(T, b3_i), convert(T, γ)]
340+
341+
ae21 = convert(T, γ)
171342
ae31 = convert(T, 3212788860 // 10000000000)
172343
ae32 = convert(T, 3966543748 // 10000000000)
173344
ae41 = -convert(T, 1058582960 // 10000000000)
@@ -182,9 +353,9 @@ function ARS343Tableau(T, T2)
182353
Ae[4, 2] = ae42
183354
Ae[4, 3] = ae43
184355

185-
be_vec = T[zero(T), b2_i, b3_i, γ]
356+
be_vec = T[zero(T), convert(T, b2_i), convert(T, b3_i), convert(T, γ)]
186357

187-
c_vec = T2[zero(T2), convert(T2, c2), c3, c4]
358+
c_vec = T2[zero(T2), convert(T2, c2), convert(T2, c3), convert(T2, c4)]
188359

189360
btilde_vec = bi_vec .- T[zero(T), zero(T), zero(T), one(T)]
190361
ebtilde_vec = be_vec .- T[zero(T), zero(T), zero(T), one(T)]

0 commit comments

Comments
 (0)