Skip to content
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

Jump to bottom

Adding ETD2RK4 exponential integrator algorithm #1661

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
+110 −4
Open

Adding ETD2RK4 exponential integrator algorithm #1661

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
739d990
Initial (v early) edits to start adding ETD2RK4 integrator.
May 6, 2022
be490ed
First draft of ETD2RK4 integrator. Just krylov method to start with.
May 7, 2022
793bf4e
ETD2RK4 algorithm working and validated!
May 7, 2022
9a243e6
Flagged ETD2RK4 as extrapolating in alg_utils.
May 8, 2022
030c449
non-krylov, oop version of ETD2RK4 coded and validated.
May 8, 2022
6515eb7
Added tests and flagged alg_order = 4 in alg_utils.
May 9, 2022
0a271d3
Changed ETD2RK4 so that first step is calulated using DP5.
May 9, 2022
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: 1 addition & 1 deletion src/OrdinaryDiffEq.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -261,7 +261,7 @@ using DocStringExtensions
RosenbrockW6S4OS, ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3

export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A, EPIRK4s3B,
EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43
EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43, ETD2RK4

export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog,
McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5,
Expand Down
4 changes: 4 additions & 0 deletions src/alg_utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ isfsal(alg::SSPRK932) = false
isfsal(alg::SSPRK54) = false
isfsal(alg::SSPRK104) = false

isfsal(alg::ETD2RK4) = false

get_current_isfsal(alg, cache) = isfsal(alg)
get_current_isfsal(alg::CompositeAlgorithm, cache) = isfsal(alg.algs[cache.current])::Bool
all_fsal(alg, cache) = isfsal(alg)
Expand Down Expand Up @@ -256,6 +258,7 @@ alg_extrapolates(alg::ABDF2) = true
alg_extrapolates(alg::SBDF) = true
alg_extrapolates(alg::MEBDF2) = true
alg_extrapolates(alg::MagnusLeapfrog) = true
alg_extrapolates(alg::ETD2RK4) = true

