Merge pull request #1047 from utkarsh530/PFRK87

Phase Fitted RK8(7) with zero dissipation
SciML · Feb 26, 2020 · d565aed · d565aed
2 parents 75ec5dd + 1f870c8
commit d565aed
Show file tree

Hide file tree

Showing 7 changed files with 383 additions and 1 deletion.
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
@@ -176,7 +176,7 @@ module OrdinaryDiffEq
 
   export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5,
          BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
-         Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65
+         Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87
 
   export SSPRK22, SSPRK33, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK63, SSPRK73, SSPRK83, SSPRK432,
          SSPRKMSVS32, SSPRKMSVS43, SSPRK932, SSPRK54, SSPRK104

diff --git a/src/alg_utils.jl b/src/alg_utils.jl
@@ -290,6 +290,7 @@ alg_order(alg::Hairer42) = 4
 alg_order(alg::Feagin10) = 10
 alg_order(alg::Feagin12) = 12
 alg_order(alg::Feagin14) = 14
+alg_order(alg::PFRK87) = 8
 
 alg_order(alg::Rosenbrock23) = 2
 alg_order(alg::Rosenbrock32) = 3

diff --git a/src/algorithms.jl b/src/algorithms.jl
@@ -403,6 +403,11 @@ struct FRK65{T} <: OrdinaryDiffEqAdaptiveAlgorithm
   FRK65(omega=0.0) = new{typeof(omega)}(omega)
 end
 
+struct PFRK87{T} <: OrdinaryDiffEqAdaptiveAlgorithm
+  omega::T
+  PFRK87(omega=0.0) = new{typeof(omega)}(omega)
+end
+
 
 ################################################################################
 

diff --git a/src/caches/high_order_rk_caches.jl b/src/caches/high_order_rk_caches.jl
@@ -123,3 +123,50 @@ function alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNo
 end
 
 alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = TsitPap8ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
