Feature: Update Full RHS Calls (#324)

Fixed a bug in ERKStep where methods with `c[s-1] = 1` but `a[s-1,j] != b[j]` were incorrectly treated as having the first same as last (FSAL) property. Fixed a bug in `MRIStepCoupling_Write` where explicit coupling tables were not written to the output file pointer. ARKStep, ERKStep, MRIStep, and SPRKStep were updated to remove a potentially unnecessary right-hand side evaluation at the end of an integration. ARKStep was additionally updated to remove extra right-hand side evaluations when using an explicit method or an implicit method with an explicit first stage. The `MRIStepInnerStepper` class in MRIStep was updated to make supplying an `MRIStepInnerFullRhsFn` optional. --------- Co-authored-by: Daniel R. Reynolds <[email protected]> Co-authored-by: Cody Balos <[email protected]>
LLNL · Dec 18, 2023 · 75e9cd2 · 75e9cd2
1 parent 671d486
commit 75e9cd2
Show file tree

Hide file tree

Showing 171 changed files with 24,262 additions and 4,277 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -59,6 +59,20 @@ the stop time. Additionally, with ARKODE, CVODE, and CVODES this fix removes an
 unnecessary interpolation of the solution at the stop time that could occur in
 this case.
 
+Fixed a bug in ERKStep where methods with `c[s-1] = 1` but `a[s-1,j] != b[j]`
+were incorrectly treated as having the first same as last (FSAL) property.
+
+Fixed a bug in `MRIStepCoupling_Write` where explicit coupling tables were not
+written to the output file pointer.
+
+ARKStep, ERKStep, MRIStep, and SPRKStep were updated to remove a potentially
+unnecessary right-hand side evaluation at the end of an integration. ARKStep was
+additionally updated to remove extra right-hand side evaluations when using an
+explicit method or an implicit method with an explicit first stage.
+
+The `MRIStepInnerStepper` class in MRIStep was updated to make supplying an
+`MRIStepInnerFullRhsFn` optional.
+
 ## Changes to SUNDIALS in release 6.6.0
 
 A new time-stepping module, `SPRKStep`, was added to ARKODE. This time-stepper

diff --git a/doc/README.md b/doc/README.md
@@ -0,0 +1,10 @@
+# SUNDIALS Documentation
+
+The SUNDIALS documentation is written using reStructuredText and
+[Sphinx](https://www.sphinx-doc.org/).
+
+To build the documentation with Sphinx you will need Python 3.9+. Sphinx and the
+necessary extensions can be installed using the requirements file i.e.,
+`pip install -r requirements.txt`. Additionally, building the developer
+documentation requires [Graphviz](https://graphviz.org/) for generating
+flowcharts.
diff --git a/doc/arkode/guide/source/ARKodeButcherTable.rst b/doc/arkode/guide/source/ARKodeButcherTable.rst
@@ -76,33 +76,35 @@ ARKodeButcherTable functions
 .. _ARKodeButcherTable.FunctionsTable:
 .. table:: ARKodeButcherTable functions
 
-   +----------------------------------------------+------------------------------------------------------------+
-   | **Function name**                            | **Description**                                            |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_LoadERK()`       | Retrieve a given explicit Butcher table by its unique ID   |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_LoadERKByName()` | Retrieve a given explicit Butcher table by its unique name |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_LoadDIRK()`      | Retrieve a given implicit Butcher table by its unique ID   |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_LoadDIRKByName()`| Retrieve a given implicit Butcher table by its unique name |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_Alloc()`         | Allocate an empty Butcher table                            |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_Create()`        | Create a new Butcher table                                 |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_Copy()`          | Create a copy of a Butcher table                           |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_Space()`         | Get the Butcher table real and integer workspace size      |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_Free()`          | Deallocate a Butcher table                                 |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_Write()`         | Write the Butcher table to an output file                  |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_CheckOrder()`    | Check the order of a Butcher table                         |
-   +----------------------------------------------+------------------------------------------------------------+
-   | :c:func:`ARKodeButcherTable_CheckARKOrder()` | Check the order of an ARK pair of Butcher tables           |
-   +----------------------------------------------+------------------------------------------------------------+
+   +--------------------------------------------------+------------------------------------------------------------+
+   | **Function name**                                | **Description**                                            |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_LoadERK()`           | Retrieve a given explicit Butcher table by its unique ID   |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_LoadERKByName()`     | Retrieve a given explicit Butcher table by its unique name |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_LoadDIRK()`          | Retrieve a given implicit Butcher table by its unique ID   |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_LoadDIRKByName()`    | Retrieve a given implicit Butcher table by its unique name |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_Alloc()`             | Allocate an empty Butcher table                            |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_Create()`            | Create a new Butcher table                                 |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_Copy()`              | Create a copy of a Butcher table                           |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_Space()`             | Get the Butcher table real and integer workspace size      |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_Free()`              | Deallocate a Butcher table                                 |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_Write()`             | Write the Butcher table to an output file                  |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_IsStifflyAccurate()` | Determine if ``A[stages - 1][i] == b[i]``                  |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_CheckOrder()`        | Check the order of a Butcher table                         |
+   +--------------------------------------------------+------------------------------------------------------------+
+   | :c:func:`ARKodeButcherTable_CheckARKOrder()`     | Check the order of an ARK pair of Butcher tables           |
+   +--------------------------------------------------+------------------------------------------------------------+
 
 .. c:function:: ARKodeButcherTable ARKodeButcherTable_LoadERK(ARKODE_ERKTableID emethod)
 
@@ -252,6 +254,19 @@ ARKodeButcherTable functions
       The *outfile* argument can be ``stdout`` or ``stderr``, or it
       may point to a specific file created using ``fopen``.
 
+.. c:function:: void ARKodeButcherTable_IsStifflyAccurate(ARKodeButcherTable B)
+
+   Determine if the table satisfies ``A[stages - 1][i] == b[i]``
+
+   **Arguments:**
+      * *B* -- the Butcher table.
+
+   **Returns**
+      * ``SUNTRUE`` if the method is "stiffly accurate", otherwise returns
+        ``SUNFALSE``
+
+   .. versionadded:: vX.X.X
+
 .. c:function:: int ARKodeButcherTable_CheckOrder(ARKodeButcherTable B, int* q, int* p, FILE* outfile)
 
    Determine the analytic order of accuracy for the specified Butcher

diff --git a/doc/arkode/guide/source/Introduction.rst b/doc/arkode/guide/source/Introduction.rst
@@ -188,6 +188,21 @@ requested output time is the same as the stop time. Additionally, this fix
 removes an unnecessary interpolation of the solution at the stop time that could
 occur in this case.
 
+Fixed a bug in ERKStep where methods with :math:`c_s = 1` but
+:math:`a_{s,j} \neq b_j` were incorrectly treated as having the first same as
+last (FSAL) property.
+
+Fixed a bug in :c:func:`MRIStepCoupling_Write` where explicit coupling tables
+were not written to the output file pointer.
+
+ARKStep, ERKStep, MRIStep, and SPRKStep were updated to remove a potentially
+unnecessary right-hand side evaluation at the end of an integration. ARKStep was
+additionally updated to remove extra right-hand side evaluations when using an
+explicit method or an implicit method with an explicit first stage.
+
+The :c:type:`MRIStepInnerStepper` class in MRIStep was updated to make supplying
+an :c:func:`MRIStepInnerFullRhsFn` optional.
+
 Changes in v5.6.0
 -----------------
 

diff --git a/doc/arkode/guide/source/Mathematics.rst b/doc/arkode/guide/source/Mathematics.rst
@@ -466,39 +466,41 @@ than the more general form :eq:`ARKODE_IVP_simple_explicit`.
 SPRKStep -- Symplectic Partitioned Runge--Kutta methods
 =======================================================
 
-The SPRKStep time-stepping module in ARKODE is designed for IVPs of the form
+The SPRKStep time-stepping module in ARKODE is designed for problems where the
+state vector is partitioned as
 
 .. math::
-   \dot{p} = f_1(t, q) = \frac{\partial V(t, q)}{\partial q}, \quad
-   \dot{q} = f_2(t, p) = \frac{\partial T(t, p)}{\partial p},
-   \qquad p(t_0) = p_0,\quad q(t_0) = q_0,
-   :label: ARKODE_IVP_Hamiltonian
+   y(t) =
+   \begin{bmatrix}
+     p(t) \\
+     q(t)
+   \end{bmatrix}
 
-where the system Hamiltonian
+and the component partitioned IVP is given by
 
 .. math::
-   H(t, p, q) = T(t, p) + V(t, q)
+   \dot{p} &= f_1(t, q), \qquad p(t_0) = p_0 \\
+   \dot{q} &= f_2(t, p), \qquad q(t_0) = q_0.
+   :label: ARKODE_IVP_SPRK
 
-**is separable**. When *H* is autonomous, then *H* is a conserved quantity.
-Often this correponds to the conservation of energy (for example, in *n*-body
-problems). For non-autonomous *H*, the invariants are no longer directly
-obtainable from the Hamiltonian :cite:p:`Struckmeier:02`.
+The right-hand side functions :math:`f_1(t,p)` and :math:`f_2(t,q)` typically
+arise from the **separable** Hamiltonian system
 
-In solving the IVP :eq:`ARKODE_IVP_Hamiltonian`, we consider the problem in the form
+.. math::
+   H(t, p, q) = T(t, p) + V(t, q)
+
+where
 
 .. math::
-   \dot{y} =
-      \begin{bmatrix}
-         f_1(t, q) \\
-         f_2(t, p)
-      \end{bmatrix}, \qquad
-      y(t_0) =
-      \begin{bmatrix}
-         p_0\\
-         q_0
-      \end{bmatrix}.
+   f_1(t, q) \equiv \frac{\partial V(t, q)}{\partial q}, \qquad
+   f_2(t, p) \equiv \frac{\partial T(t, p)}{\partial p}.
+
+When *H* is autonomous, then *H* is a conserved quantity. Often this corresponds
+to the conservation of energy (for example, in *n*-body problems). For
+non-autonomous *H*, the invariants are no longer directly obtainable from the
+Hamiltonian :cite:p:`Struckmeier:02`.
 
-In practice, the ordering  of the variables does not matter and is determined by the user.
+In practice, the ordering of the variables does not matter and is determined by the user.
 SPRKStep utilizes Symplectic Partitioned Runge-Kutta (SPRK) methods represented by the pair
 of explicit and diagonally implicit Butcher tableaux,
 
@@ -529,17 +531,18 @@ schemes with order of accuracy and conservation equal to
 the default methods used are given in the section :numref:`Butcher.sprk`.
 
 In the default case, the algorithm for a single time-step is as follows
-(for autonomous Hamiltonian systems the times provided to :math:`f1` and :math:`f2`
+(for autonomous Hamiltonian systems the times provided to :math:`f_1` and
+:math:`f_2`
 can be ignored).
 
-#. Set :math:`P_0 = p_n, Q_1 = q_n`
+#. Set :math:`P_0 = p_{n-1}, Q_1 = q_{n-1}`
 
 #. For :math:`i = 1,\ldots,s` do:
 
-   #. :math:`P_i = P_{i-1} + h_{n+1} \hat{a}_i f_1(t_n + \hat{c}_i h, Q_i)`
-   #. :math:`Q_{i+1} = Q_i + h_{n+1} a_i f_2(t_n + c_i h, P_i)`
+   #. :math:`P_i = P_{i-1} + h_n \hat{a}_i f_1(t_{n-1} + \hat{c}_i h_n, Q_i)`
+   #. :math:`Q_{i+1} = Q_i + h_n a_i f_2(t_{n-1} + c_i h_n, P_i)`
 
-#. Set :math:`p_{n+1} = P_s, q_{n+1} = Q_{s+1}`
+#. Set :math:`p_n = P_s, q_n = Q_{s+1}`
 
 .. _ARKODE.Mathematics.SPRKStep.Compensated:
 
@@ -554,12 +557,12 @@ form is used to compute a time step:
 
 #. For :math:`i = 1,\ldots,s` do:
 
-   #. :math:`\Delta P_i = \Delta P_{i-1} + h_{n+1} \hat{a}_i f_1(t_n + \hat{c}_i h, q_n + \Delta Q_i)`
-   #. :math:`\Delta Q_{i+1} = \Delta Q_i + h_{n+1} a_i f_2(t_n + c_i h, p_n + \Delta P_i)`
+   #. :math:`\Delta P_i = \Delta P_{i-1} + h_n \hat{a}_i f_1(t_{n-1} + \hat{c}_i h_n, q_{n-1} + \Delta Q_i)`
+   #. :math:`\Delta Q_{i+1} = \Delta Q_i + h_n a_i f_2(t_{n-1} + c_i h_n, p_{n-1} + \Delta P_i)`
 
-#. Set :math:`\Delta p_{n+1} = \Delta P_s, \Delta q_{n+1} = \Delta Q_{s+1}`
+#. Set :math:`\Delta p_n = \Delta P_s, \Delta q_n = \Delta Q_{s+1}`
 
-#. Using compensated summation, set :math:`p_{n+1} = p_n + \Delta p_{n+1}, q_{n+1} = q_n + \Delta Q_{s+1}`
+#. Using compensated summation, set :math:`p_n = p_{n-1} + \Delta p_n, q_n = q_{n-1} + \Delta q_n`
 
 Since temporal error based adaptive time-stepping is known to ruin the
 conservation property :cite:p:`HaWa:06`,  SPRKStep employs a fixed time-step size.
@@ -706,7 +709,7 @@ characterized by nonzero values on or above the diagonal of the matrices
 :math:`\Gamma^{\{k\}}`. Typically, MRI-GARK and IMEX-MRI-GARK methods are at
 most diagonally-implicit (i.e., :math:`\gamma_{i,j}^{\{k\}}=0` for all
 :math:`j>i`). Furthermore, diagonally-implicit stages are characterized as being
-"solve-decoupled" if :math:`\Delta c_i^S = 0` when `\gamma_{i,i}^{\{k\}} \ne 0`,
+"solve-decoupled" if :math:`\Delta c_i^S = 0` when :math:`\gamma_{i,i}^{\{k\}} \ne 0`,
 in which case the stage is computed as standard ARK or DIRK update. Alternately,
 a diagonally-implicit stage :math:`i` is considered "solve-coupled" if
 :math:`\Delta c^S_i \gamma_{i,j}^{\{k\}} \ne 0`, in which

diff --git a/...ode/guide/source/Usage/MRIStep_c_interface/Custom_Inner_Stepper/Description.rst b/...ode/guide/source/Usage/MRIStep_c_interface/Custom_Inner_Stepper/Description.rst
@@ -365,6 +365,11 @@ member functions:
    **Example codes:**
       * ``examples/arkode/CXX_parallel/ark_diffusion_reaction_p.cpp``
 
+Optional Member Functions
+"""""""""""""""""""""""""
+
+An :c:type:`MRIStepInnerStepper` *may* provide implementations of any of the
+following member functions:
 
 .. c:type:: int (*MRIStepInnerFullRhsFn)(MRIStepInnerStepper stepper, realtype t, N_Vector v, N_Vector f, int mode)
 
@@ -393,11 +398,9 @@ member functions:
    **Example codes:**
       * ``examples/arkode/CXX_parallel/ark_diffusion_reaction_p.cpp``
 
-Optional Member Functions
-"""""""""""""""""""""""""
+   .. versionchanged:: vX.X.X
 
-An :c:type:`MRIStepInnerStepper` *may* provide implementations of any of the
-following member functions:
+      Supplying a full right-hand side function was made optional.
 
 .. c:type:: int (*MRIStepInnerResetFn)(MRIStepInnerStepper stepper, realtype tR, N_Vector vR)
 

diff --git a/doc/arkode/guide/source/Usage/SPRKStep_c_interface/index.rst b/doc/arkode/guide/source/Usage/SPRKStep_c_interface/index.rst
@@ -17,19 +17,19 @@ Using the SPRKStep time-stepping module
 ==========================================
 
 This chapter is concerned with the use of the SPRKStep time-stepping module for
-the solution of Hamiltonian initial value problems (IVPs) of the form
-:eq:`ARKODE_IVP_Hamiltonian` in a C or C++ language setting. The following sections
+the solution of initial value problems (IVPs) of the form
+:eq:`ARKODE_IVP_SPRK` in a C or C++ language setting. The following sections
 discuss the header files and the layout of the user's main program, and provide
 descriptions of the SPRKStep user-callable functions and user-supplied functions.
 
-The example programs located in the source code ``examples/arkode`` folder, may 
+The example programs located in the source code ``examples/arkode`` folder, may
 be helpful as templates for new codes. In particular,
 
 * ``examples/arkode/C_serial/ark_harmonic_symplectic.c``
 * ``examples/arkode/C_serial/ark_damped_harmonic_symplectic.c``, and
 * ``examples/arkode/C_serial/ark_kepler.c``
 
-demonstrate ``SPRKStep`` usage. 
+demonstrate ``SPRKStep`` usage.
 
 SPRKStep uses the input and output constants from the shared ARKODE infrastructure.
 These are defined as needed in this chapter, but for convenience the full list is

diff --git a/doc/sundials_developers/source/conf.py b/doc/sundials_developers/source/conf.py
@@ -34,7 +34,8 @@ def setup(app: Sphinx):
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
 extensions = ['sphinx_rtd_theme', 'sphinx.ext.ifconfig', 'sphinx.ext.mathjax',
-              'sphinxfortran.fortran_domain', 'sphinxcontrib.bibtex', 'sphinx_copybutton']
+              'sphinxfortran.fortran_domain', 'sphinxcontrib.bibtex',
+              'sphinx_copybutton', 'sphinx.ext.graphviz']
 
 # References
 bibtex_bibfiles = ['../../shared/sundials.bib']

diff --git a/doc/sundials_developers/source/index.rst b/doc/sundials_developers/source/index.rst
@@ -28,5 +28,6 @@ SUNDIALS Developers Guide
    benchmarks/index
    pull_requests/index
    releases/index
+   packages/index
    appendix/index
    references