Add BHR 553 methods

docs/src/theseus.md

@@ -29,6 +29,7 @@ Theseus.DIRK43
 Theseus.CooperSayfy5
 Theseus.CrouzeixRaviart34
 Theseus.HairerWannerSDIRK4
+Theseus.ESDIRK43SA2
 ```
 
 ## Implicit–Explicit (IMEX) Runge–Kutta Methods
@@ -58,6 +59,8 @@ Theseus.ARS233
 Theseus.ARS443
 Theseus.KenCarpARK437
 Theseus.KenCarpARK548
+Theseus.BHR553G1
+Theseus.BHR553G2
 ```
 
 ## Rosenbrock-W Methods

libs/Theseus/src/dirk/tableau.jl

@@ -100,33 +100,91 @@ function RKTableau(alg::DIRK43, RealT)
     return DIRKButcher(a, b, c)
 end
 
+"""
+    Theseus.ESDIRK43SA2()
+
+A fourth-order, six-stage, stiffly accurate ESDIRK method
+with an embedded third-order method for error estimation.
+
+## References
+See Table 9 on p. 242 for the Butcher tableau.
+- Christopher A. Kennedy and Mark H. Carpenter (2019)
+  *Diagonally Implicit Runge–Kutta Methods for stiff ODEs*
+  *Applied Numerical Mathematics* 146:221–244.
+  [DOI: 10.1016/j.apnum.2019.07.008] (https://doi.org/10.1016/j.apnum.2019.07.008)
+"""
+struct ESDIRK43SA2 <: DIRK{6} end
+function RKTableau(alg::ESDIRK43SA2, RealT)
+    nstage = 6
+    a = zeros(RealT, nstage, nstage)
+    a[2, 1] = 31 // 125
+    a[2, 2] = a[2, 1]
+    a[3, 1] = -360286518617 // 7014585480527
+    a[3, 2] = a[3, 1]
+    a[3, 3] = a[2, 2]
+    a[4, 1] = -506388693497 // 5937754990171
+    a[4, 2] = a[4, 1]
+    a[4, 3] = 7149918333491 // 13390931526268
+    a[4, 4] = a[2, 2]
+    a[5, 1] = -7628305438933 // 11061539393788
+    a[5, 2] = a[5, 1]
+    a[5, 3] = 21592626537567 // 14352247503901
+    a[5, 4] = 11630056083252 // 17263101053231
+    a[5, 5] = a[2, 2]
+    a[6, 1] = -12917657251 // 5222094901039
+    a[6, 2] = a[6, 1]
+    a[6, 3] = 5602338284630 // 15643096342197
+    a[6, 4] = 9002339615474 // 18125249312447
+    a[6, 5] = -2420307481369 // 24731958684496
+    a[6, 6] = a[2, 2]
+
+    b = a[6, :]
+    #weights for the embedded method:
+    #b_e = zeros(RealT, nstage)
+    #b_e[1] = -1007911106287 // 12117826057527
+    #b_e[2] = b_e[1]
+    #b_e[3] = 17694008993113 // 35931961998873
+    #b_e[4] = 5816803040497 // 11256217655929
+    #b_e[5] = -538664890905 // 7490061179786
+    #b_e[6] = 2032560730450 // 8872919773257
+
+    c = zeros(RealT, nstage)
+    c[2] = 62 // 125
+    c[3] = 486119545908 // 3346201505189
+    c[4] = 1043 // 1706
+    c[5] = 1361 // 1300
+    c[6] = 1
+    return DIRKButcher(a, b, c)
+end
+
 """
     Theseus.CooperSayfy5()
+
 A fifth-order, five-stage, A-stable DIRK method.
 
 ## References
 - E. Hairer, G. Wanner. Solving ordinary differential equations II: Stiff and Differential-Algebraic Problems.
-  Springer, 1996. 
+  Springer, 1996.
   page.101
-- Cooper, G. J., and A. Sayfy. Semiexplicit Runge-Kutta methods for stiff differential equations. 
-  Mathematics of Computation 33, 
-  no. 146 (1979): 541-556. 
+- Cooper, G. J., and A. Sayfy. Semiexplicit Runge-Kutta methods for stiff differential equations.
+  Mathematics of Computation 33,
+  no. 146 (1979): 541-556.
   doi:10.1090/S0025-5718-1979-0521275-1.
 """
 struct CooperSayfy5 <: DIRK{5} end
 function RKTableau(alg::CooperSayfy5, RealT)
     nstage = 5
     sqrt6 = sqrt(convert(RealT, 6))