+
+@cache struct PFRK87Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
+  u::uType
+  uprev::uType
+  fsalfirst::rateType
+  k2::rateType
+  k3::rateType
+  k4::rateType
+  k5::rateType
+  k6::rateType
+  k7::rateType
+  k8::rateType
+  k9::rateType
+  k10::rateType
+  k11::rateType
+  k12::rateType
+  k13::rateType
+  utilde::uType
+  tmp::uType
+  atmp::uNoUnitsType
+  k::rateType
+  tab::TabType
+end
+
+function alg_cache(alg::PFRK87,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
+  tab = PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
+  k1 = zero(rate_prototype) 
+  k2 = zero(rate_prototype)
+  k3 = zero(rate_prototype) 
+  k4 = zero(rate_prototype)
+  k5 = zero(rate_prototype)
+  k6 = zero(rate_prototype)
+  k7 = zero(rate_prototype)
+  k8 = zero(rate_prototype)
+  k9 = zero(rate_prototype)
+  k10 = zero(rate_prototype)
+  k11 = zero(rate_prototype)
+  k12 = zero(rate_prototype)
+  k13 = zero(rate_prototype)
+  utilde = similar(u)
+  k = zero(rate_prototype)
+  tmp = similar(u)
+  atmp = similar(u,uEltypeNoUnits)
+  PFRK87Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k,tab)
+end
+
+alg_cache(alg::PFRK87,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
diff --git a/src/perform_step/high_order_rk_perform_step.jl b/src/perform_step/high_order_rk_perform_step.jl
@@ -601,3 +601,176 @@ end
   integrator.destats.nf += 1
 end
 =#
+
+function initialize!(integrator, cache::PFRK87ConstantCache)
+  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+  integrator.destats.nf += 1
+  integrator.kshortsize = 2
+  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
+  integrator.k[2] = integrator.fsallast
+end
+
+@muladd function perform_step!(integrator, cache::PFRK87ConstantCache, repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α0908, α1008, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache
+  alg = unwrap_alg(integrator, true)
+  ν = alg.omega*dt
+  νsq = ν^2
+
+  c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(1/100000000)*(νsq^5))
+  d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(1/10000000)*(νsq^5))
+
+  α0807 = c_ν/d_ν
+  α0801 = 0.026876256 + 0.0576576*α0807
+  α0804 = 0.22464336 - 0.944944*α0807
+  α0805 = 0.000369024 - 0.2061696*α0807
+  α0806 = 0.21311136 + 0.093456*α0807
+  α0901 = 0.07239997637512857 + 0.01913119863380767*α0807
+  α0904 = -0.688400520601143 - 0.3135390887207368*α0807
+  α0905 = -0.17301267570583073 - 0.06840852844816077*α0807
+  α0906 = 0.1440060555560846 + 0.031009360422930017*α0807
+  α0907 = 0.9982362892760762 + 0.33180705811215994*α0807
+  α1001 = 0.16261514523236525 - 0.12125171966747463*α0807
+  α1004 = -2.1255544052061124 + 1.9871809612169453*α0807
+  α1005 = -0.216403903283323 + 0.43356675517460624*α0807
+  α1006 = -0.060417230254934076 - 0.1965343807796979*α0807
+  α1007 = 2.4846281621788395 - 2.102961615944379*α0807
+  α1101 = -1.0320124180911034 + 1.061943768952537*α0807
+  α1104 = 13.666683232895137 - 17.40407843561103*α0807
+  α1105 = 0.25990355211486116 - 3.797253476860588*α0807
+  α1106 = -5.759316475814002 + 1.7212824826428488*α0807
+  α1107 = -12.822511612651839 + 18.41810566087623*α0807
+  α1201 = 0.2478349764611783 - 0.06383934946543009*α0807
+  α1204 = -4.593782880309185 + 1.046256005127882*α0807
+  α1205 = -0.39566692537411896 + 0.22827403748244698*α0807
+  α1206 = -3.0673550479691665 - 0.10347586863902129*α0807
+  α1207 = 5.386688702227177 - 1.1072148245058775*α0807
+  α1301 = 0.7332242174431163 - 0.5164807626867616*α0807
+  α1304 = -10.196728938160977 + 8.464545832921925*α0807
+  α1305 = -0.43865244706547707 + 1.846809999910238*α0807
+  α1306 = 0.5693856884667226 - 0.8371528845746959*α0807
+  α1307 = 10.52865228002416 - 8.957722185570706*α0807
+
+  k1 = integrator.fsalfirst
+  k2 = f(uprev+dt*α0201*k1, p, t + c2*dt)
+  k3 = f(uprev+dt*(α0301*k1+α0302*k2), p, t + c3*dt)
+  k4 = f(uprev+dt*(α0401*k1       +α0403*k3), p, t + c4*dt)
+  k5 = f(uprev+dt*(α0501*k1       +α0503*k3+α0504*k4), p, t + c5*dt)
+  k6 = f(uprev+dt*(α0601*k1                +α0604*k4+α0605*k5), p, t + c6*dt)
+  k7 = f(uprev+dt*(α0701*k1                +α0704*k4+α0705*k5+α0706*k6), p, t + c7*dt)
+  k8 = f(uprev+dt*(α0801*k1                +α0804*k4+α0805*k5+α0806*k6+α0807*k7), p, t + c8*dt)
+  k9 = f(uprev+dt*(α0901*k1                +α0904*k4+α0905*k5+α0906*k6+α0907*k7+α0908*k8), p, t + c9*dt)
+  k10 =f(uprev+dt*(α1001*k1                +α1004*k4+α1005*k5+α1006*k6+α1007*k7+α1008*k8+α1009*k9), p, t + c10*dt)
+  k11= f(uprev+dt*(α1101*k1                +α1104*k4+α1105*k5+α1106*k6+α1107*k7+α1108*k8+α1109*k9+α1110*k10), p, t + c11*dt)
+  k12= f(uprev+dt*(α1201*k1                +α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11), p, t+c12*dt)
+  k13= f(uprev+dt*(α1301*k1                +α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10+α1311*k11), p, t+c13*dt)
+  integrator.destats.nf += 12
+  u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k13)
+  if integrator.opts.adaptive
+    utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + β12tilde*k12)
+    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
+    integrator.EEst = integrator.opts.internalnorm(atmp,t)
+  end
+  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
+
+
+function initialize!(integrator, cache::PFRK87Cache)
+  integrator.fsalfirst = cache.fsalfirst
+  integrator.fsallast = cache.k
+  integrator.kshortsize = 2
+  resize!(integrator.k, integrator.kshortsize)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.k[2] = integrator.fsallast
+  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+  integrator.destats.nf += 1
+end
+
+@muladd function perform_step!(integrator, cache::PFRK87Cache, repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack α0201, α0301, α0401, α0501, α0601, α0701, α0302, α0403, α0503, α0504, α0604, α0704, α0605, α0705, α0706, α0908, α1008, α1108, α1208, α1308, α1009, α1109, α1209, α1309, α1110, α1210, α1310, α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde, β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13 = cache.tab
+  @unpack k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k = cache
+
+  alg = unwrap_alg(integrator, true)
+  ν = alg.omega*dt
+  νsq = ν^2
+
+  c_ν = -0.19781108078634084 - (0.164050909125528499*νsq) + (0.042578310088756321*(νsq^2)) - (0.002300513610963998*(νsq^3)) + (0.000033467244551879287*(νsq^4)) - (7.8661142036921924*(10^(-8))*(νsq^5))
+  d_ν = 1-(0.296457092123567400*νsq) + (0.0015793885907465726*(νsq^2)) - (0.00018913011771688527*(νsq^3)) + (0.000017089234650765179*(νsq^4)) - (1.2705211682518626*(10^(-7))*(νsq^5))
+
+  α0807 = c_ν/d_ν
+  α0801 = 0.026876256 + 0.0576576*α0807
+  α0804 = 0.22464336 - 0.944944*α0807
+  α0805 = 0.000369024 - 0.2061696*α0807
+  α0806 = 0.21311136 + 0.093456*α0807
+  α0901 = 0.07239997637512857 + 0.01913119863380767*α0807
+  α0904 = -0.688400520601143 - 0.3135390887207368*α0807
+  α0905 = -0.17301267570583073 - 0.06840852844816077*α0807
+  α0906 = 0.1440060555560846 + 0.031009360422930017*α0807
+  α0907 = 0.9982362892760762 + 0.33180705811215994*α0807
+  α1001 = 0.16261514523236525 - 0.12125171966747463*α0807
+  α1004 = -2.1255544052061124 + 1.9871809612169453*α0807
+  α1005 = -0.216403903283323 + 0.43356675517460624*α0807
+  α1006 = -0.060417230254934076 - 0.1965343807796979*α0807
+  α1007 = 2.4846281621788395 - 2.102961615944379*α0807
+  α1101 = -1.0320124180911034 + 1.061943768952537*α0807
+  α1104 = 13.666683232895137 - 17.40407843561103*α0807
+  α1105 = 0.25990355211486116 - 3.797253476860588*α0807
+  α1106 = -5.759316475814002 + 1.7212824826428488*α0807
+  α1107 = -12.822511612651839 + 18.41810566087623*α0807
+  α1201 = 0.2478349764611783 - 0.06383934946543009*α0807
+  α1204 = -4.593782880309185 + 1.046256005127882*α0807
+  α1205 = -0.39566692537411896 + 0.22827403748244698*α0807
+  α1206 = -3.0673550479691665 - 0.10347586863902129*α0807
+  α1207 = 5.386688702227177 - 1.1072148245058775*α0807
+  α1301 = 0.7332242174431163 - 0.5164807626867616*α0807
+  α1304 = -10.196728938160977 + 8.464545832921925*α0807
+  α1305 = -0.43865244706547707 + 1.846809999910238*α0807
+  α1306 = 0.5693856884667226 - 0.8371528845746959*α0807
+  α1307 = 10.52865228002416 - 8.957722185570706*α0807
+
+  k1 = cache.fsalfirst
+  f(k1, uprev, p, t) 
+  @.. tmp = uprev+dt*α0201*k1
+  f(k2, tmp, p, t + c2*dt)
+  @.. tmp = uprev+dt*(α0301*k1+α0302*k2)
+  f(k3, tmp, p, t + c3*dt)
+  @.. tmp = uprev+dt*(α0401*k1+α0403*k3)
+  f(k4, tmp, p, t + c4*dt)
+  @.. tmp = uprev+dt*(α0501*k1+α0503*k3+α0504*k4)
+  f(k5, tmp, p, t + c5*dt)
+  @.. tmp = uprev+dt*(α0601*k1+α0604*k4+α0605*k5)
+  f(k6, tmp, p, t + c6*dt)
+  @.. tmp = uprev+dt*(α0701*k1+α0704*k4+α0705*k5+α0706*k6)
+  f(k7, tmp, p, t + c7*dt)
+  @.. tmp = uprev+dt*(α0801*k1+α0804*k4+α0805*k5+α0806*k6+α0807*k7)
+  f(k8, tmp, p, t + c8*dt)
+  @.. tmp = uprev+dt*(α0901*k1+α0904*k4+α0905*k5+α0906*k6+α0907*k7+α0908*k8)
+  f(k9, tmp, p, t + c9*dt)
+  @.. tmp = uprev+dt*(α1001*k1+α1004*k4+α1005*k5+α1006*k6+α1007*k7+α1008*k8+α1009*k9)
+  f(k10, tmp, p, t + c10*dt)
+  @.. tmp = uprev+dt*(α1101*k1+α1104*k4+α1105*k5+α1106*k6+α1107*k7+α1108*k8+α1109*k9+α1110*k10)
+  f(k11, tmp, p, t + c11*dt)
+  @.. tmp = uprev+dt*(α1201*k1+α1204*k4+α1205*k5+α1206*k6+α1207*k7+α1208*k8+α1209*k9+α1210*k10+α1211*k11)
+  f(k12, tmp, p, t+c12*dt)
+  @.. tmp = uprev+dt*(α1301*k1+α1304*k4+α1305*k5+α1306*k6+α1307*k7+α1308*k8+α1309*k9+α1310*k10+α1311*k11)
+  f(k13, tmp, p, t+c13*dt)
+  @.. u = uprev + dt*(β1*k1+β6*k6+β7*k7+β8*k8+β9*k9+β10*k10+β11*k11+β12*k12+β13*k13)
+  integrator.destats.nf += 13
+  if integrator.opts.adaptive
+    @.. utilde = dt*(β1tilde*k1 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9 + β10tilde*k10 + β11tilde*k11 + β12tilde*k12)
+    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
+    integrator.EEst = integrator.opts.internalnorm(atmp,t)
+  end
+  f(k, u, p, t+dt)
+  integrator.destats.nf += 1
+  return nothing
+end