Feature/operator splitting

CHANGELOG.md

@@ -4,8 +4,15 @@
 
 ### Major Features
 
+Added an operator splitting module, SplittingStep, and forcing method module,
+ForcingStep, to ARKODE. These modules support a broad range of operator-split
+time integration methods for multiphysics applications.
+
 ### New Features and Enhancements
 
+Added the `ARKodeSetStepDirection` and `ARKodeGetStepDirection` functions to
+change and query the direction of integration.
+
 Added the `SUNStepper` base class to represent a generic solution procedure for
 IVPs. This is used by the SplittingStep and ForcingStep modules of ARKODE. A
 SUNStepper can be created from an ARKstep memory block with the new function

doc/arkode/guide/source/Constants.rst

@@ -317,6 +317,30 @@ contains the ARKODE output constants.
    +-----------------------------------------------+------------------------------------------------------------+
    |                                               |                                                            |
    +-----------------------------------------------+------------------------------------------------------------+
+   | **Splitting Coefficients specification**      |                                                            |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SPLITTING_LIE_TROTTER_1_1_2`   | 1st order Lie-Trotter splitting for problems with two      |
+   |                                               | partitions.                                                |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SPLITTING_STRANG_2_2_2`        | 2nd order Strang splitting for problems with two           |
+   |                                               | partitions.                                                |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SPLITTING_BEST_2_2_2`          | 2nd order splitting with optimal error for problems with   |
+   |                                               | two partitions.                                            |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SPLITTING_SUZUKI_3_3_2`        | 3rd order Suzuki splitting for problems with two           |
+   |                                               | partitions.                                                |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SPLITTING_RUTH_3_3_2`          | 3rd order Ruth splitting for problems with two partitions. |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SPLITTING_YOSHIDA_4_4_2`       | 4th order Yoshida splitting for problems with two          |
+   |                                               | partitions.                                                |
+   +-----------------------------------------------+------------------------------------------------------------+
+   | :index:`ARKODE_SPLITTING_YOSHIDA_8_6_2`       | 6th order Yoshida splitting for problems with two          |
+   |                                               | partitions.                                                |
+   +-----------------------------------------------+------------------------------------------------------------+
+   |                                               |                                                            |
+   +-----------------------------------------------+------------------------------------------------------------+
    | **MRI method types**                          |                                                            |
    +-----------------------------------------------+------------------------------------------------------------+
    | :index:`MRISTEP_EXPLICIT`                     | Use an explicit (at the slow time scale) MRI method.       |

doc/arkode/guide/source/Mathematics.rst

