Fifth order EPIRK methods

src/OrdinaryDiffEq.jl

@@ -164,7 +164,8 @@ module OrdinaryDiffEq
 
   export GenericIIF1, GenericIIF2
 
-  export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A, EPIRK4s3B, ETD2
+  export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A, EPIRK4s3B, 
+         EPIRK5s3, EXPRB53s3, ETD2
 
   export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog,
          McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5,

src/alg_utils.jl

@@ -112,6 +112,8 @@ alg_order(alg::HochOst4) = 4
 alg_order(alg::Exp4) = 4
 alg_order(alg::EPIRK4s3A) = 4
 alg_order(alg::EPIRK4s3B) = 4
+alg_order(alg::EPIRK5s3) = 5
+alg_order(alg::EXPRB53s3) = 5
 alg_order(alg::SplitEuler) = 1
 alg_order(alg::ETD2) = 2
 

src/algorithms.jl

@@ -732,7 +732,7 @@ for Alg in [:LawsonEuler, :NorsettEuler, :ETDRK2, :ETDRK3, :ETDRK4, :HochOst4]
   @eval Base.@pure $Alg(;krylov=false, m=30, iop=0) = $Alg(krylov, m, iop)
 end
 ETD1 = NorsettEuler # alias
-for Alg in [:Exp4, :EPIRK4s3A, :EPIRK4s3B]
+for Alg in [:Exp4, :EPIRK4s3A, :EPIRK4s3B, :EPIRK5s3, :EXPRB53s3]
   @eval struct $Alg <: OrdinaryDiffEqExponentialAlgorithm
     m::Int
     iop::Int

src/caches/linear_nonlinear_caches.jl

