Merge pull request #1033 from utkarsh530/FRK65

Zero Dissipation ERK6(5)

Author: ChrisRackauckas

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
+         Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65
 
   export SSPRK22, SSPRK33, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK63, SSPRK73, SSPRK83, SSPRK432,
          SSPRKMSVS32, SSPRKMSVS43, SSPRK932, SSPRK54, SSPRK104

src/alg_utils.jl

@@ -25,6 +25,7 @@ isfsal(alg::NDBLSRK144) = false
 isfsal(alg::PDIRK44) = false
 isfsal(alg::DImplicitEuler) = false
 isfsal(alg::RKO65) = false
+isfsal(alg::FRK65) = true
 
 get_current_isfsal(alg, cache) = isfsal(alg)
 get_current_isfsal(alg::CompositeAlgorithm, cache) = isfsal(alg.algs[cache.current])
@@ -164,7 +165,7 @@ alg_order(alg::Anas5) = 5
 alg_order(alg::RK46NL) = 4
 alg_order(alg::KuttaPRK2p5) = 5
 alg_order(alg::RKO65) = 5
-
+alg_order(alg::FRK65) = 6
 
 alg_order(alg::SymplecticEuler) = 1
 alg_order(alg::VelocityVerlet) = 2

src/algorithms.jl

@@ -398,6 +398,11 @@ struct Vern9 <: OrdinaryDiffEqAdaptiveAlgorithm
   lazy::Bool
   Vern9(;lazy=true) = new(lazy)
 end
+struct FRK65{T} <: OrdinaryDiffEqAdaptiveAlgorithm
+  omega::T
+  FRK65(omega=0.0) = new{typeof(omega)}(omega)
+end
+
 
 ################################################################################
 

src/caches/low_order_rk_caches.jl