@@ -61,13 +61,15 @@ for interpolated solution output.  We then discuss the current suite
 of time-stepping modules supplied with ARKODE, including the ARKStep
 module for :ref:`additive Runge--Kutta methods <ARKODE.Mathematics.ARK>`,
 the ERKStep module that is optimized for :ref:`explicit Runge--Kutta
-methods <ARKODE.Mathematics.ERK>`, and the MRIStep module for :ref:`multirate
-infinitesimal step (MIS), multirate infinitesimal GARK (MRI-GARK), and
-implicit-explicit MRI-GARK (IMEX-MRI-GARK) methods <ARKODE.Mathematics.MRIStep>`.
-We then discuss the :ref:`adaptive temporal error controllers
-<ARKODE.Mathematics.Adaptivity>` shared by the time-stepping modules, including
-discussion of our choice of norms for measuring errors within various components
-of the solver.
+methods <ARKODE.Mathematics.ERK>`, SplittingStep for :ref:`operator splitting
+methods <ARKODE.Mathematics.SplittingStep>`, ForcingStep for :ref:`a forcing
+method <ARKODE.Mathematics.ForcingStep>`, and the MRIStep module for
+:ref:`multirate infinitesimal step (MIS), multirate infinitesimal GARK
+(MRI-GARK), and implicit-explicit MRI-GARK (IMEX-MRI-GARK) methods
+<ARKODE.Mathematics.MRIStep>`. We then discuss the :ref:`adaptive temporal error
+controllers <ARKODE.Mathematics.Adaptivity>` shared by the time-stepping
+modules, including discussion of our choice of norms for measuring errors within
+various components of the solver.
 
 We then discuss the nonlinear and linear solver strategies used by
 ARKODE for solving implicit algebraic systems that arise in computing each
@@ -572,6 +574,110 @@ using a fixed time-step size.
 .. The `ark_kepler.c` example demonstrates an implementation of such controller.
 
 
+.. _ARKODE.Mathematics.SplittingStep:
+
+SplittingStep -- Operator splitting methods
+================================================
+
+The SplittingStep time-stepping module in ARKODE is designed for IVPs of the
+form
+
+.. math::
+   \dot{y} = f_1(t,y) + f_2(t,y) + \dots + f_P(t,y), \qquad y(t_0) = y_0,
+
+with :math:`P > 1` additive partitions. Operator splitting methods, such as
+those implemented in SplittingStep, allow each partition to be integrated
+separately, possibly with different numerical integrators or exact solution
+procedures. Coupling is only performed though initial conditions which are
+passed from the flow of one partition to the next.
+
+The following algorithmic procedure is used in the Splitting-Step module:
+
+#. For :math:`i = 1, \dots, r` do:
+
+   #. Set :math:`y_{n, i} = y_{n - 1}`.
+
+   #. For :math:`j = 0, \dots, s` do:
+
+      #. For :math:`k = 1, \dots, P` do:
+
+         #. Let :math:`t_{\text{start}} = t_{n-1} + \beta_{i,j-1,k} h_n` and
+            :math:`t_{\text{end}} = t_{n-1} + \beta_{i,j,k} h_n`.
+
+         #. Let :math:`v(t_{\text{start}}) = y_{n,i}`.
+
+         #. For :math:`t \in [t_{\text{start}}, t_{\text{end}}]` solve
+            :math:`\dot{v} = f_{k}(t, v)`.
+
+         #. Set :math:`y_{n, i} = v(t_{\text{end}})`.
+
+#. Set :math:`y_n = \sum_{i=1}^r \alpha_i y_{n,i}`
+
+Here, :math:`s` denotes the number of stages, while :math:`r` denotes the number
+of sequential methods within the overall operator splitting scheme. The
+sequential methods have independent flows which are linearly combined to produce
+the next step. The real coefficients :math:`\alpha_i` and :math:`\beta_{i,j,k}`
+determine the particular scheme and properties such as the order of accuracy.
+
+An alternative representation of the SplittingStep solution is
+
+.. math::
+   y_n = \sum_{i=1}^P \alpha_i \left( \phi^P_{\gamma_{i,1,P} h} \circ
+   \phi^{P-1}_{\gamma_{i,1,P-1} h} \circ \dots \circ \phi^1_{\gamma_{i,1,1} h}
+   \circ \phi^P_{\gamma_{i,2,P} h} \circ \dots \circ \phi^P_{\gamma_{i,s,P} h}
+   \circ \dots \circ \phi^1_{\gamma_{i,s,1} h} \right),
+   (y_{n-1})
+
+where :math:`\gamma_{i,j,k} = \beta_{i,j,k} - \beta_{i,j-1,k}` and
+:math:`\phi^k_{h}` is the flow map for partition :math:`k`:
+
+.. math::
+   \phi^k_{h}(y_{n_1}) = v(t_n),
+   \quad \begin{cases}
+      v(t_{n-1}) = y_{n-1}, \\ \dot{v} = f_k(t, v).
+   \end{cases}
+
+SplittingStep provides standard operator splitting methods such as Lie-Trotter
+and Strang splitting, as well as schemes of arbitrarily high order.
+Alternatively, users may construct their own coefficients (see
+:numref:`ARKODE.Usage.SplittingStep.SplittingStepCoefficients`). Generally,
+methods of order three and higher with real coefficients require backward
+integration, i.e., there exist negative :math:`\gamma_{i,j,k}` coefficients.
+Currently, a fixed time step must be specified for the outer SplittingStep, but
+inner integrators are free to use adaptive time steps.
+
+
+.. _ARKODE.Mathematics.ForcingStep:
+
+ForcingStep -- Forcing method
+=============================
+
+The ForcingStep time-stepping module in ARKODE is designed for IVPs of the form
+
+.. math::
+   \dot{y} = f_1(t,y) + f_2(t,y), \qquad y(t_0) = y_0,
+
+with two additive partitions. One step of the forcing method implemented in
+ForcingStep is given by
+
+.. math::
+   v_1(t_{n-1}) &= y_{n-1}, \\
+   \dot{v}_1 &= f_1(t, v_1), \\
+   f_1^* &= \frac{v_1(t_n) - y_{n-1}}{h_n}, \\
+   v_2(t_{n-1}) &= y_{n-1}, \\
+   \dot{v}_2 &= f_1^* + f_2(t, v_2), \\
+   y_n &= v_2(t_n).
+
+Like a Lie-Trotter method from SplittingStep, the partitions are evolved through
+a sequence of inner IVPs which can be solved with an arbitrary integrator or
+exact solution procedure. However, the IVP for partition two includes a
+"forcing" or "tendency" term :math:`f_1^*` to strengthen the coupling. This
+coupling leads to a first order method provided :math:`v_1` and :math:`v_2` are
+integrated to at least first order accuracy. Currently, a fixed time step must
+be specified for the outer ForcingStep, but inner integrators are free to use
+adaptive time steps.
+
+
 .. _ARKODE.Mathematics.MRIStep:
 
 MRIStep -- Multirate infinitesimal step methods

doc/arkode/guide/source/Usage/ForcingStep/User_callable.rst

@@ -0,0 +1,114 @@
+.. ----------------------------------------------------------------
+   Programmer(s): Steven B. Roberts @ LLNL
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _ARKODE.Usage.ForcingStep.UserCallable:
+
+ForcingStep User-callable functions
+===================================
+
+This section describes the ForcingStep-specific functions that may be called
+by the user to setup and then solve an IVP using the ForcingStep time-stepping
+module.
+
+As discussed in the main :ref:`ARKODE user-callable function introduction
+<ARKODE.Usage.UserCallable>`, each of ARKODE's time-stepping modules
+clarifies the categories of user-callable functions that it supports.
+ForcingStep does not support any of the categories beyond the functions that
+apply for all time-stepping modules.
+
+
+.. _ARKODE.Usage.ForcingStep.Initialization:
+
+ForcingStep initialization functions
+------------------------------------
+
+.. c:function:: void* ForcingStepCreate(SUNStepper stepper1, SUNStepper stepper2, sunrealtype t0, N_Vector y0, SUNContext sunctx)
+
+   This function allocates and initializes memory for a problem to be solved
+   using the ForcingStep time-stepping module in ARKODE.
+
+   **Arguments:**
+      * *stepper1* -- a :c:type:`SUNStepper` to integrate partition one.
+      * *stepper2* -- a :c:type:`SUNStepper` to integrate partition two
+        including the forcing from partition one.
+      * *t0* -- the initial value of :math:`t`.
+      * *y0* -- the initial condition vector :math:`y(t_0)`.
+      * *sunctx* -- the :c:type:`SUNContext` object (see :numref:`SUNDIALS.SUNContext`)
+
+   **Return value:**
+      If successful, a pointer to initialized problem memory of type ``void*``,
+      to be passed to all user-facing ForcingStep routines listed below. If
+      unsuccessful, a ``NULL`` pointer will be returned, and an error message
+      will be printed to ``stderr``.
+
+   **Example usage:**
+
+      .. code-block:: C
+
+         /* inner ARKODE objects for integrating individual partitions */
+         void *partition_mem[] = {NULL, NULL};
+
+         /* SUNSteppers to wrap the inner ARKStep objects */
+         SUNStepper steppers[] = {NULL, NULL};
+
+         /* create ARKStep objects, setting right-hand side functions and the
+            initial condition */
+         partition_mem[0] = ARKStepCreate(fe1, fi1, t0, y0, sunctx);
+         partition_mem[1] = ARKStepCreate(fe2, fi2, t0, y0, sunctx);
+
+         /* setup ARKStep */
+         . . .
+
+         /* create SUNStepper wrappers for the ARKStep memory blocks */
+         flag = ARKStepCreateSUNStepper(partition_mem[0], &stepper[0]);
+         flag = ARKStepCreateSUNStepper(partition_mem[1], &stepper[1]);
+
+         /* create a ForcingStep object */
+         arkode_mem = ForcingStepCreate(steppers[0], steppers[1], t0, y0, sunctx);
+
+   **Example codes:**
+      * ``examples/arkode/C_serial/ark_analytic_partitioned.c``
+
+   .. versionadded:: x.y.z
+
+
+.. _ARKODE.Usage.ForcingStep.OptionalOutputs:
+
+
+Optional output functions
+------------------------------
+
+.. c:function:: int ForcingStep_GetNumEvolves(void* arkode_mem, int partition, long int *evolves)
+
+   Returns the number of times the :c:type:`SUNStepper` for the given partition
+   index has been evolved (so far).
+
+   **Arguments:**
+
+   * *arkode_mem* -- pointer to the ForcingStep memory block.
+
+   * *partition* -- index of the partition (0 or 1) or a negative number to
+     indicate the total number across both partitions.
+
+   * *evolves* -- number of :c:type:`SUNStepper` evolves.
+
+   **Return value:**
+
+   * *ARK_SUCCESS* if successful
+
+   * *ARK_MEM_NULL* if the ForcingStep memory was ``NULL``
+
+   * *ARK_ILL_INPUT* if *partition* was out of bounds
+
+   .. versionadded:: x.y.z

doc/arkode/guide/source/Usage/ForcingStep/index.rst

@@ -0,0 +1,32 @@
+.. ----------------------------------------------------------------
+   Programmer(s): Steven B. Roberts @ LLNL
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _ARKODE.Usage.ForcingStep:
+
+==========================================
+Using the ForcingStep time-stepping module
+==========================================
+
+This section is concerned with the use of the ForcingStep time-stepping module
+for the solution of initial value problems (IVPs) in a C or C++ language
+setting.  Usage of ForcingStep follows that of the rest of ARKODE, and so in
+this section we primarily focus on those usage aspects that are specific to
+ForcingStep. A skeleton of a program using ForcingStep follows essentially the
+same structure as SplittingStep
+(see :numref:`ARKODE.Usage.SplittingStep.Skeleton`).
+
+.. toctree::
+   :maxdepth: 1
+
+   User_callable

doc/arkode/guide/source/Usage/General.rst

@@ -51,11 +51,13 @@ to the SUNDIALS core header file.
 
 .. code:: c
 
-  #include <sundials/sundials_core.h> // Provides core SUNDIALS types
-  #include <arkode/arkode_erkstep.h>  // ERKStep provides explicit RK methods.
-  #include <arkode/arkode_arkstep.h>  // ARKStep provides explicit, implicit, IMEX additive RK methods.
-  #include <arkode/arkode_mristep.h>  // MRIStep provides mutlirate RK methods.
-  #include <arkode/arkode_sprkstep.h> // SPRKStep provides symplectic partition RK methods.
+  #include <sundials/sundials_core.h>      // Provides core SUNDIALS types
+  #include <arkode/arkode_arkstep.h>       // ARKStep provides explicit, implicit, IMEX additive RK methods.
+  #include <arkode/arkode_erkstep.h>       // ERKStep provides explicit RK methods.
+  #include <arkode/arkode_sprkstep.h>      // SPRKStep provides symplectic partitioned RK methods.
+  #include <arkode/arkode_splittingstep.h> // SplittingStep provides operator splitting methods.
+  #include <arkode/arkode_forcingstep.h>   // ForcingStep provides a forcing method.
+  #include <arkode/arkode_mristep.h>       // MRIStep provides multirate RK methods.
 
 Each of these define several types and various constants, include function
 prototypes, and include the shared ``arkode/arkode.h`` and