alg_order(alg::Union{OrdinaryDiffEqAlgorithm,DAEAlgorithm}) = error("Order is not defined for this algorithm")
get_current_alg_order(alg::Union{OrdinaryDiffEqAlgorithm,DAEAlgorithm},cache) = alg_order(alg)
Expand Down Expand Up @@ -300,6 +303,7 @@ alg_order(alg::CayleyEuler) = 2
alg_order(alg::ETDRK2) = 2
alg_order(alg::ETDRK3) = 3
alg_order(alg::ETDRK4) = 4
alg_order(alg::ETD2RK4) = 4
alg_order(alg::HochOst4) = 4
alg_order(alg::Exp4) = 4
alg_order(alg::EPIRK4s3A) = 4
Expand Down
3 changes: 2 additions & 1 deletion src/algorithms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3854,11 +3854,12 @@ RosenbrockW6S4OS(;chunk_size=Val{0}(),autodiff=true, standardtag = Val{true}(),

######################################

for Alg in [:LawsonEuler, :NorsettEuler, :ETDRK2, :ETDRK3, :ETDRK4, :HochOst4]
for Alg in [:LawsonEuler, :NorsettEuler, :ETDRK2, :ETDRK3, :ETDRK4, :HochOst4, :ETD2RK4]

"""
Hochbruck, Marlis, and Alexander Ostermann. “Exponential Integrators.” Acta
Numerica 19 (2010): 209–86. doi:10.1017/S0962492910000048.
ETD2RK4 from Krogstad, 2005
"""
@eval struct $Alg{FDT,ST,CJ} <: OrdinaryDiffEqExponentialAlgorithm{FDT,ST,CJ}
krylov::Bool
Expand Down
10 changes: 9 additions & 1 deletion src/caches/linear_nonlinear_caches.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,14 +58,21 @@ function expRK_operators(::HochOst4, dt, A)
B5 = 4P[3] - 8P[4]
return A21, A31, A32, A41, A42, A51, A52, A54, B1, B4, B5
end
function expRK_operators(::ETD2RK4, dt, A)
P = phi(dt * A, 2)
Phalf = phi(dt/2 * A, 2)

return P,Phalf
end

# Unified constructor for constant caches
for (Alg, Cache) in [(:LawsonEuler, :LawsonEulerConstantCache),
(:NorsettEuler, :NorsettEulerConstantCache),
(:ETDRK2, :ETDRK2ConstantCache),
(:ETDRK3, :ETDRK3ConstantCache),
(:ETDRK4, :ETDRK4ConstantCache),
(:HochOst4, :HochOst4ConstantCache)]
(:HochOst4, :HochOst4ConstantCache),
(:ETD2RK4, :ETD2RK4ConstantCache)]
@eval struct $Cache{opType,FType} <: ExpRKConstantCache
ops::opType # precomputed operators
uf::FType # derivative wrapper
Expand Down Expand Up @@ -309,6 +316,7 @@ function alg_cache(alg::HochOst4,u,rate_prototype,::Type{uEltypeNoUnits},::Type{
HochOst4Cache(u,uprev,tmp,dz,rtmp,rtmp2,Au,F2,F3,F4,F5,du1,jac_config,uf,J,ops,KsCache)
end


####################################
# EPIRK method caches
function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int)
Expand Down
88 changes: 88 additions & 0 deletions src/perform_step/exponential_rk_perform_step.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -646,6 +646,94 @@ function perform_step!(integrator, cache::HochOst4Cache, repeat_step=false)
# integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::ETD2RK4ConstantCache, repeat_step=false)

#Draft code for new ETD2RK4 integrator
#Variable names follow those in Krogstad, 2005

@unpack t,dt,uprev,uprev2,f,p = integrator
A = isa(f, SplitFunction) ? f.f1.f : calc_J(integrator, cache) # get linear operator
alg = unwrap_alg(integrator, true)
#F1 = integrator.fsalfirst #not sure if FSAL, so ignore for now
halfdt = dt/2

if integrator.iter == 1
# Initialize the first step using Dorm-Prince-5 (also 4th-order)
aij = [
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We normally write the tableau in https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/tableaus/low_order_rk_tableaus.jl so it is only initialized once.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @YingboMa , sorry for the slow reply. So, I should use one of the tableaus for Dorm-Prince-5 that is already in that file for the first-step initialization? Ok, I will try to make that change and resubmit over the next few days.

0 0 0 0 0 0 0;
1/5 0 0 0 0 0 0;
3/40 9/40 0 0 0 0 0;
44/45 −56/15 32/9 0 0 0 0;
19372/6561 −25360/2187 64448/6561 −212/729 0 0 0;
9017/3168 −355/33 46732/5247 49/176 −5103/18656 0 0;
35/384 0 500/1113 125/192 −2187/6784 11/84 0;
]
k = zeros(length(uprev),7)
u_temp = similar(uprev)

k[:,1] = f(uprev,p,t)

for j = 2:7
u_temp = uprev .+ dt.*(k*aij[j,:])
#kview = @view k[:,j]
k[:,j] = f(u_temp,p,t)
end

u = u_temp

else
if alg.krylov
kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
Nn = _compute_nl(f,uprev,p,t,A)
Nnprev = _compute_nl(f,uprev2,p,t,A)
PhiFbar = phiv_timestep([halfdt,dt],A,[uprev Nn (Nn.-Nnprev)./dt])
Na = _compute_nl(f,PhiFbar[:,1],p,t+halfdt,A)
c = PhiFbar[:,1] .+ (halfdt).*(Na .- (3/2).*Nn .+ (1/2).*Nnprev)
Nc = _compute_nl(f,c,p,t+halfdt,A)
d = PhiFbar[:,2] .+ dt.*expv(halfdt,A,(Nc .- (3/2).*Nn .+ (1/2).*Nnprev))
Nd = _compute_nl(f,d,p,t+dt,A)

if isa(f, SplitFunction)
integrator.destats.nf2 += 3
else
integrator.destats.nf += 3
end

# update u

u = PhiFbar[:,2] .+ (dt/3).*expv(halfdt,A,(Na .+ Nc .- 3 .*Nn .+ Nnprev)) .+ (dt/6).*(Nd .- 2 .*Nn .+ Nnprev)

else

P, Phalf = cache.ops

Nn = _compute_nl(f,uprev,p,t,A)
Nnprev = _compute_nl(f,uprev2,p,t,A)
a = (Phalf[1]*uprev) .+ (halfdt .* Phalf[2] * Nn) .+ (halfdt/2 .* Phalf[3] * (Nn .- Nnprev))
b = (P[1]*uprev) .+ (dt .* P[2] * Nn) .+ (dt .* P[3] * (Nn .- Nnprev))
Na = _compute_nl(f,a,p,t+halfdt,A)
c = a .+ (halfdt).*(Na .- (3/2).*Nn .+ (1/2).*Nnprev)
Nc = _compute_nl(f,c,p,t+halfdt,A)
d = b .+ dt .* Phalf[1] * (Nc .- (3/2).*Nn .+ (1/2).*Nnprev)
Nd = _compute_nl(f,d,p,t+dt,A)

integrator.destats.nf2 += 3

# update u

u = b .+ (dt/3) .* Phalf[1] * (Na .+ Nc .- 3 .*Nn .+ Nnprev) .+ (dt/6).*(Nd .- 2 .*Nn .+ Nnprev)

end
end
# Update integrator state
integrator.fsallast = f(u, p, t + dt)
integrator.destats.nf += 1
integrator.k[1] = integrator.fsalfirst
integrator.k[2] = integrator.fsallast
integrator.u = u
end


#############################################
# EPIRK integrators
function perform_step!(integrator, cache::Exp4ConstantCache, repeat_step=false)
Expand Down
2 changes: 1 addition & 1 deletion test/algconvergence/linear_nonlinear_convergence_tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ using OrdinaryDiffEq: alg_order

Random.seed!(100)
dts = 1 ./2 .^(7:-1:4) #14->7 good plot
for Alg in [LawsonEuler,NorsettEuler,ETDRK2,ETDRK3,ETDRK4,HochOst4,ETD2,KenCarp3,CFNLIRK3]
for Alg in [LawsonEuler,NorsettEuler,ETDRK2,ETDRK3,ETDRK4,HochOst4,ETD2,KenCarp3,CFNLIRK3,ETD2RK4]
sim = test_convergence(dts,prob,Alg())
@test sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.2
end
Expand Down
5 changes: 5 additions & 0 deletions test/algconvergence/linear_nonlinear_krylov_tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,11 @@ end

println(Alg) # prevent Travis hanging
end
for Alg = ETD2RK4 #iip version not implemened yet
sol = solve(prob, Alg(krylov=true, m=20); dt=dt, reltol=tol)
sol_ref = solve(prob, Tsit5(); reltol=tol)
@test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)
end
end

@testset "EPIRK" begin
Expand Down