Merge branch 'master' into myb/radau

REQUIRE

@@ -7,7 +7,6 @@ GenericSVD 0.0.2
 NLsolve 0.14.1
 RecursiveArrayTools 0.13.0
 DiffEqDiffTools 0.4.0
-DataStructures 0.4.6
 Reexport
 MuladdMacro 0.2.1
 StaticArrays

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, LDDRK46, 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,7 +186,8 @@ alg_order(alg::OwrenZen3) = 3
 alg_order(alg::OwrenZen4) = 4
 alg_order(alg::OwrenZen5) = 5
 alg_order(alg::LDDRK64) = 4
-alg_order(alg::LDDRK46) = 4
+alg_order(alg::CFRLDDRK64) = 4
+alg_order(alg::NDBLSRK144) = 4
 
 alg_order(alg::DP5) = 5
 alg_order(alg::DP5Threaded) = 5

src/algorithms.jl

@@ -46,7 +46,8 @@ struct OwrenZen3 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct LDDRK64 <: OrdinaryDiffEqAlgorithm end
-struct LDDRK46 <: 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

@@ -536,7 +536,7 @@ end
 
 alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
 
-@cache struct LDDRK46Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
+@cache struct CFRLDDRK64Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
   u::uType
   uprev::uType
   k::rateType
@@ -545,7 +545,7 @@ alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeN
   tab::TabType
 end
 
-struct LDDRK46ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
+struct CFRLDDRK64ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
   α1::T
   α2::T
   α3::T
@@ -563,7 +563,7 @@ struct LDDRK46ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
   c5::T2
   c6::T2
 
-  function LDDRK46ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
+  function CFRLDDRK64ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
     α1 = T(0.17985400977138)
     α2 = T(0.14081893152111)
     α3 = T(0.08255631629428)
@@ -584,14 +584,122 @@ struct LDDRK46ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
   end
 end
 
-function alg_cache(alg::LDDRK46,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+function alg_cache(alg::CFRLDDRK64,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 = LDDRK46ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
-  LDDRK46Cache(u,uprev,k,tmp,fsalfirst,tab)
+  tab = CFRLDDRK64ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
+  CFRLDDRK64Cache(u,uprev,k,tmp,fsalfirst,tab)
 end
 
-function alg_cache(alg::LDDRK46,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
-  LDDRK46ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
+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

@@ -804,7 +804,7 @@ end
   end
 end
 
-function initialize!(integrator,cache::LDDRK46ConstantCache)
+function initialize!(integrator,cache::CFRLDDRK64ConstantCache)
   integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
   integrator.kshortsize = 1
   integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
@@ -814,7 +814,7 @@ function initialize!(integrator,cache::LDDRK46ConstantCache)
   integrator.k[1] = integrator.fsalfirst
 end
 
-@muladd function perform_step!(integrator,cache::LDDRK46ConstantCache,repeat_step=false)
+@muladd function perform_step!(integrator,cache::CFRLDDRK64ConstantCache,repeat_step=false)
   @unpack t,dt,uprev,u,f,p = integrator
   @unpack α1,α2,α3,α4,α5,β1,β2,β3,β4,β5,β6,c2,c3,c4,c5,c6 = cache
 
@@ -844,7 +844,7 @@ end
   integrator.u = u
 end
 
-function initialize!(integrator,cache::LDDRK46Cache)
+function initialize!(integrator,cache::CFRLDDRK64Cache)
   @unpack k,fsalfirst = cache
   integrator.fsalfirst = fsalfirst
   integrator.fsallast = k
@@ -854,7 +854,7 @@ function initialize!(integrator,cache::LDDRK46Cache)
   integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
 end
 
-@muladd function perform_step!(integrator,cache::LDDRK46Cache,repeat_step=false)
+@muladd function perform_step!(integrator,cache::CFRLDDRK64Cache,repeat_step=false)
   @unpack t,dt,uprev,u,f,p = integrator
   @unpack k,fsalfirst,tmp = cache
   @unpack α1,α2,α3,α4,α5,β1,β2,β3,β4,β5,β6,c2,c3,c4,c5,c6 = cache.tab
@@ -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