-    
+
     γ = (6 - sqrt6) / 10
-    
+
     a = zeros(RealT, nstage, nstage)
-    
+
     a[1, 1] = γ
 
     a[2, 1] = (6 + 5 * sqrt6) / 14
     a[2, 2] = γ
-    
+
     a[3, 1] = (888 + 607 * sqrt6) / 2850
     a[3, 2] = (126 - 161 * sqrt6) / 1425
     a[3, 3] = γ
@@ -141,55 +199,56 @@ function RKTableau(alg::CooperSayfy5, RealT)
     a[5, 3] = (1329 - 544 * sqrt6) / 2500
     a[5, 4] = (-96 + 131 * sqrt6) / 625
     a[5, 5] = γ
-    
+
     b = zeros(RealT, nstage)
     b[1] = 0
     b[2] = 0
     b[3] = 1 // 9
     b[4] = (16 - sqrt6) / 36
     b[5] = (16 + sqrt6) / 36
-    
+
     c = zeros(RealT, nstage)
     c[1] = γ
     c[2] = (6 + 9 * sqrt6) / 35
-    c[3] = 1 
+    c[3] = 1
     c[4] = (4 - sqrt6) / 10
     c[5] = (4 + sqrt6) / 10
-    
+
     return DIRKButcher(a, b, c)
 end
 
 
 """
     Theseus.HairerWannerSDIRK4()
+
 A fourth-order, five-stage, L-stable SDIRK method.
 
 ## References
 - E. Hairer, G. Wanner. Solving ordinary differential equations II: Stiff and Differential-Algebraic Problems.
-  Springer, 1996. 
+  Springer, 1996.
   page. 100
 """
 struct HairerWannerSDIRK4 <: DIRK{5} end
 function RKTableau(alg::HairerWannerSDIRK4, RealT)
     nstage = 5
     a = zeros(RealT, nstage, nstage)
-    
+
     γ = 1 // 4
-    
+
     a[1, 1] = γ
-    
+
     a[2, 1] = 1 // 2
     a[2, 2] = γ
-    
+
     a[3, 1] = 17 // 50
     a[3, 2] = -1 // 25
     a[3, 3] = γ
-    
+
     a[4, 1] = 371 // 1360
     a[4, 2] = -137 // 2720
     a[4, 3] = 15 // 544
     a[4, 4] = γ
-    
+
     a[5, 1] = 25 // 24
     a[5, 2] = -49 // 48
     a[5, 3] = 125 // 16
@@ -202,37 +261,38 @@ function RKTableau(alg::HairerWannerSDIRK4, RealT)
     b[3] = 125 // 16
     b[4] = -85 // 12
     b[5] = 1 // 4
-    
+
     c = zeros(RealT, nstage)
     c[1] = 1 // 4
     c[2] = 3 // 4
     c[3] = 11 // 20
     c[4] = 1 // 2
     c[5] = 1
-    
+
     return DIRKButcher(a, b, c)
 end
 
 
 """
     Theseus.CrouzeixRaviart34()
+
 A fourth-order, three-stage, L-stable SDIRK method.
 
 ## References
 - E. Hairer, G. Wanner. Solving ordinary differential equations II: Stiff and Differential-Algebraic Problems.
-  Springer, 1996. 
+  Springer, 1996.
   page.100
-- M. Crouzeix. Sur l’approximation des équations différentielles opérationnelles linéaires par des méthodes de Runge-Kutta. 
+- M. Crouzeix. Sur l’approximation des équations différentielles opérationnelles linéaires par des méthodes de Runge-Kutta.
   Thèse d'état, Univ. Paris 6 192pp, 1975.
 """
 struct CrouzeixRaviart34 <: DIRK{3} end
 function RKTableau(alg::CrouzeixRaviart34, RealT)
     nstage = 3