@@ -617,3 +617,248 @@ function alg_cache(alg::RKO65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUni
   tab = RKO65ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
   RKO65Cache(u,uprev,k, k1,k2,k3,k4,k5,k6, tmp, fsalfirst, tab)
 end
+
+@cache struct FRK65Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
+  u::uType
+  uprev::uType
+  utilde::uType
+  k1::rateType
+  k2::rateType
+  k3::rateType
+  k4::rateType
+  k5::rateType
+  k6::rateType
+  k7::rateType
+  k8::rateType
+  k9::rateType
+  tmp::uType
+  atmp::uNoUnitsType
+  tab::TabType
+end
+
+struct FRK65ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache
+  
+  α21::T1
+  α31::T1
+  α41::T1
+  α51::T1
+  α61::T1
+  α71::T1
+  α81::T1
+  α91::T1
+
+  α32::T1
+
+  α43::T1
+  α53::T1
+  α63::T1
+  α73::T1
+  α83::T1
+
+  α54::T1
+  α64::T1
+  α74::T1
+  α84::T1
+  α94::T1
+  
+  α65::T1
+  α75::T1
+  α85::T1
+  α95::T1
+
+  α76::T1
+  α86::T1
+  α96::T1
+
+  α87::T1
+  α97::T1
+  
+  α98::T1
+  
+  β1::T1
+  # β4::T1
+  # β5::T1
+  # β6::T1
+  β7::T1
+  β8::T1
+
+  β1tilde::T1
+  β4tilde::T1
+  β5tilde::T1
+  β6tilde::T1
+  β7tilde::T1
+  β8tilde::T1
+  β9tilde::T1
+
+  c2::T2
+  c3::T2
+  c4::T2
+  c5::T2
+  c6::T2
+  c7::T2
+  c8::T2
+  c9::T2
+
+  d1::T1
+  d2::T1
+  d3::T1
+  d4::T1
+  d5::T1
+  d6::T1
+  d7::T1
+  d8::T1
+  d9::T1
+  d10::T1
+  d11::T1
+  d12::T1
+  d13::T1
+
+  e1::T1
+  e2::T1
+  e3::T1
+  e4::T1
+  e5::T1
+  e6::T1
+  e7::T1
+  e8::T1
+  e9::T1
+  e10::T1
+  e11::T1
+
+  f1::T1
+  f2::T1
+  f3::T1
+  f4::T1
+  f5::T1
+  f6::T1
+  f7::T1
+  f8::T1
+  f9::T1
+  f10::T1
+  f11::T1
+  
+  function FRK65ConstantCache(T1,T2)
+
+    # elements of Butcher Table
+    α21 = T1(1//89)
+    α31 = T1(-38624//142129)
+    α41 = T1(51//1508)
+    α51 = T1(3259284578//3517556363)
+    α61 = T1(-108363632681//45875676369)
+    α71 = T1(7137368591//11299833148)
+    α81 = T1(8898824396//9828950919)
+    α91 = T1(1026331676//33222204855)
+
+    α32 = T1(51442//142129)
+
+    α43 = T1(153//1508)
+    α53 = T1(-69727055112//19553806387)
+    α63 = T1(80902506271//8700424616)
+    α73 = T1(-33088067061//10572251159)
+    α83 = T1(25673454973//11497947835)
+
+
+    α54 = T1(36230363390//11788838981)
+    α64 = T1(-120088218786//17139312481)
+    α74 = T1(11481363823//3650030081)
+    α84 = T1(-74239028301//15737704666)
+    α94 = T1(1450675392//5936579813)
+   
+    α65 = T1(4533285649//6676940598)
+    α75 = T1(-4096673444//7349814937)
+    α85 = T1(222688842816//44196813415)
+    α95 = T1(4617877550//16762182457)
+    
+    α76 = T1(9911918171//12847192605)
+    α86 = T1(-105204445705//30575217706)
+    α96 = T1(1144867463//6520294355)
+   
+    α87 = T1(8799291910//8966990271)
+    α97 = T1(1822809703//7599996644)
+
+    α98 = T1(79524953//2351253316)
+    
+    β1 = T1(1026331676//33222204855)
+    #β4 = T1(1450675392//5936579813)
+    #β5 = T1(4617877550//16762182457)
+    #β6 = T1(1144867463//6520294355)
+    β7 = T1(1822809703//7599996644)
+    β8 = T1(79524953//2351253316)
+
+    β1tilde = T1(413034411//13925408836)
+    β4tilde = T1(1865954212//7538591735)
+    β5tilde = T1(4451980162//16576017119)
+    β6tilde = T1(1157843020//6320223511)
+    β7tilde = T1(802708729//3404369569)
+    β8tilde = T1(-251398161//17050111121)
+    β9tilde = T1(1//20)
+
+    c2 = T2(1//89)
+    c3 = T2(34//377)
+    c4 = T2(51//377)
+    c5 = T2(14497158//33407747)
+    c6 = T2(9744566553//16002998914)
+    c7 = T2(330//383)
+    c8 = T2(1//1)
+    c9 = T2(1//1)
+
+    d1 = T1(140209127//573775965)
+    d2 = T1(-8530039//263747097)
+    d3 = T1(-308551//104235790)
+    d4 = T1(233511//333733259)
+    d5 = T1(9126//184950985)
+    d6 = T1(22//50434083)
+    d7 = T1(19//424427471)
+    d8 = T1(-28711//583216934059)
+    d9 = T1(-3831531//316297807)
+    d10 = T1(551767//187698280)
+    d11 = T1(9205//210998423)
+    d12 = T1(-250//519462673)
+    d13 = T1(67//327513887)
+  
+    e1=T1(437217689//1587032700)
+    e2=T1(-15824413//592362279)
+    e3=T1(-1563775//341846569)
+    e4=T1(270497//369611210)
+    e5=T1(-26623//453099487)
+    e6=T1(-616297487849)
+    e7=T1(-47682337//491732789)
+    e8=T1(-4778275//287766311)
+    e9=T1(641177//265376522)
+    e10=T1(44633//291742143)
+    e11=T1(611//223639880)
+  
+    f1=T1(44861261//255495624)
+    f2=T1(-11270940//352635157)
+    f3=T1(-182222//232874507)
+    f4=T1(164263//307215200)
+    f5=T1(32184//652060417)
+    f6=T1(-352//171021903)
+    f7=T1(-18395427//101056291)
+    f8=T1(-621686//139501937)
+    f9=T1(2030024//612171255)
+    f10=T1(-711049//7105160932)
+    f11=T1(267//333462710)
+
+    new{T1,T2}(α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11)
+  end
+end
+
+alg_cache(alg::FRK65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = FRK65ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
+
+function alg_cache(alg::FRK65,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
+  tab = FRK65ConstantCache(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)
+  utilde = similar(u)
+  atmp = similar(u,uEltypeNoUnits)
+  tmp = similar(u)
+  FRK65Cache(u, uprev, utilde, k1, k2, k3, k4, k5, k6, k7, k8, k9, tmp, atmp, tab)
+end