Add Sofroniou-Spaletta-5-3-4 method only

CHANGELOG.md

@@ -27,6 +27,8 @@ Fixed scaling bug in `SUNMatScaleAddI_Sparse` for non-square matrices.
 Fixed missing soversions in some `SUNLinearSolver` and `SUNNonlinearSolver`
 CMake targets.
 
+Added the fourth order ERK method `ARKODE_SOFRONIOU_SPALETTA_5_3_4`.
+
 ## Changes to SUNDIALS in release 6.6.1
 
 Updated the Tpetra NVector interface to support Trilinos 14.

doc/arkode/guide/source/Butcher.rst

@@ -331,6 +331,43 @@ This is the default 3th order slow and fast MRIStep method (from
    Linear stability region for the Knoth-Wolke method
 
 
+.. _Butcher.Sofroniou_Spaletta:
+
+Sofroniou-Spaletta-5-3-4
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. index:: Sofroniou-Spaletta-5-3-4 ERK method
+
+Accessible via the constant ``ARKODE_SOFRONIOU_SPALETTA_5_3_4`` to
+:c:func:`ARKStepSetTableNum`, :c:func:`ERKStepSetTableNum`
+or :c:func:`ARKodeButcherTable_LoadERK`.
+Accessible via the string ``"ARKODE_SOFRONIOU_SPALETTA_5_3_4"`` to
+:c:func:`ARKStepSetTableName`, :c:func:`ERKStepSetTableName` or
+:c:func:`ARKodeButcherTable_LoadERKByName`.
+(from :cite:p:`Sof:04`).
+
+.. math::
+
+   \renewcommand{\arraystretch}{1.5}
+   \begin{array}{r|ccccc}
+     0 & 0 & 0 & 0 & 0 & 0 \\
+     \frac{2}{5} & \frac{2}{5} & 0 & 0 & 0 & 0 \\
+     \frac{3}{5} & -\frac{3}{20} & \frac{3}{4} & 0 & 0 & 0 \\
+     1 & \frac{19}{44} & -\frac{15}{44} & \frac{10}{11} & 0 & 0 \\
+     1 & \frac{11}{72} & \frac{25}{72} & \frac{25}{72} & \frac{11}{72} & 0 \\
+     4 & \frac{11}{72} & \frac{25}{72} & \frac{25}{72} & \frac{11}{72} & 0 \\
+     3 & \frac{1251515}{8970912} & \frac{3710105}{8970912} & \frac{2519695}{8970912} & \frac{61105}{8970912} & \frac{119041}{747576} \\
+   \end{array}
+
+.. figure:: /figs/arkode/sofroniou_spaletta_erk_stab_region.png
+   :scale: 50 %
+   :align: center
+
+   Linear stability region for the Sofroniou-Spaletta method.  The method's
+   region is outlined in blue; the embedding's region is in red.
+
+
+
 
 .. _Butcher.Zonneveld:
 

doc/arkode/guide/source/Constants.rst

@@ -80,6 +80,8 @@ contains the ARKODE output constants.
    +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ARKODE_ARK324L2SA_ERK_4_2_3`          | Use the ARK-4-2-3 ERK method.                              |
    +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SOFRONIOU_SPALETTA_5_3_4`      | Use the Sofroniou-Spaletta-5-3-4 ERK method.               |
+   +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ARKODE_ZONNEVELD_5_3_4`               | Use the Zonneveld-5-3-4 ERK method.                        |
    +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ARKODE_ARK436L2SA_ERK_6_3_4`          | Use the ARK-6-3-4 ERK method.                              |

doc/arkode/guide/source/Introduction.rst

@@ -153,6 +153,8 @@ Fixed scaling bug in ``SUNMatScaleAddI_Sparse`` for non-square matrices.
 Fixed missing soversions in some ``SUNLinearSolver`` and ``SUNNonlinearSolver``
 CMake targets.
 
+Added the fourth order ERK method ``ARKODE_SOFRONIOU_SPALETTA_5_3_4``.
+
 Changes in v5.6.1
 -----------------
 

doc/shared/sundials.bib

@@ -2053,6 +2053,18 @@ @article{Sof:02
   pages   = {506--515}
 }
 
+@article{Sof:04,
+  author = {M. Sofroniou and G. Spaletta},
+  title = {Construction of explicit {Runge-Kutta} pairs with stiffness detection},
+  journal = {Mathematical and Computer Modelling},
+  volume = {40},
+  number = {11},
+  pages = {1157-1169},
+  year = {2004},
+  issn = {0895-7177},
+  doi = {10.1016/j.mcm.2005.01.010}
+}
+
 @article{Ver:78,
   author  = {Verner, J.H},
   title   = {Explicit Runge-Kutta methods with estimates of the local truncation error},

include/arkode/arkode_butcher_erk.h

@@ -85,7 +85,8 @@ typedef enum {
   ARKODE_ARK437L2SA_ERK_7_3_4,
   ARKODE_ARK548L2SAb_ERK_8_4_5,
   ARKODE_ARK2_ERK_3_1_2,
-  ARKODE_MAX_ERK_NUM = ARKODE_ARK2_ERK_3_1_2
+  ARKODE_SOFRONIOU_SPALETTA_5_3_4,
+  ARKODE_MAX_ERK_NUM = ARKODE_SOFRONIOU_SPALETTA_5_3_4
 } ARKODE_ERKTableID;
 
 /* Accessor routine to load built-in ERK table */

src/arkode/arkode_butcher_erk.def

@@ -41,6 +41,7 @@
      ARKODE_HEUN_EULER_2_1_2             Y
      ARKODE_BOGACKI_SHAMPINE_4_2_3       Y
      ARKODE_ARK324L2SA_ERK_4_2_3*        N
+     ARKODE_SOFRONIOU_SPALETTA_5_3_4     Y
      ARKODE_ZONNEVELD_5_3_4              Y
      ARKODE_ARK436L2SA_ERK_6_3_4*        N
      ARKODE_ARK437L2SA_ERK_7_3_4*        N
@@ -165,6 +166,39 @@ ARK_BUTCHER_TABLE(ARKODE_ARK324L2SA_ERK_4_2_3, { /* ARK3(2)4L[2]SA-ERK */
     return B;
   })
 
+ARK_BUTCHER_TABLE(ARKODE_SOFRONIOU_SPALETTA_5_3_4, { /* Sofroniou-Spaletta-ERK */
+    ARKodeButcherTable B = ARKodeButcherTable_Alloc(5, SUNTRUE);
+    B->q = 4;
+    B->p = 3;
+    B->A[1][0] = RCONST(2.0) / RCONST(5.0);
+    B->A[2][0] = RCONST(-3.0) / RCONST(20.0);
+    B->A[2][1] = RCONST(3.0) / RCONST(4.0);
+    B->A[3][0] = RCONST(19.0)/RCONST(44.0);
+    B->A[3][1] = RCONST(-15.0)/RCONST(44.0);
+    B->A[3][2] = RCONST(10.0)/RCONST(11.0);
+    B->A[4][0] = RCONST(11.0)/RCONST(72.0);
+    B->A[4][1] = RCONST(25.0)/RCONST(72.0);
+    B->A[4][2] = RCONST(25.0)/RCONST(72.0);
+    B->A[4][3] = RCONST(11.0)/RCONST(72.0);
+
+    B->b[0] = RCONST(11.0)/RCONST(72.0);
+    B->b[1] = RCONST(25.0)/RCONST(72.0);
+    B->b[2] = RCONST(25.0)/RCONST(72.0);
+    B->b[3] = RCONST(11.0)/RCONST(72.0);
+
+    B->d[0] = RCONST(1251515.0)/RCONST(8970912.0);
+    B->d[1] = RCONST(3710105.0)/RCONST(8970912.0);
+    B->d[2] = RCONST(2519695.0)/RCONST(8970912.0);
+    B->d[3] = RCONST(61105.0)/RCONST(8970912.0);
+    B->d[4] = RCONST(119041.0)/RCONST(747576.0);
+
+    B->c[1] = RCONST(2.0) / RCONST(5.0);
+    B->c[2] = RCONST(3.0) / RCONST(5.0);
+    B->c[3] = RCONST(1.0);
+    B->c[4] = RCONST(1.0);
+    return B;
+  })
+
 ARK_BUTCHER_TABLE(ARKODE_ZONNEVELD_5_3_4, { /* Zonneveld */
     ARKodeButcherTable B = ARKodeButcherTable_Alloc(5, SUNTRUE);
     B->q = 4;