Merge pull request #611 from deeepeshthakur/NDBLSRK124

Implemented NDBLSRK124 from Niegemann, Diehl, Busch (2012)

src/OrdinaryDiffEq.jl

@@ -157,7 +157,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, NDBLSRK144, BS3, BS5, CarpenterKennedy2N54,
+         LDDRK64, CFRLDDRK64, NDBLSRK124, 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

@@ -187,6 +187,7 @@ alg_order(alg::OwrenZen4) = 4
 alg_order(alg::OwrenZen5) = 5
 alg_order(alg::LDDRK64) = 4
 alg_order(alg::CFRLDDRK64) = 4
+alg_order(alg::NDBLSRK124) = 4
 alg_order(alg::NDBLSRK144) = 4
 
 alg_order(alg::DP5) = 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 NDBLSRK124 <: OrdinaryDiffEqAlgorithm end
 struct NDBLSRK144 <: OrdinaryDiffEqAlgorithm end
 struct CarpenterKennedy2N54 <: OrdinaryDiffEqAlgorithm end
 struct ORK256 <: OrdinaryDiffEqAlgorithm end

src/caches/low_order_rk_caches.jl

@@ -596,6 +596,102 @@ function alg_cache(alg::CFRLDDRK64,u,rate_prototype,uEltypeNoUnits,uBottomEltype
   CFRLDDRK64ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
 end
 
+@cache struct NDBLSRK124Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
+  u::uType
+  uprev::uType
+  k::rateType
+  tmp::uType
+  fsalfirst::rateType
+  tab::TabType
+end
+
+struct NDBLSRK124ConstantCache{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
+  β1::T
+  β2::T
+  β3::T
+  β4::T
+  β5::T
+  β6::T
+  β7::T
+  β8::T
+  β9::T
+  β10::T
+  β11::T
+  β12::T
+  c2::T2
+  c3::T2
+  c4::T2
+  c5::T2
+  c6::T2
+  c7::T2
+  c8::T2
+  c9::T2
+  c10::T2
+  c11::T2
+  c12::T2
+
+  function NDBLSRK124ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
+    α2  = T(-0.0923311242368072)
+    α3  = T(-0.9441056581158819)
+    α4  = T(-4.3271273247576394)
+    α5  = T(-2.1557771329026072)
+    α6  = T(-0.9770727190189062)
+    α7  = T(-0.7581835342571139)
+    α8  = T(-1.7977525470825499)
+    α9  = T(-2.6915667972700770)
+    α10 = T(-4.6466798960268143)
+    α11 = T(-0.1539613783825189)
+    α12 = T(-0.5943293901830616)
+    β1  = T(0.0650008435125904)
+    β2  = T(0.0161459902249842)
+    β3  = T(0.5758627178358159)
+    β4  = T(0.1649758848361671)
+    β5  = T(0.3934619494248182)
+    β6  = T(0.0443509641602719)
+    β7  = T(0.2074504268408778)
+    β8  = T(0.6914247433015102)
+    β9  = T(0.3766646883450449)
+    β10 = T(0.0757190350155483)
+    β11 = T(0.2027862031054088)
+    β12 = T(0.2167029365631842)
+    c2  = T2(0.0650008435125904)
+    c3  = T2(0.0796560563081853)
+    c4  = T2(0.1620416710085376)
+    c5  = T2(0.2248877362907778)
+    c6  = T2(0.2952293985641261)
+    c7  = T2(0.3318332506149405)
+    c8  = T2(0.4094724050198658)
+    c9  = T2(0.6356954475753369)
+    c10 = T2(0.6806551557645497)
+    c11 = T2(0.7143773712418350)
+    c12 = T2(0.9032588871651854)
+    new{T,T2}(α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12)
+  end
+end
+
+function alg_cache(alg::NDBLSRK124,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 = NDBLSRK124ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
+  NDBLSRK124Cache(u,uprev,k,tmp,fsalfirst,tab)
+end
+
+function alg_cache(alg::NDBLSRK124,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
+  NDBLSRK124ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
+end
+
 @cache struct NDBLSRK144Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
   u::uType
   uprev::uType

src/perform_step/low_order_rk_perform_step.jl

@@ -888,6 +888,131 @@ end
   f( k,  u, p, t+dt)
 end
 
+function initialize!(integrator,cache::NDBLSRK124ConstantCache)
+  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::NDBLSRK124ConstantCache,repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack α2,α3,α4,α5,α6,α7,α8,α9,α10,α11,α12,β1,β2,β3,β4,β5,β6,β7,β8,β9,β10,β11,β12,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12 = 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
+  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::NDBLSRK124Cache)
+  @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::NDBLSRK124Cache,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,β1,β2,β3,β4,β5,β6,β7,β8,β9,β10,β11,β12,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12 = 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
+
+  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

test/ode/ode_ssprk_tests.jl

@@ -311,7 +311,6 @@ for prob in test_problems_nonlinear
   sim = test_convergence(dts, prob, alg)
   @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()
@@ -327,6 +326,23 @@ for prob in test_problems_nonlinear
   sim = test_convergence(dts, prob, alg)
   @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
 end
+
+# for NDBLSRK124 to be in asymptotic range
+dts = 1 .//2 .^(7:-1:3)
+alg = NDBLSRK124()
+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
+
 #reverting back to original dts
 dts = 1 .//2 .^(8:-1:4)