Skip to content

Commit cb95437

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #2308 from ParamThakkar123/master
Added ExponentialRK, IMEXMultistep and Linear Solvers
2 parents d27d963 + c6282ce commit cb95437

29 files changed

+404
-235
lines changed

‎.github/workflows/CI.yml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ jobs:
3333
- SSPRK
3434
- LowStorageRK
3535
- QPRK
36+
- Linear
3637
version:
3738
- '1'
3839
steps:

‎lib/OrdinaryDiffEqExponentialRK/Project.toml

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
name = "OrdinaryDiffEqExponentialRK"
2+
uuid = "e0540318-69ee-4070-8777-9e2de6de23de"
3+
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
4+
version = "1.0.0"
5+
6+
[deps]
7+
FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
8+
MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
9+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
10+
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
11+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
12+
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
13+
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
14+
15+
[compat]
16+
julia = "1.10"
17+
18+
[extras]
19+
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
20+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
21+
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
22+
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
23+
24+
[targets]
25+
test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

‎lib/OrdinaryDiffEqExponentialRK/src/OrdinaryDiffEqExponentialRK.jl

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
module OrdinaryDiffEqExponentialRK
2+
3+
import OrdinaryDiffEq: alg_order, alg_adaptive_order, ismultistep, OrdinaryDiffEqExponentialAlgorithm,
4+
_unwrap_val, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
5+
build_jac_config, UJacobianWrapper, @cache, alg_cache, UDerivativeWrapper,
6+
initialize!, perform_step!, @unpack, unwrap_alg, calc_J, calc_J!,
7+
OrdinaryDiffEqAdaptiveExponentialAlgorithm, CompositeAlgorithm,
8+
ExponentialAlgorithm, fsal_typeof, isdtchangeable, calculate_residuals, calculate_residuals!
9+
using RecursiveArrayTools
10+
using MuladdMacro, FastBroadcast
11+
using LinearAlgebra: axpy!, mul!
12+
using DiffEqBase, SciMLBase
13+
using ExponentialUtilities
14+
import RecursiveArrayTools: recursivecopy!
15+
16+
include("algorithms.jl")
17+
include("alg_utils.jl")
18+
include("exponential_rk_caches.jl")
19+
include("exponential_rk_perform_step.jl")
20+
21+
export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A,
22+
EPIRK4s3B,
23+
EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43
24+
end

‎lib/OrdinaryDiffEqExponentialRK/src/alg_utils.jl

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
function isdtchangeable(alg::Union{LawsonEuler, NorsettEuler, ETDRK2, ETDRK3, ETDRK4, HochOst4, ETD2})
2+
false
3+
end # due to caching
4+
5+
alg_order(alg::LawsonEuler) = 1
6+
alg_order(alg::NorsettEuler) = 1
7+
alg_order(alg::ETDRK2) = 2
8+
alg_order(alg::ETDRK3) = 3
9+
alg_order(alg::ETDRK4) = 4
10+
alg_order(alg::HochOst4) = 4
11+
alg_order(alg::Exp4) = 4
12+
alg_order(alg::EPIRK4s3A) = 4
13+
alg_order(alg::EPIRK4s3B) = 4
14+
alg_order(alg::EPIRK5s3) = 5
15+
alg_order(alg::EPIRK5P1) = 5
16+
alg_order(alg::EPIRK5P2) = 5
17+
alg_order(alg::EXPRB53s3) = 5
18+
alg_order(alg::ETD2) = 2
19+
alg_order(alg::Exprb32) = 3
20+
alg_order(alg::Exprb43) = 4
21+
22+
alg_adaptive_order(alg::Exprb32) = 2
23+
alg_adaptive_order(alg::Exprb43) = 4
24+
25+
function DiffEqBase.prepare_alg(
26+
alg::ETD2,
27+
u0::AbstractArray,
28+
p, prob)
29+
alg
30+
end
31+
32+
fsal_typeof(alg::ETD2, rate_prototype) = ETD2Fsal{typeof(rate_prototype)}
33+
function fsal_typeof(alg::CompositeAlgorithm, rate_prototype)
34+
fsal = map(x -> fsal_typeof(x, rate_prototype), alg.algs)
35+
@assert length(unique(fsal))==1 "`fsal_typeof` must be consistent"
36+
return fsal[1]
37+
end
38+
39+
ismultistep(alg::ETD2) = true

