Skip to content

Commit 4fd5afb

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request SciML#2286 from ParamThakkar123/FIRK
Added FIRK solvers
2 parents acce9cc + 66d3a86 commit 4fd5afb

File tree

19 files changed

+308
-236
lines changed

19 files changed

+308
-236
lines changed

‎lib/OrdinaryDiffEqExtrapolation/test/runtests.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
11
using SafeTestsets
22

3-
@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")
3+
@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")

‎lib/OrdinaryDiffEqFIRK/Project.toml‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
name = "OrdinaryDiffEqFIRK"
2+
uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
3+
authors = ["ParamThakkar123 <[email protected]>"]
4+
version = "1.0.0"
5+
6+
[deps]
7+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
8+
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
9+
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
10+
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
11+
SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
12+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
13+
14+
[compat]
15+
julia = "1.10"
16+
17+
[extras]
18+
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
19+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
20+
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
21+
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
22+
23+
[targets]
24+
test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

‎lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
module OrdinaryDiffEqFIRK
2+
3+
import OrdinaryDiffEq: alg_order, calculate_residuals!,
4+
initialize!, perform_step!, @unpack, unwrap_alg,
5+
calculate_residuals,
6+
OrdinaryDiffEqAlgorithm, OrdinaryDiffEqNewtonAdaptiveAlgorithm,
7+
OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
8+
OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
9+
alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
10+
constvalue, _unwrap_val, du_alias_or_new,
11+
explicit_rk_docstring, trivial_limiter!,
12+
_ode_interpolant!, _ode_addsteps!, AbstractController,
13+
NLStatus, qmax_default, alg_adaptive_order, DEFAULT_PRECS,
14+
UJacobianWrapper, build_J_W, build_jac_config, UDerivativeWrapper,
15+
Convergence, calc_J!, dolinsolve, FastConvergence, calc_J,
16+
stepsize_controller!, islinearfunction, step_accept_controller!, step_reject_controller!,
17+
PredictiveController, alg_can_repeat_jac, NewtonAlgorithm, fac_default_gamma,
18+
get_new_W_γdt_cutoff, get_current_adaptive_order, VerySlowConvergence,
19+
Divergence, isfirk
20+
using MuladdMacro, DiffEqBase, RecursiveArrayTools
21+
using SciMLOperators: AbstractSciMLOperator
22+
using LinearAlgebra: I, UniformScaling, mul!, lu
23+
import LinearSolve
24+
import FastBroadcast: @..
25+
include("algorithms.jl")
26+
include("alg_utils.jl")
27+
include("controllers.jl")
28+
include("firk_caches.jl")
29+
include("firk_tableaus.jl")
30+
include("firk_perform_step.jl")
31+
include("integrator_interface.jl")
32+
33+
export RadauIIA3, RadauIIA5, RadauIIA7
34+
35+
end

‎lib/OrdinaryDiffEqFIRK/src/alg_utils.jl‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA7}) = 8
2+
3+
alg_order(alg::RadauIIA3) = 3
4+
alg_order(alg::RadauIIA5) = 5
5+
alg_order(alg::RadauIIA7) = 7
6+
7+
isfirk(alg::RadauIIA3) = true
8+
isfirk(alg::RadauIIA5) = true
9+
isfirk(alg::RadauIIA7) = true
10+
11+
alg_adaptive_order(alg::RadauIIA3) = 1
12+
alg_adaptive_order(alg::RadauIIA5) = 3
13+
alg_adaptive_order(alg::RadauIIA7) = 5

