Merge pull request #616 from deeepeshthakur/NDBLSRK144

Implemented NDBLSRK144 from Niegemann, Diehl, Busch (2012)

src/OrdinaryDiffEq.jl

@@ -154,7 +154,7 @@ module OrdinaryDiffEq
   export FunctionMap, Euler, Heun, Ralston, Midpoint, SSPRK22,
          SSPRK33, SSPRK53, SSPRK53_2N1,SSPRK53_2N2, SSPRK63, SSPRK73, SSPRK83, SSPRK432, SSPRK932,
          SSPRK54, SSPRK104, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5,
-         LDDRK64, CFRLDDRK64, BS3, BS5, CarpenterKennedy2N54,
+         LDDRK64, CFRLDDRK64, NDBLSRK144, BS3, BS5, CarpenterKennedy2N54,
          ORK256, RK46NL, DP5, DP5Threaded, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
          Vern9,Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5
 

src/alg_utils.jl

@@ -186,6 +186,7 @@ alg_order(alg::OwrenZen4) = 4
 alg_order(alg::OwrenZen5) = 5
 alg_order(alg::LDDRK64) = 4
 alg_order(alg::CFRLDDRK64) = 4
+alg_order(alg::NDBLSRK144) = 4
 
 alg_order(alg::DP5) = 5
 alg_order(alg::DP5Threaded) = 5

src/algorithms.jl

@@ -47,6 +47,7 @@ struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct LDDRK64 <: OrdinaryDiffEqAlgorithm end
 struct CFRLDDRK64 <: OrdinaryDiffEqAlgorithm end
