Feature/new tables (#359)

Additional embedded explicit Runge--Kutta methods to ARKODE. --------- Co-authored-by: Steven Roberts <[email protected]> Co-authored-by: David J. Gardner <[email protected]>
LLNL · Oct 28, 2023 · 6c60a56 · 6c60a56
1 parent 9717d6d
commit 6c60a56
Show file tree

Hide file tree

Showing 20 changed files with 901 additions and 128 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -32,7 +32,11 @@ Added Fortran support for the LAPACK dense `SUNLinearSolver` implementations.
 Fixed the build system support for MAGMA when using a NVIDIA HPC SDK installation of CUDA
 and fixed the targets used for rocBLAS and rocSPARSE.
 
-Added the fourth order ERK method `ARKODE_SOFRONIOU_SPALETTA_5_3_4`.
+Added the third order ERK method `ARKODE_SHU_OSHER_3_2_3`, the fourth order
+ERK method `ARKODE_SOFRONIOU_SPALETTA_5_3_4`, the sixth order ERK method
+`ARKODE_VERNER_9_5_6`, the seventh order ERK method `ARKODE_VERNER_10_6_7`,
+the eighth order ERK method `ARKODE_VERNER_13_7_8`, and the ninth order ERK
+method `ARKODE_VERNER_16_8_9`.
 
 
 ## Changes to SUNDIALS in release 6.6.1

diff --git a/doc/arkode/guide/source/Butcher.rst b/doc/arkode/guide/source/Butcher.rst
diff --git a/doc/arkode/guide/source/Constants.rst b/doc/arkode/guide/source/Constants.rst
@@ -106,6 +106,16 @@ contains the ARKODE output constants.
    +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ARKODE_ARK548L2SAb_ERK_8_4_5`         | Use the ARK-8-4-5b ERK method.                             |
    +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SHU_OSHER_3_2_3`               | Use the Shu-Osher-3-2-3 ERK method.                        |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_VERNER_9_5_6`                  | Use the Verner-9-5-6 ERK method.                           |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_VERNER_10_6_7`                 | Use the Verner-10-6-7 ERK method.                          |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_VERNER_13_7_8`                 | Use the Verner-13-7-8 ERK method.                          |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_VERNER_16_8_9`                 | Use the Verner-16-8-9 ERK method.                          |
+   +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ARKSTEP_DEFAULT_ERK_2`                | Use ARKStep's default second-order ERK method              |
    |                                               | (ARKODE_HEUN_EULER_2_1_2).                                 |
    +-----------------------------------------------+------------------------------------------------------------+
@@ -121,9 +131,15 @@ contains the ARKODE output constants.
    | :index:`ARKSTEP_DEFAULT_ERK_6`                | Use ARKStep's default sixth-order ERK method               |
    |                                               | (ARKODE_VERNER_8_5_6).                                     |
    +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKSTEP_DEFAULT_ERK_7`                | Use ARKStep's default seventh-order ERK method             |
+   |                                               | (ARKODE_VERNER_10_6_7).                                    |
+   +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ARKSTEP_DEFAULT_ERK_8`                | Use ARKStep's default eighth-order ERK method              |
    |                                               | (ARKODE_FEHLBERG_13_7_8).                                  |
    +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKSTEP_DEFAULT_ERK_9`                | Use ARKStep's default ninth-order ERK method               |
+   |                                               | (ARKODE_VERNER_16_8_9).                                    |
+   +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ERKSTEP_DEFAULT_2`                    | Use ERKStep's default second-order ERK method              |
    |                                               | (ARKODE_HEUN_EULER_2_1_2).                                 |
    +-----------------------------------------------+------------------------------------------------------------+
@@ -139,9 +155,15 @@ contains the ARKODE output constants.
    | :index:`ERKSTEP_DEFAULT_6`                    | Use ERKStep's default sixth-order ERK method               |
    |                                               | (ARKODE_VERNER_8_5_6).                                     |
    +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ERKSTEP_DEFAULT_7`                    | Use ERKStep's default seventh-order ERK method             |
+   |                                               | (ARKODE_VERNER_10_6_7).                                    |
+   +-----------------------------------------------+------------------------------------------------------------+
    | :index:`ERKSTEP_DEFAULT_8`                    | Use ERKStep's default eighth-order ERK method              |
    |                                               | (ARKODE_FEHLBERG_13_7_8).                                  |
    +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ERKSTEP_DEFAULT_9`                    | Use ERKStep's default ninth-order ERK method               |
+   |                                               | (ARKODE_VERNER_16_8_9).                                    |
+   +-----------------------------------------------+------------------------------------------------------------+
    |                                               |                                                            |
    +-----------------------------------------------+------------------------------------------------------------+
    | **Implicit Butcher table specification**      |                                                            |

diff --git a/doc/arkode/guide/source/Introduction.rst b/doc/arkode/guide/source/Introduction.rst
@@ -64,8 +64,8 @@ for splittings tuned for use with optimal implicit solver algorithms.
 This framework allows for significant freedom over the constitutive
 methods used for each component, and ARKODE is packaged with a wide
 array of built-in methods for use.  These built-in Butcher tables
-include adaptive explicit methods of orders 2-8, adaptive implicit
-methods of orders 2-5, and adaptive ImEx methods of orders 3-5.
+include adaptive explicit methods of orders 2-9, adaptive implicit
+methods of orders 2-5, and adaptive ImEx methods of orders 2-5.
 
 
 *ERKStep* focuses specifically on problems posed in explicit form,
@@ -77,7 +77,7 @@ methods of orders 2-5, and adaptive ImEx methods of orders 3-5.
 allowing for increased computational efficiency and memory savings.
 The algorithms used in ERKStep are adaptive- and fixed-step explicit
 Runge--Kutta methods.   As with ARKStep, the ERKStep module is packaged
-with adaptive explicit methods of orders 2-8.
+with adaptive explicit methods of orders 2-9.
 
 *SPRKStep* focuses on Hamiltonian systems posed in the form,
 
@@ -156,10 +156,13 @@ Fixed scaling bug in ``SUNMatScaleAddI_Sparse`` for non-square matrices.
 Fixed missing soversions in some ``SUNLinearSolver`` and ``SUNNonlinearSolver``
 CMake targets.
 
-Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
-
-Added the fourth order ERK method ``ARKODE_SOFRONIOU_SPALETTA_5_3_4``.
+Added the third order ERK method ``ARKODE_SHU_OSHER_3_2_3``, the fourth order
+ERK method ``ARKODE_SOFRONIOU_SPALETTA_5_3_4``, the sixth order ERK method
+``ARKODE_VERNER_9_5_6``, the seventh order ERK method ``ARKODE_VERNER_10_6_7``,
+the eighth order ERK method ``ARKODE_VERNER_13_7_8``, and the ninth order ERK
+method ``ARKODE_VERNER_16_8_9``.
 
+Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
 
 Changes in v5.6.1
 -----------------

diff --git a/doc/arkode/guide/source/Mathematics.rst b/doc/arkode/guide/source/Mathematics.rst
@@ -375,10 +375,11 @@ described in the section :numref:`Butcher`.
 For mixed stiff/nonstiff problems, a user should provide both of the
 functions :math:`f^E` and :math:`f^I` that define the IVP system.  For
 such problems, ARKStep currently implements the ARK methods proposed in
-:cite:p:`KenCarp:03`, allowing for methods having order of accuracy :math:`q =
-\{3,4,5\}` and embeddings with orders :math:`p = \{2, 3, 4\}`;
-the tables for these methods are given in section :numref:`Butcher.additive`.
-Additionally, user-defined ARK tables are supported.
+:cite:p:`KenCarp:03,KenCarp:19,giraldo2013implicit`, allowing for methods having
+order of accuracy :math:`q = \{2,3,4,5\}` and embeddings with orders :math:`p =
+\{1,2,3,4\}`; the tables for these methods are given in section
+:numref:`Butcher.additive`.  Additionally, user-defined ARK tables are
+supported.
 
 For nonstiff problems, a user may specify that :math:`f^I = 0`,
 i.e. the equation :eq:`ARKODE_IMEX_IVP` reduces to the non-split IVP
@@ -391,8 +392,8 @@ In this scenario, the coefficients :math:`A^I=0`, :math:`c^I=0`,
 :math:`b^I=0` and :math:`\tilde{b}^I=0` in :eq:`ARKODE_ARK`, and the ARK
 methods reduce to classical :index:`explicit Runge--Kutta methods`
 (ERK).  For these classes of methods, ARKODE provides coefficients
-with orders of accuracy :math:`q = \{2,3,4,5,6,8\}`, with embeddings
-of orders :math:`p = \{1,2,3,4,5,7\}`.  These default to the methods in
+with orders of accuracy :math:`q = \{2,3,4,5,6,7,8,9\}`, with embeddings
+of orders :math:`p = \{1,2,3,4,5,6,7,8\}`.  These default to the methods in
 sections
 :numref:`Butcher.Heun_Euler`,
 :numref:`Butcher.Bogacki_Shampine`, :numref:`Butcher.Zonneveld`,

diff --git a/doc/shared/figs/arkode/shu_osher_erk_stab_region.png b/doc/shared/figs/arkode/shu_osher_erk_stab_region.png
diff --git a/doc/shared/figs/arkode/v65b_erk_stab_region.png b/doc/shared/figs/arkode/v65b_erk_stab_region.png
diff --git a/doc/shared/figs/arkode/v76_erk_stab_region.png b/doc/shared/figs/arkode/v76_erk_stab_region.png
diff --git a/doc/shared/figs/arkode/v87_erk_stab_region.png b/doc/shared/figs/arkode/v87_erk_stab_region.png
diff --git a/doc/shared/figs/arkode/v98_erk_stab_region.png b/doc/shared/figs/arkode/v98_erk_stab_region.png
diff --git a/doc/shared/sundials.bib b/doc/shared/sundials.bib
@@ -445,9 +445,9 @@ @book{HaWa:91
 address   = {Berlin},
 year      = 1991,
 }
-@book{HaWa:06, 
-author    = {Hairer, Ernst and Wanner, Gerhard and Lubich, Christian}, 
-title     = {{Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations}}, 
+@book{HaWa:06,
+author    = {Hairer, Ernst and Wanner, Gerhard and Lubich, Christian},
+title     = {{Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations}},
 publisher = {Springer Series in Computational Mathematics},
 doi       = {10.1007/3-540-30666-8},
 year      = {2006}
@@ -1806,6 +1806,17 @@ @article{HaSo:05
   publisher={SIAM}
 }
 
+@article{FCS:21,
+  author = {Imre Fekete and Sidafa Conde and John N. Shadid},
+  title = {Embedded pairs for optimal explicit strong stability preserving {Runge--Kutta} methods},
+  journal = {Journal of Computational and Applied Mathematics},
+  volume = {412},
+  pages = {114325},
+  year = {2022},
+  issn = {0377-0427},
+  doi = {10.1016/j.cam.2022.114325}
+}
+
 @article{FFKMS:14,
   author  = {Falgout, R.D. and Friedhoff, S. and Kolev, TZ.V. and MacLachlan, S.P. and Schroder, J.B.},
   title   = {Parallel Time Integration with Multigrid},
@@ -1998,6 +2009,18 @@ @article{Shampine:80
   doi     = {10.1137/0901005}
 }
 
+@article{ShOs:88,
+  author = {Chi-Wang Shu and Stanley Osher},
+  title = {Efficient implementation of essentially non-oscillatory shock-capturing schemes},
+  journal = {Journal of Computational Physics},
+  volume = {77},
+  number = {2},
+  pages = {439-471},
+  year = {1988},
+  issn = {0021-9991},
+  doi = {10.1016/0021-9991(88)90177-5}
+}
+
 @article{Sod:98,
   author  = {Soderlind, G.},
   title   = {The automatic control of numerical integration},
@@ -2076,6 +2099,17 @@ @article{Ver:78
   doi     = {10.1137/0715051}
 }
 
+@article{Ver:10,
+  author  = {Verner, J.H},
+  title   = {Numerically optimal {Runge–Kutta} pairs with interpolants},
+  journal = {Numerical Algorithms},
+  volume  = {53},
+  number  = {2},
+  pages   = {383-396},
+  year    = {2010},
+  doi     = {10.1007/s11075-009-9290-3}
+}
+
 @article{Ruth:93,
   title   = {A canonical integration technique},
   author  = {Ruth, Ronald D},

diff --git a/include/arkode/arkode_arkstep.h b/include/arkode/arkode_arkstep.h
@@ -41,7 +41,9 @@ static const int ARKSTEP_DEFAULT_ERK_3 = ARKODE_BOGACKI_SHAMPINE_4_2_3;
 static const int ARKSTEP_DEFAULT_ERK_4 = ARKODE_ZONNEVELD_5_3_4;
 static const int ARKSTEP_DEFAULT_ERK_5 = ARKODE_CASH_KARP_6_4_5;
 static const int ARKSTEP_DEFAULT_ERK_6 = ARKODE_VERNER_8_5_6;
+static const int ARKSTEP_DEFAULT_ERK_7 = ARKODE_VERNER_10_6_7;
 static const int ARKSTEP_DEFAULT_ERK_8 = ARKODE_FEHLBERG_13_7_8;
+static const int ARKSTEP_DEFAULT_ERK_9 = ARKODE_VERNER_16_8_9;
 
 /*    implicit */
 static const int ARKSTEP_DEFAULT_DIRK_2 = ARKODE_SDIRK_2_1_2;

diff --git a/include/arkode/arkode_butcher_erk.h b/include/arkode/arkode_butcher_erk.h
@@ -64,7 +64,7 @@ extern "C" {
 /* DEPRECATED MIN_ERK_NUM: use ARKODE_MIN_ERK_NUM */
 #define MIN_ERK_NUM              0
 /* DEPRECATED MAX_ERK_NUM: use ARKODE_MAX_ERK_NUM */
-#define MAX_ERK_NUM             14
+#define MAX_ERK_NUM             21
 
 typedef enum {
   ARKODE_ERK_NONE = -1, /* ensure enum is signed int */
@@ -86,7 +86,12 @@ typedef enum {
   ARKODE_ARK548L2SAb_ERK_8_4_5,
   ARKODE_ARK2_ERK_3_1_2,
   ARKODE_SOFRONIOU_SPALETTA_5_3_4,
-  ARKODE_MAX_ERK_NUM = ARKODE_SOFRONIOU_SPALETTA_5_3_4
+  ARKODE_SHU_OSHER_3_2_3,
+  ARKODE_VERNER_9_5_6,
+  ARKODE_VERNER_10_6_7,
+  ARKODE_VERNER_13_7_8,
+  ARKODE_VERNER_16_8_9,
+  ARKODE_MAX_ERK_NUM = ARKODE_VERNER_16_8_9
 } ARKODE_ERKTableID;
 
 /* Accessor routine to load built-in ERK table */

diff --git a/include/arkode/arkode_erkstep.h b/include/arkode/arkode_erkstep.h
@@ -36,7 +36,9 @@ static const int ERKSTEP_DEFAULT_3 = ARKODE_BOGACKI_SHAMPINE_4_2_3;
 static const int ERKSTEP_DEFAULT_4 = ARKODE_ZONNEVELD_5_3_4;
 static const int ERKSTEP_DEFAULT_5 = ARKODE_CASH_KARP_6_4_5;
 static const int ERKSTEP_DEFAULT_6 = ARKODE_VERNER_8_5_6;
+static const int ERKSTEP_DEFAULT_7 = ARKODE_VERNER_10_6_7;
 static const int ERKSTEP_DEFAULT_8 = ARKODE_FEHLBERG_13_7_8;
+static const int ERKSTEP_DEFAULT_9 = ARKODE_VERNER_16_8_9;
 
 #ifndef DEFAULT_ERK_2
 /* DEPRECATED DEFAULT_ERK_2: use ERKSTEP_DEFAULT_2 */

diff --git a/src/arkode/arkode_arkstep.c b/src/arkode/arkode_arkstep.c
@@ -2031,14 +2031,19 @@ int arkStep_SetButcherTables(ARKodeMem ark_mem)
       etable = ARKSTEP_DEFAULT_ERK_6;
       break;
     case(7):
+      etable = ARKSTEP_DEFAULT_ERK_7;
+      break;
     case(8):
       etable = ARKSTEP_DEFAULT_ERK_8;
       break;
+    case(9):
+      etable = ARKSTEP_DEFAULT_ERK_9;
+      break;
     default:    /* no available method, set default */
       arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKODE::ARKStep",
                       "arkStep_SetButcherTables",
-                      "No explicit method at requested order, using q=6.");
-      etable = ARKSTEP_DEFAULT_ERK_6;
+                      "No explicit method at requested order, using q=9.");
+      etable = ARKSTEP_DEFAULT_ERK_9;
       break;
     }