@@ -424,6 +424,65 @@ function alg_cache(alg::EPIRK4s3B,u,rate_prototype,uEltypeNoUnits,uBottomEltypeN
   EPIRK4s3BCache(u,uprev,tmp,rtmp,rtmp2,K,A,B,KsCache)
 end
 
+struct EPIRK5s3ConstantCache <: ExpRKConstantCache end
+alg_cache(alg::EPIRK5s3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
+  uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = EPIRK5s3ConstantCache()
+
+struct EPIRK5s3Cache{uType,rateType,matType,KsType} <: ExpRKCache
+  u::uType
+  uprev::uType
+  tmp::uType
+  k::uType
+  rtmp::rateType
+  rtmp2::rateType
+  A::matType
+  B::matType
+  KsCache::KsType
+end
+function alg_cache(alg::EPIRK5s3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
+  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+  tmp, k = (similar(u) for i = 1:2)                   # uType caches
+  rtmp, rtmp2 = (zeros(rate_prototype) for i = 1:2)   # rateType caches
+  # Allocate matrices
+  n = length(u); T = eltype(u)
+  A = Matrix{T}(n, n) # TODO: sparse Jacobian support
+  B = zeros(T, n, 5)
+  # Allocate caches for phiv_timestep
+  maxiter = min(alg.m, n)
+  KsCache = _phiv_timestep_caches(u, maxiter, 4)
+  EPIRK5s3Cache(u,uprev,tmp,k,rtmp,rtmp2,A,B,KsCache)
+end
+
+struct EXPRB53s3ConstantCache <: ExpRKConstantCache end
+alg_cache(alg::EXPRB53s3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
+  uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = EXPRB53s3ConstantCache()
+
+struct EXPRB53s3Cache{uType,rateType,matType,KsType} <: ExpRKCache
+  u::uType
+  uprev::uType
+  tmp::uType
+  rtmp::rateType
+  rtmp2::rateType
+  K::matType
+  A::matType
+  B::matType
+  KsCache::KsType
+end
+function alg_cache(alg::EXPRB53s3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
+  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+  tmp = similar(u)                                    # uType caches
+  rtmp, rtmp2 = (zeros(rate_prototype) for i = 1:2)   # rateType caches
+  # Allocate matrices
+  n = length(u); T = eltype(u)
+  K = Matrix{T}(n, 2)
+  A = Matrix{T}(n, n) # TODO: sparse Jacobian support
+  B = zeros(T, n, 5)
+  # Allocate caches for phiv_timestep
+  maxiter = min(alg.m, n)
+  KsCache = _phiv_timestep_caches(u, maxiter, 4)
+  EXPRB53s3Cache(u,uprev,tmp,rtmp,rtmp2,K,A,B,KsCache)
+end
+
 ####################################
 # Multistep exponential method caches
 

src/perform_step/exponential_rk_perform_step.jl

@@ -799,6 +799,167 @@ function perform_step!(integrator, cache::EPIRK4s3BCache, repeat_step=false)
   # integrator.k is automatically set due to aliasing
 end
 
+function perform_step!(integrator, cache::EPIRK5s3ConstantCache, repeat_step=false)
+  @unpack t,dt,uprev,f,p = integrator
+  A = f.jac(uprev, p, t)
+  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
+  f0 = integrator.fsalfirst # f(uprev) is fsaled
+  kwargs = [(:tol, integrator.opts.reltol), (:iop, alg.iop), (:norm, integrator.opts.internalnorm), (:adaptive, true)]
+
+  # Compute U2 horizontally
+  B = zeros(eltype(f0), length(f0), 4)
+  B[:, 3] = (55 / (8 * dt)) * f0
+  B[:, 4] = (-3025 / (192 * dt^2)) * f0
+  k = phiv_timestep(48dt/55, A, B; kwargs...)
+  U2 = uprev + k
+  R2 = f(U2, p, t + 48dt/55)  - f0 - A*k # remainder of U2
+
+  # Compute U3 horizontally
+  B[:, 2] = (53/5) * f0
+  B[:, 3] = (-648 / (5 * dt)) * f0
+  B[:, 4] = (2916 / (5 * dt^2)) * f0 + (32065 / (1152 * dt^2)) * R2
+  k = phiv_timestep(4dt/9, A, B; kwargs...)
+  U3 = uprev + k
+  R3 = f(U3, p, t + 4dt/9)  - f0 - A*k # remainder of U3
+
+  # Update u (horizontally)
+  B = zeros(eltype(f0), length(f0), 5)
+  B[:, 2] = f0
+  B[:, 4] = (-166375 / (61056 * dt^2)) * R2 + (2187 / (106 * dt^2)) * R3
+  B[:, 5] = (499125 / (27136 * dt^3)) * R2 - (2187 / (106 * dt^3)) * R3
+  u = uprev + phiv_timestep(dt, A, B; kwargs...)
+
+  # Update integrator state
+  integrator.fsallast = f(u, p, t + dt)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.k[2] = integrator.fsallast
+  integrator.u = u
+end
+
+function perform_step!(integrator, cache::EPIRK5s3Cache, repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack tmp,k,rtmp,rtmp2,A,B,KsCache = cache
+  f.jac(A, uprev, p, t)
+  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
+  f0 = integrator.fsalfirst # f(u0) is fsaled
+  kwargs = [(:tol, integrator.opts.reltol), (:iop, alg.iop), (:norm, integrator.opts.internalnorm), 
+    (:adaptive, true), (:caches, KsCache)]
+
+  # Compute U2 horizontally
+  fill!(@view(B[:, 2]), zero(eltype(B)))
+  B[:, 3] .= (55 / (8 * dt)) .* f0
+  B[:, 4] .= (-3025 / (192 * dt^2)) .* f0
+  phiv_timestep!(k, 48dt/55, A, @view(B[:, 1:4]); kwargs...)
+  ## Compute R2
+  @. tmp = uprev + k # tmp is now U2
+  f(rtmp, tmp, p, t + 48dt/55); A_mul_B!(rtmp2, A, k)
+  @. rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
+
+  # Compute U3 horizontally
+  B[:, 2] .= (53/5) .* f0
+  B[:, 3] .= (-648 / (5 * dt)) .* f0
+  B[:, 4] .= (2916 / (5 * dt^2)) .* f0 + (32065 / (1152 * dt^2)) .* rtmp
+  phiv_timestep!(k, 4dt/9, A, @view(B[:, 1:4]); kwargs...)
+  ## Update B matrix using R2
+  B[:, 2] .= f0
+  fill!(@view(B[:, 3]), zero(eltype(B)))
+  B[:, 4] .= (-166375 / (61056 * dt^2)) .* rtmp
+  B[:, 5] .= (499125 / (27136 * dt^3)) .* rtmp
+  ## Compute R3 and update B
+  @. tmp = uprev + k # tmp is now U3
+  f(rtmp, tmp, p, t + 4dt/9); A_mul_B!(rtmp2, A, k)
+  @. rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
+  B[:, 4] .+= (2187 / (106 * dt^2)) .* rtmp
+  B[:, 5] .-= (2187 / (106 * dt^3)) .* rtmp
+
+  # Update u
+  phiv_timestep!(k, dt, A, B; kwargs...)
+  @. u = uprev + k
+
+  # Update integrator state
+  f(integrator.fsallast, u, p, t + dt)
+  # integrator.k is automatically set due to aliasing
+end
+
+function perform_step!(integrator, cache::EXPRB53s3ConstantCache, repeat_step=false)
+  @unpack t,dt,uprev,f,p = integrator
+  A = f.jac(uprev, p, t)
+  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
+  f0 = integrator.fsalfirst # f(uprev) is fsaled
+  kwargs = [(:tol, integrator.opts.reltol), (:iop, alg.iop), (:norm, integrator.opts.internalnorm), (:adaptive, true)]
+
+  # Compute the first group for U2 and U3
+  B = [zeros(f0) f0]
+  K = phiv_timestep([dt/2, 9dt/10], A, B; kwargs...)
+  U2 = uprev + K[:, 1]
+  R2 = f(U2, p, t + dt/2)  - f0 - A*K[:, 1] # remainder of U2
+  U3 = uprev + K[:, 2] # partially
+
+  # Compute the second group for U3
+  B = [zeros(R2) zeros(R2) zeros(R2) R2]
+  K = phiv_timestep([dt/2, 9dt/10], A, B; kwargs...)
+  U3 .+= 216/(25*dt^2) .* K[:, 1] + 8/dt^2 .* K[:, 2]
+  R3 = f(U3, p, t + 9dt/10)  - f0 - A*(U3 - uprev) # remainder of U3
+
+  # Compute the third group for u
+  B = zeros(eltype(f0), length(f0), 5)
+  B[:, 2] = f0
+  B[:, 4] = (18 / dt^2) * R2 - (250 / (81 * dt^2)) * R3
+  B[:, 5] = (-60 / dt^3) * R2 + (500 / (27 * dt^3)) * R3
+  u = uprev + phiv_timestep(dt, A, B; kwargs...)
+
+  # Update integrator state
+  integrator.fsallast = f(u, p, t + dt)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.k[2] = integrator.fsallast
+  integrator.u = u
+end
+
+function perform_step!(integrator, cache::EXPRB53s3Cache, repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack tmp,rtmp,rtmp2,K,A,B,KsCache = cache
+  f.jac(A, uprev, p, t)
+  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg
+  f0 = integrator.fsalfirst # f(u0) is fsaled
+  kwargs = [(:tol, integrator.opts.reltol), (:iop, alg.iop), (:norm, integrator.opts.internalnorm), 
+    (:adaptive, true), (:caches, KsCache)]
+
+  # Compute the first group for U2 and U3
+  B[:, 2] .= f0
+  phiv_timestep!(K, [dt/2, 9dt/10], A, @view(B[:, 1:2]); kwargs...)
+  ## U2 and R2
+  @. tmp = uprev + @view(K[:, 1]) # tmp is now U2
+  f(rtmp, tmp, p, t + dt/2); A_mul_B!(rtmp2, A, @view(K[:, 1]))
+  @. rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
+  @. tmp = uprev + @view(K[:, 2]) # tmp is now U3 (partially)
+
+  # Compute the second group for U3
+  fill!(@view(B[:, 2]), zero(eltype(B)))
+  B[:, 4] .= rtmp
+  phiv_timestep!(K, [dt/2, 9dt/10], A, @view(B[:, 1:4]); kwargs...)
+  ## Update B using R2
+  B[:, 2] .= f0
+  B[:, 4] .= (18 / dt^2) .* rtmp
+  B[:, 5] .= (-60 / dt^3) .* rtmp
+  ## U3 and R3
+  @views tmp .+= 216/(25*dt^2) .* K[:, 1] + 8/dt^2 .* K[:, 2] # tmp is now U3
+  f(rtmp, tmp, p, t + 9dt/10)
+  tmp .-= uprev; A_mul_B!(rtmp2, A, tmp)
+  @. rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
+  ## Update B using R3
+  B[:, 4] .-= (250 / (81 * dt^2)) * rtmp
+  B[:, 5] .+= (500 / (27 * dt^3)) * rtmp
+
+  # Update u
+  du = @view(K[:, 1])
+  phiv_timestep!(du, dt, A, B; kwargs...)
+  @. u = uprev + du
+
+  # Update integrator state
+  f(integrator.fsallast, u, p, t + dt)
+  # integrator.k is automatically set due to aliasing
+end
+
 ######################################################
 # Multistep exponential integrators
 function initialize!(integrator,cache::ETD2ConstantCache)

test/linear_nonlinear_krylov_tests.jl

@@ -50,7 +50,7 @@ end
     prob_ip = ODEProblem{true}(f_ip, u0, (0.0, 1.0))
 
     dt = 0.05; tol=1e-5
-    Algs = [Exp4, EPIRK4s3A, EPIRK4s3B]
+    Algs = [Exp4, EPIRK4s3A, EPIRK4s3B, EXPRB53s3]
     for Alg in Algs
         gc()
         sol = solve(prob, Alg(); dt=dt, internalnorm=Base.norm, reltol=tol)
@@ -62,4 +62,14 @@ end
         @test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)
         println(Alg) # prevent Travis hanging
     end
+
+    gc()
+    sol = solve(prob, EPIRK5s3(); dt=dt, internalnorm=Base.norm, reltol=tol)
+    sol_ref = solve(prob, Tsit5(); reltol=tol)
+    @test_broken isapprox(sol(1.0), sol_ref(1.0); rtol=tol)
+
+    sol = solve(prob_ip, EPIRK5s3(); dt=dt, internalnorm=Base.norm, reltol=tol)
+    sol_ref = solve(prob_ip, Tsit5(); reltol=tol)
+    @test_broken isapprox(sol(1.0), sol_ref(1.0); rtol=tol)
+    println(EPIRK5s3) # prevent Travis hanging
 end