-    
+
     sqrt3 = sqrt(convert(RealT, 3))
     γ = (1 / sqrt3) * cospi(one(RealT) / 18) + 1 // 2
     δ = 1 / (6 * (2 * γ - 1)^2)
-    
+
     a = zeros(RealT, nstage, nstage)
 
     a[1, 1] = γ
@@ -243,18 +303,16 @@ function RKTableau(alg::CrouzeixRaviart34, RealT)
     a[3, 1] = 2 * γ
     a[3, 2] = 1 - 4 * γ
     a[3, 3] = γ
-    
+
     b = zeros(RealT, nstage)
     b[1] = δ
     b[2] = 1 - 2 * δ
     b[3] = δ
-    
+
     c = zeros(RealT, nstage)
     c[1] = γ
     c[2] = 1 // 2
     c[3] = 1 - γ
-    
+
     return DIRKButcher(a, b, c)
 end
-
-

libs/Theseus/src/imex/tableau.jl

@@ -531,6 +531,135 @@ function RKTableau(alg::ARS443, RealT)
     return IMEXButcher(a, b, c, a_im, b_im, c_im)
 end
 
+"""
+    Theseus.BHR553G1()
+
+A third order, stiffly accurate, L-stable type II IMEX method developed by
+Boscarino and Russo (2009).
+
+## References
+- Sebastiano Boscarino and Giovanni Russo (2009)
+  On a class of uniformly accurate IMEX Runge-Kutta schemes and 
+  applications to hyperbolic systems with relaxation
+  [DOI: 10.1137/080713562], (https://doi.org/10.1137/080713562)
+"""
+struct BHR553G1 <: RKIMEX{5} end
+function RKTableau(alg::BHR553G1, RealT)
+    # BHR(5,5,3)_g1 IMEX Runge-Kutta - Third order
+    nstage = 5
+    gamma = 424782 // 974569
+    a = zeros(RealT, nstage, nstage)
+    a[2, 1] = 2*gamma
+    a[3, 1] = gamma
+    a[3, 2] = gamma
+    a[4, 1] = -475883375220285986033264 // 594112726933437845704163
+    a[4, 3] = 1866233449822026827708736//594112726933437845704163
+    a[5, 1] = 62828845818073169585635881686091391737610308247 //176112910684412105319781630311686343715753056000
+    a[5, 2] = -302987763081184622639300143137943089 //1535359944203293318639180129368156500
+    a[5, 3] = 262315887293043739337088563996093207 // 297427554730376353252081786906492000
+    a[5, 4] = -987618231894176581438124717087 // 23877337660202969319526901856000
+    b = zeros(RealT, nstage)
+    b[1] = 487698502336740678603511//1181159636928185920260208
+    b[3] = 302987763081184622639300143137943089//1535359944203293318639180129368156500
+    b[4] = -105235928335100616072938218863//2282554452064661756575727198000
+    #
+    b[5] = gamma
+    c = zeros(RealT, nstage)
+    c[2] = 2* gamma
+    c[3] = 902905985686//1035759735069
+    c[4] = 2684624//1147171
+    c[5] = 1
+    a_im = zeros(RealT, nstage, nstage)
+    a_im[2, 1] = gamma
+    a_im[2, 2] = gamma
+    a_im[3, 1] = gamma
+    a_im[4, 1] = -3012378541084922027361996761794919360516301377809610//45123394056585269977907753045030512597955897345819349
+    a_im[5, 1] = b[1]   
+    a_im[3, 2] = -31733082319927313//455705377221960889379854647102
+    a_im[3, 3] = gamma
+    a_im[4, 2] = -62865589297807153294268//102559673441610672305587327019095047
+    a_im[4, 3] = 418769796920855299603146267001414900945214277000//212454360385257708555954598099874818603217167139
+    a_im[4, 4] = gamma
+    a_im[5, 3] = b[3]
+    a_im[5, 4] = b[4]
+    a_im[5, 5] = gamma
+    b_im = zeros(RealT, nstage)
+    b_im[1] = b[1]
+    b_im[3] = b[3]
+    b_im[4] = b[4]
+    b_im[5] = b[5]
+    c_im = zeros(RealT, nstage)
+    c_im[2] = c[2]
+    c_im[3] = c[3]
+    c_im[4] = c[4]
+    c_im[5] = c[5]
+    return IMEXButcher(a, b, c, a_im, b_im, c_im)
+end
+
+"""
+    Theseus.BHR553G2()
+
+A third order, stiffly accurate, L-stable type II IMEX method developed by
+Boscarino and Russo (2009).
+
+## References
+- Sebastiano Boscarino and Giovanni Russo (2009)
+  On a class of uniformly accurate IMEX Runge-Kutta schemes and 
+  applications to hyperbolic systems with relaxation
+  [DOI: 10.1137/080713562], (https://doi.org/10.1137/080713562)
+"""
+struct BHR553G2 <: RKIMEX{5} end
+function RKTableau(alg::BHR553G2, RealT)
+    # BHR(5,5,3)_g2 IMEX Runge-Kutta - Third order
+    nstage = 5
+    gamma2= 2051948 // 3582211
+    a = zeros(RealT, nstage, nstage)
+    a[2, 1] = 2*gamma2
+    a[3, 1] = 473447115440655855452482357894373//1226306256343706154920072735579148
+    a[4, 1] = 37498105210828143724516848//172642583546398006173766007
+    a[5, 1] =  -3409975860212064612303539855622639333030782744869519//5886704102363745137792385361113084313351870216475136
+    a[3, 2] = 129298766034131882323069978722019//1226306256343706154920072735579148
+    a[5, 2] =  -237416352433826978856941795734073//554681702576878342891447163499456
+    a[4, 3] = 76283359742561480140804416//172642583546398006173766007
+    a[5, 3] = 4298159710546228783638212411650783228275//2165398513352098924587211488610407046208
+    a[5, 4] = 6101865615855760853571922289749//272863973025878249803640374568448
+    b = zeros(RealT, nstage)
+    b[1] =  -2032971420760927701493589//38017147656515384190997416
+    b[3] = 2197602776651676983265261109643897073447//945067123279139583549933947379097184164
+    b[4] = -128147215194260398070666826235339//69468482710687503388562952626424
+    b[5] = gamma2
+    c = zeros(RealT, nstage)
+    c[2] = 2* gamma2
+    c[3] = 12015769930846//24446477850549
+    c[4] = 3532944//5360597
+    c[5] = 1
+    a_im = zeros(RealT, nstage, nstage)
+    a_im[2, 1] = gamma2
+    a_im[2, 2] = gamma2
+    a_im[3, 1] = 259252258169672523902708425780469319755//4392887760843243968922388674191715336228
+    a_im[4, 1] = 1103202061574553405285863729195740268785131739395559693754//9879457735937277070641522414590493459028264677925767305837
+    a_im[5, 1] = b[1]   
+    a_im[3, 2] = -172074174703261986564706189586177//1226306256343706154920072735579148
+    a_im[3, 3] = gamma2
+    a_im[4, 2] =  -103754520567058969566542556296087324094//459050363888246734833121482275319954529
+    a_im[4, 3] = 3863207083069979654596872190377240608602701071947128//19258690251287609765240683320611425745736762681950551
+    a_im[4, 4] = gamma2
+    a_im[5, 3] = b[3]
+    a_im[5, 4] = b[4]
+    a_im[5, 5] = gamma2
+    b_im = zeros(RealT, nstage)
+    b_im[1] = b[1]
+    b_im[3] = b[3]
+    b_im[4] = b[4]
+    b_im[5] = b[5]
+    c_im = zeros(RealT, nstage)
+    c_im[2] = c[2]
+    c_im[3] = c[3]
+    c_im[4] = c[4]
+    c_im[5] = c[5]
+    return IMEXButcher(a, b, c, a_im, b_im, c_im)
+end
+
 """
     Theseus.KenCarpARK437()