‎lib/OrdinaryDiffEqFIRK/src/algorithms.jl‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
# FIRK Methods
2+
3+
"""
4+
@article{hairer1999stiff,
5+
title={Stiff differential equations solved by Radau methods},
6+
author={Hairer, Ernst and Wanner, Gerhard},
7+
journal={Journal of Computational and Applied Mathematics},
8+
volume={111},
9+
number={1-2},
10+
pages={93--111},
11+
year={1999},
12+
publisher={Elsevier}
13+
}
14+
15+
RadauIIA3: Fully-Implicit Runge-Kutta Method
16+
An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
17+
"""
18+
struct RadauIIA3{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
19+
OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
20+
linsolve::F
21+
precs::P
22+
extrapolant::Symbol
23+
κ::Tol
24+
maxiters::Int
25+
fast_convergence_cutoff::C1
26+
new_W_γdt_cutoff::C2
27+
controller::Symbol
28+
step_limiter!::StepLimiter
29+
end
30+
31+
function RadauIIA3(; chunk_size = Val{0}(), autodiff = Val{true}(),
32+
standardtag = Val{true}(), concrete_jac = nothing,
33+
diff_type = Val{:forward},
34+
linsolve = nothing, precs = DEFAULT_PRECS,
35+
extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
36+
new_W_γdt_cutoff = 1 // 5,
37+
controller = :Predictive, κ = nothing, maxiters = 10,
38+
step_limiter! = trivial_limiter!)
39+
RadauIIA3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
40+
typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
41+
typeof(κ), typeof(fast_convergence_cutoff),
42+
typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
43+
precs,
44+
extrapolant,
45+
κ,
46+
maxiters,
47+
fast_convergence_cutoff,
48+
new_W_γdt_cutoff,
49+
controller,
50+
step_limiter!)
51+
end
52+
53+
"""
54+
@article{hairer1999stiff,
55+
title={Stiff differential equations solved by Radau methods},
56+
author={Hairer, Ernst and Wanner, Gerhard},
57+
journal={Journal of Computational and Applied Mathematics},
58+
volume={111},
59+
number={1-2},
60+
pages={93--111},
61+
year={1999},
62+
publisher={Elsevier}
63+
}
64+
65+
RadauIIA5: Fully-Implicit Runge-Kutta Method
66+
An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
67+
"""
68+
struct RadauIIA5{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
69+
OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
70+
linsolve::F
71+
precs::P
72+
smooth_est::Bool
73+
extrapolant::Symbol
74+
κ::Tol
75+
maxiters::Int
76+
fast_convergence_cutoff::C1
77+
new_W_γdt_cutoff::C2
78+
controller::Symbol
79+
step_limiter!::StepLimiter
80+
end
81+
82+
function RadauIIA5(; chunk_size = Val{0}(), autodiff = Val{true}(),
83+
standardtag = Val{true}(), concrete_jac = nothing,
84+
diff_type = Val{:forward},
85+
linsolve = nothing, precs = DEFAULT_PRECS,
86+
extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
87+
new_W_γdt_cutoff = 1 // 5,
88+
controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
89+
step_limiter! = trivial_limiter!)
90+
RadauIIA5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
91+
typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
92+
typeof(κ), typeof(fast_convergence_cutoff),
93+
typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
94+
precs,
95+
smooth_est,
96+
extrapolant,
97+
κ,
98+
maxiters,
99+
fast_convergence_cutoff,
100+
new_W_γdt_cutoff,
101+
controller,
102+
step_limiter!)
103+
end
104+
105+
"""
106+
@article{hairer1999stiff,
107+
title={Stiff differential equations solved by Radau methods},
108+
author={Hairer, Ernst and Wanner, Gerhard},
109+
journal={Journal of Computational and Applied Mathematics},
110+
volume={111},
111+
number={1-2},
112+
pages={93--111},
113+
year={1999},
114+
publisher={Elsevier}
115+
}
116+
117+
RadauIIA7: Fully-Implicit Runge-Kutta Method
118+
An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
119+
"""
120+
struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
121+
OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
122+
linsolve::F
123+
precs::P
124+
smooth_est::Bool
125+
extrapolant::Symbol
126+
κ::Tol
127+
maxiters::Int
128+
fast_convergence_cutoff::C1
129+
new_W_γdt_cutoff::C2
130+
controller::Symbol
131+
step_limiter!::StepLimiter
132+
end
133+
134+
function RadauIIA7(; chunk_size = Val{0}(), autodiff = Val{true}(),
135+
standardtag = Val{true}(), concrete_jac = nothing,
136+
diff_type = Val{:forward},
137+
linsolve = nothing, precs = DEFAULT_PRECS,
138+
extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
139+
new_W_γdt_cutoff = 1 // 5,
140+
controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
141+
step_limiter! = trivial_limiter!)
142+
RadauIIA7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
143+
typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
144+
typeof(κ), typeof(fast_convergence_cutoff),
145+
typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
146+
precs,
147+
smooth_est,
148+
extrapolant,
149+
κ,
150+
maxiters,
151+
fast_convergence_cutoff,
152+
new_W_γdt_cutoff,
153+
controller,
154+
step_limiter!)
155+
end