+struct NDBLSRK144 <: OrdinaryDiffEqAlgorithm end
 struct CarpenterKennedy2N54 <: OrdinaryDiffEqAlgorithm end
 struct ORK256 <: OrdinaryDiffEqAlgorithm end
 struct SSPRK22{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm

src/caches/low_order_rk_caches.jl

@@ -595,3 +595,111 @@ end
 function alg_cache(alg::CFRLDDRK64,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
   CFRLDDRK64ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
 end
+
+@cache struct NDBLSRK144Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
+  u::uType
+  uprev::uType
+  k::rateType
+  tmp::uType
+  fsalfirst::rateType
+  tab::TabType
+end
+
+struct NDBLSRK144ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
+  α2::T
+  α3::T
+  α4::T
+  α5::T
+  α6::T
+  α7::T
+  α8::T
+  α9::T
+  α10::T
+  α11::T
+  α12::T
+  α13::T
+  α14::T
+  β1::T
+  β2::T
+  β3::T
+  β4::T
+  β5::T
+  β6::T
+  β7::T
+  β8::T
+  β9::T
+  β10::T
+  β11::T
+  β12::T
+  β13::T
+  β14::T
+  c2::T2
+  c3::T2
+  c4::T2
+  c5::T2
+  c6::T2
+  c7::T2
+  c8::T2
+  c9::T2
+  c10::T2
+  c11::T2
+  c12::T2
+  c13::T2
+  c14::T2
+
+  function NDBLSRK144ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
+    α2  = T(-0.7188012108672410)
+    α3  = T(-0.7785331173421570)
+    α4  = T(-0.0053282796654044)
+    α5  = T(-0.8552979934029281)
+    α6  = T(-3.9564138245774565)
+    α7  = T(-1.5780575380587385)
+    α8  = T(-2.0837094552574054)
+    α9  = T(-0.7483334182761610)
+    α10 = T(-0.7032861106563359)
+    α11 = T(0.0013917096117681)
+    α12 = T(-0.0932075369637460)
+    α13 = T(-0.9514200470875948)
+    α14 = T(-7.1151571693922548)
+    β1  = T(0.0367762454319673)
+    β2  = T(0.3136296607553959)
+    β3  = T(0.1531848691869027)
+    β4  = T(0.0030097086818182)
+    β5  = T(0.3326293790646110)
+    β6  = T(0.2440251405350864)
+    β7  = T(0.3718879239592277)
+    β8  = T(0.6204126221582444)
+    β9  = T(0.1524043173028741)
+    β10 = T(0.0760894927419266)
+    β11 = T(0.0077604214040978)
+    β12 = T(0.0024647284755382)
+    β13 = T(0.0780348340049386)
+    β14 = T(5.5059777270269628)
+    c2  = T2(0.0367762454319673)
+    c3  = T2(0.1249685262725025)
+    c4  = T2(0.2446177702277698)
+    c5  = T2(0.2476149531070420)
+    c6  = T2(0.2969311120382472)
+    c7  = T2(0.3978149645802642)
+    c8  = T2(0.5270854589440328)
+    c9  = T2(0.6981269994175695)
+    c10 = T2(0.8190890835352128)
+    c11 = T2(0.8527059887098624)
+    c12 = T2(0.8604711817462826)
+    c13 = T2(0.8627060376969976)
+    c14 = T2(0.8734213127600976)
+    new{T,T2}(α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, α14, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, β14, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14)
+  end
+end
+
+function alg_cache(alg::NDBLSRK144,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+  tmp = similar(u)
+  k = zero(rate_prototype)
+  fsalfirst = zero(rate_prototype)
+  tab = NDBLSRK144ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
+  NDBLSRK144Cache(u,uprev,k,tmp,fsalfirst,tab)
+end
+
+function alg_cache(alg::NDBLSRK144,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
+  NDBLSRK144ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
+end

src/perform_step/low_order_rk_perform_step.jl

@@ -887,3 +887,142 @@ end
 
   f( k,  u, p, t+dt)
 end
+
+function initialize!(integrator,cache::NDBLSRK144ConstantCache)
+  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+  integrator.kshortsize = 1
+  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+
+  # Avoid undefined entries if k is an array of arrays
+  integrator.fsallast = zero(integrator.fsalfirst)
+  integrator.k[1] = integrator.fsalfirst
+end
+
+@muladd function perform_step!(integrator,cache::NDBLSRK144ConstantCache,repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, α14, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, β14, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14 = cache
+
+  # u0
+  u   = uprev
+  #u1
+  tmp = dt*integrator.fsalfirst
+  u   = u + β1*tmp
+  #u2
+  tmp = α2*tmp + dt*f(u, p, t+c2*dt)
+  u   = u + β2*tmp
+  #u3
+  tmp = α3*tmp + dt*f(u, p, t+c3*dt)
+  u   = u + β3*tmp
+  #u4
+  tmp = α4*tmp + dt*f(u, p, t+c4*dt)
+  u   = u + β4*tmp
+  #u5
+  tmp = α5*tmp + dt*f(u, p, t+c5*dt)
+  u   = u + β5*tmp
+  #u6
+  tmp = α6*tmp + dt*f(u, p, t+c6*dt)
+  u   = u + β6*tmp
+  #u7
+  tmp = α7*tmp + dt*f(u, p, t+c7*dt)
+  u   = u + β7*tmp
+  #u8
+  tmp = α8*tmp + dt*f(u, p, t+c8*dt)
+  u   = u + β8*tmp
+  #u9
+  tmp = α9*tmp + dt*f(u, p, t+c9*dt)
+  u   = u + β9*tmp
+  #u10
+  tmp = α10*tmp + dt*f(u, p, t+c10*dt)
+  u   = u + β10*tmp
+  #u11
+  tmp = α11*tmp + dt*f(u, p, t+c11*dt)
+  u   = u + β11*tmp
+  #u12
+  tmp = α12*tmp + dt*f(u, p, t+c12*dt)
+  u   = u + β12*tmp
+  #u13
+  tmp = α13*tmp + dt*f(u, p, t+c13*dt)
+  u   = u + β13*tmp
+  #u14
+  tmp = α14*tmp + dt*f(u, p, t+c14*dt)
+  u   = u + β14*tmp
+  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
+  integrator.k[1] = integrator.fsalfirst
+  integrator.u = u
+end
+
+function initialize!(integrator,cache::NDBLSRK144Cache)
+  @unpack k,fsalfirst = cache
+  integrator.fsalfirst = fsalfirst
+  integrator.fsallast = k
+  integrator.kshortsize = 1
+  resize!(integrator.k, integrator.kshortsize)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
+end
+
+@muladd function perform_step!(integrator,cache::NDBLSRK144Cache,repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack k,fsalfirst,tmp = cache
+  @unpack α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, α14, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, β14, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14 = cache.tab
+
+  #u0
+  @. u   = uprev
+  #u1
+  @. tmp = dt*fsalfirst
+  @. u   = u + β1*tmp
+  #u2
+  f(k, u, p, t+c2*dt)
+  @. tmp = α2*tmp + dt*k
+  @. u   = u + β2*tmp
+  #u3
+  f(k, u, p, t+c3*dt)
+  @. tmp = α3*tmp + dt*k
+  @. u   = u + β3*tmp
+  #u4
+  f(k, u, p, t+c4*dt)
+  @. tmp = α4*tmp + dt*k
+  @. u   = u + β4*tmp
+  #u5
+  f(k, u, p, t+c5*dt)
+  @. tmp = α5*tmp + dt*k
+  @. u   = u + β5*tmp
+  #u6
+  f(k, u, p, t+c6*dt)
+  @. tmp = α6*tmp + dt*k
+  @. u   = u + β6*tmp
+  #u7
+  f(k, u, p, t+c7*dt)
+  @. tmp = α7*tmp + dt*k
+  @. u   = u + β7*tmp
+  #u8
+  f(k, u, p, t+c8*dt)
+  @. tmp = α8*tmp + dt*k
+  @. u   = u + β8*tmp
+  #u9
+  f(k, u, p, t+c9*dt)
+  @. tmp = α9*tmp + dt*k
+  @. u   = u + β9*tmp
+  #u10
+  f(k, u, p, t+c10*dt)
+  @. tmp = α10*tmp + dt*k
+  @. u   = u + β10*tmp
+  #u11
+  f(k, u, p, t+c11*dt)
+  @. tmp = α11*tmp + dt*k
+  @. u   = u + β11*tmp
+  #u12
+  f(k, u, p, t+c12*dt)
+  @. tmp = α12*tmp + dt*k
+  @. u   = u + β12*tmp
+  #u13
+  f(k, u, p, t+c13*dt)
+  @. tmp = α13*tmp + dt*k
+  @. u   = u + β13*tmp
+  #u14
+  f(k, u, p, t+c14*dt)
+  @. tmp = α14*tmp + dt*k
+  @. u   = u + β14*tmp
+
+  f( k,  u, p, t+dt)
+end

test/ode/ode_ssprk_tests.jl

@@ -312,6 +312,8 @@ for prob in test_problems_nonlinear
   @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
 end
 
+# for CFRLDDRK64 to be in asymptotic range
+dts = 1 .//2 .^(7:-1:4)
 alg = CFRLDDRK64()
 for prob in test_problems_only_time
 	sim = test_convergence(dts, prob, alg)
@@ -325,3 +327,19 @@ for prob in test_problems_nonlinear
 	sim = test_convergence(dts, prob, alg)
 	@test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
 end
+#reverting back to original dts
+dts = 1 .//2 .^(8:-1:4)
+
+alg = NDBLSRK144()
+for prob in test_problems_only_time
+  sim = test_convergence(dts, prob, alg)
+  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
+end
+for prob in test_problems_linear
+  sim = test_convergence(dts, prob, alg)
+  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
+end
+for prob in test_problems_nonlinear
+  sim = test_convergence(dts, prob, alg)
+  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
+end