‎lib/OrdinaryDiffEqExponentialRK/src/algorithms.jl

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
for Alg in [:LawsonEuler, :NorsettEuler, :ETDRK2, :ETDRK3, :ETDRK4, :HochOst4]
2+
"""
3+
Hochbruck, Marlis, and Alexander Ostermann. “Exponential Integrators.” Acta
4+
Numerica 19 (2010): 209–286. doi:10.1017/S0962492910000048.
5+
"""
6+
@eval struct $Alg{CS, AD, FDT, ST, CJ} <:
7+
OrdinaryDiffEqExponentialAlgorithm{CS, AD, FDT, ST, CJ}
8+
krylov::Bool
9+
m::Int
10+
iop::Int
11+
end
12+
@eval function $Alg(; krylov = false, m = 30, iop = 0, autodiff = true,
13+
standardtag = Val{true}(), concrete_jac = nothing,
14+
chunk_size = Val{0}(),
15+
diff_type = Val{:forward})
16+
$Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff),
17+
diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(krylov,
18+
m,
19+
iop)
20+
end
21+
end
22+
23+
const ETD1 = NorsettEuler # alias
24+
25+
for Alg in [:Exprb32, :Exprb43]
26+
@eval struct $Alg{CS, AD, FDT, ST, CJ} <:
27+
OrdinaryDiffEqAdaptiveExponentialAlgorithm{CS, AD, FDT, ST, CJ}
28+
m::Int
29+
iop::Int
30+
end
31+
@eval function $Alg(; m = 30, iop = 0, autodiff = true, standardtag = Val{true}(),
32+
concrete_jac = nothing, chunk_size = Val{0}(),
33+
diff_type = Val{:forward})
34+
$Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff),
35+
diff_type, _unwrap_val(standardtag),
36+
_unwrap_val(concrete_jac)}(m,
37+
iop)
38+
end
39+
end
40+
for Alg in [:Exp4, :EPIRK4s3A, :EPIRK4s3B, :EPIRK5s3, :EXPRB53s3, :EPIRK5P1, :EPIRK5P2]
41+
@eval struct $Alg{CS, AD, FDT, ST, CJ} <:
42+
OrdinaryDiffEqExponentialAlgorithm{CS, AD, FDT, ST, CJ}
43+
adaptive_krylov::Bool
44+
m::Int
45+
iop::Int
46+
end
47+
@eval function $Alg(; adaptive_krylov = true, m = 30, iop = 0, autodiff = true,
48+
standardtag = Val{true}(), concrete_jac = nothing,
49+
chunk_size = Val{0}(), diff_type = Val{:forward})
50+
$Alg{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type,
51+
_unwrap_val(standardtag), _unwrap_val(concrete_jac)}(adaptive_krylov,
52+
m,
53+
iop)
54+
end
55+
end
56+
57+
"""
58+
ETD2: Exponential Runge-Kutta Method
59+
Second order Exponential Time Differencing method (in development).
60+
"""
61+
struct ETD2 <:
62+
OrdinaryDiffEqExponentialAlgorithm{0, false, Val{:forward}, Val{true}, nothing} end

‎src/caches/linear_nonlinear_caches.jl renamed to ‎lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_caches.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -901,4 +901,4 @@ function alg_cache(alg::ETD2, u, rate_prototype, ::Type{uEltypeNoUnits},
901901
ETD2Cache(
902902
u, uprev, zero(u), zero(rate_prototype), zero(rate_prototype), Phi[1], Phi[2],
903903
Phi[2] + Phi[3], -Phi[3])
904-
end
904+
end

‎src/perform_step/exponential_rk_perform_step.jl renamed to ‎lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_perform_step.jl

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
using LinearAlgebra: axpy!
2-
31
# Helper function to compute the G_nj factors for the classical ExpRK methods
42
@inline _compute_nl(f::SplitFunction, u, p, t, A) = f.f2(u, p, t)
53
@inline _compute_nl(f, u, p, t, A) = f(u, p, t) - A * u
@@ -1635,4 +1633,4 @@ function perform_step!(integrator, cache::ETD2Cache, repeat_step = false)
16351633
integrator.stats.nf += 1
16361634
integrator.stats.nf2 += 1
16371635
@.. broadcast=false integrator.k[2]=fsallast.lin + fsallast.nl
1638-
end
1636+
end

‎test/algconvergence/linear_nonlinear_krylov_tests.jl renamed to ‎lib/OrdinaryDiffEqExponentialRK/test/linear_nonlinear_krylov_tests.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -136,4 +136,4 @@ end
136136
prob = ODEProblem(exp_fun, u0, (0.0, 1.0))
137137
sol = solve(prob, LawsonEuler(krylov = true, m = N); dt = 0.1)
138138
@test sol(1.0) ≈ exp_fun.analytic(u0, nothing, 1.0)
139-
end
139+
end

‎lib/OrdinaryDiffEqExponentialRK/test/runtests.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
using SafeTestsets
2+
3+
@time @safetestset "Linear-Nonlinear Krylov Methods Tests" include("linear_nonlinear_krylov_tests.jl")

‎lib/OrdinaryDiffEqIMEXMultistep/Project.toml

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
name = "OrdinaryDiffEqIMEXMultistep"
2+
uuid = "9f002381-b378-40b7-97a6-27a27c83f129"
3+
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
4+
version = "1.0.0"
5+
6+
[deps]
7+
FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
8+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
9+
10+
[compat]
11+
julia = "1.10"
12+
13+
[extras]
14+
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
15+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
16+
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
17+
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
18+
19+
[targets]
20+
test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

0 commit comments

Comments
 (0)