‎lib/OrdinaryDiffEqFIRK/src/controllers.jl‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
2+
@unpack qmin, qmax, gamma = integrator.opts
3+
EEst = DiffEqBase.value(integrator.EEst)
4+
5+
if iszero(EEst)
6+
q = inv(qmax)
7+
else
8+
if fac_default_gamma(alg)
9+
fac = gamma
10+
else
11+
if isfirk(alg)
12+
@unpack iter = integrator.cache
13+
@unpack maxiters = alg
14+
else
15+
@unpack iter, maxiters = integrator.cache.nlsolver
16+
end
17+
fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
18+
end
19+
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
20+
qtmp = DiffEqBase.fastpow(EEst, expo) / fac
21+
@fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
22+
integrator.qold = q
23+
end
24+
q
25+
end
26+
27+
function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
28+
@unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
29+
EEst = DiffEqBase.value(integrator.EEst)
30+
31+
if integrator.success_iter > 0
32+
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
33+
qgus = (integrator.dtacc / integrator.dt) *
34+
DiffEqBase.fastpow((EEst^2) / integrator.erracc, expo)
35+
qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
36+
qacc = max(q, qgus)
37+
else
38+
qacc = q
39+
end
40+
if qsteady_min <= qacc <= qsteady_max
41+
qacc = one(qacc)
42+
end
43+
integrator.dtacc = integrator.dt
44+
integrator.erracc = max(1e-2, EEst)
45+
return integrator.dt / qacc
46+
end
47+
48+
function step_reject_controller!(integrator, controller::PredictiveController, alg)
49+
@unpack dt, success_iter, qold = integrator
50+
integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
51+
end

‎src/caches/firk_caches.jl‎ renamed to ‎lib/OrdinaryDiffEqFIRK/src/firk_caches.jl‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -200,7 +200,6 @@ mutable struct RadauIIA5Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
200200
status::NLStatus
201201
step_limiter!::StepLimiter
202202
end
203-
TruncatedStacktraces.@truncate_stacktrace RadauIIA5Cache 1
204203

205204
function alg_cache(alg::RadauIIA5, u, rate_prototype, ::Type{uEltypeNoUnits},
206205
::Type{uBottomEltypeNoUnits},
@@ -370,7 +369,6 @@ mutable struct RadauIIA7Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
370369
status::NLStatus
371370
step_limiter!::StepLimiter
372371
end
373-
TruncatedStacktraces.@truncate_stacktrace RadauIIA7Cache 1
374372

375373
function alg_cache(alg::RadauIIA7, u, rate_prototype, ::Type{uEltypeNoUnits},
376374
::Type{uBottomEltypeNoUnits},

‎src/perform_step/firk_perform_step.jl‎ renamed to ‎lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl‎

File renamed without changes.

‎src/tableaus/firk_tableaus.jl‎ renamed to ‎lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl‎

File renamed without changes.

‎lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA7},
2+
cache::OrdinaryDiffEqMutableCache)
3+
(cache.tmp, cache.atmp)
4+
end

0 commit comments

Comments
 (0)