merged with upstream changes

doc/arkode/guide/source/Mathematics.rst

@@ -641,9 +641,11 @@ algorithm for a single step:
    diagonally-implicit, or additive Runge--Kutta stage update,
 
    .. math::
-      z_i - \theta_{i,i} h^S f^I(t_{n-1}+c_i^S h^S, z_i) = a_i.
+      z_i - \theta_{i,i} h^S f^I(t_{n,i}^S, z_i) = a_i.
       :label: MRI_implicit_solve
 
+   where :math:`t_{n,j}^S = t_{n-1}*h^S c^S_j`.
+
 #. Set :math:`y_{n} = z_{s}`.
 
 #. If the method has an embedding, compute the embedded solution,
@@ -682,19 +684,19 @@ MIS, MRI-GARK, and IMEX-MRI-GARK Methods
 
 The methods in IMEX-MRI-GARK family, which includes MIS and MRI-GARK methods,
 are defined by a vector of slow stage time abscissae, :math:`c^S \in
-\mathbb{R}^{s+1}`, and a set of coupling tensors,
-:math:`\Omega\in\mathbb{R}^{(s+1)\times(s+1)\times k}` and
-:math:`\Gamma\in\mathbb{R}^{(s+1)\times(s+1)\times k}`, that specify the
+\mathbb{R}^{s}`, and a set of coupling tensors,
+:math:`\Omega\in\mathbb{R}^{(s+1)\times s \times k}` and
+:math:`\Gamma\in\mathbb{R}^{(s+1)\times s \times k}`, that specify the
 slow-to-fast coupling for the explicit and implicit components, respectively.
 
 The fast stage IVPs, :eq:`MRI_fast_IVP`, are evolved over non-overlapping
-intervals :math:`[t_{0,i},t_{F,i}] = [t_{n-1}+c_{i-1}h^S, t_{n-1}+c_ih^S]` with
+intervals :math:`[t_{0,i},t_{F,i}] = [t_{n,i-1}^S, t_{n,i}^S]` with
 the initial condition :math:`v_{0,i}=z_{i-1}`. The fast IVP forcing function is
 given by
 
 .. math::
-   r_i(t) = \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^{i-1} \omega_{i,j}(\tau) f^E(t_{n,j}^I, z_j)
-   + \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^i \gamma_{i,j}(\tau) f^I(t_{n,j}^I, z_j)
+   r_i(t) = \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^{i-1} \omega_{i,j}(\tau) f^E(t_{n,j}^S, z_j)
+   + \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^i \gamma_{i,j}(\tau) f^I(t_{n,j}^S, z_j)
 
 where :math:`\Delta c_i^S=\left(c^S_i - c^S_{i-1}\right)`, :math:`\tau = (t -
 t_{n,i-1}^S)/(h^S \Delta c_i^S)` is the normalized time, the coefficients
@@ -721,7 +723,7 @@ stage is computed as
    :label: ARKODE_MRI_delta_c_zero
 
 Similarly, the embedded solution IVP, :eq:`MRI_embedding_fast_IVP`, is evolved
-over the interval :math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n-1}+c_{s-1}h^S, t_{n}]`
+over the interval :math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n,s-1}^S, t_{n}]`
 with the initial condition :math:`\tilde{v}_0=z_{s-1}`.
 
 As with standard ARK and DIRK methods, implicitness at the slow time scale is
@@ -746,26 +748,87 @@ The IMEX-MRI-SR family of methods perform *both* the fast IVP evolution,
 :eq:`MRI_fast_IVP` or :eq:`MRI_embedding_fast_IVP`, *and* stage update,
 :eq:`MRI_implicit_solve` or :eq:`MRI_embedding_implicit_solve`, in every stage
 (but these methods typically have far fewer stages than implicit MRI-GARK or
-IMEX-MRI-GARK methods).
+IMEX-MRI-GARK methods).  These methods are defined by a vector of slow stage
+time abscissae :math:`c^S \in \mathbb{R}^{s}`, a set of coupling tensors
+:math:`\Omega\in\mathbb{R}^{(s+1)\times s\times k}`, and a Butcher table of
+slow-implicit coefficients, :math:`\Gamma\in\mathbb{R}^{(s+1) \times s}`.
+
+The fast stage IVPs, :eq:`MRI_fast_IVP`, are evolved on overlapping
+intervals :math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n,i}^S]` with
+the initial condition :math:`v_{0,i}=y_{n-1}`. The fast IVP forcing function is
+given by
+
+.. math::
+   r_i(t) = \frac{1}{c_i^S} \sum\limits_{j=1}^{i-1} \omega_{i,j}(\tau) \left( f^E(t_{n,j}^S, z_j) + f^I(t_{n,j}^S, z_j)\right),
+   :label: IMEXMRISR_forcing
+
+where :math:`\tau = (t - t_n)/(h^S c_i^S)` is the normalized time, and the coefficients
+:math:`\omega_{i,j}` are polynomials in time of degree :math:`k` that are also given by
+:eq:`ARKODE_MRI_coupling`.  The solution of these fast IVPs defines an intermediate stage
+solution, :math:`\tilde{z}_i`.
+
+The implicit solve that follows each fast IVP must solve the algebraic equation for :math:`z_i`
+
+.. math::
+   z_i = \tilde{z}_i + h^S \sum_{j=1}^{i} \gamma_{i,j} f^I(t_{n,j}^S, z_j).
+   :label: ARKODE_MRISR_implicit
+
+We note that IMEX-MRI-SR methods are solve-decoupled by construction, and thus the structure
+of a given stage never needs to be deduced based on :math:`\Delta c_i^S`.  However, ARKODE
+still checks the value of :math:`\gamma_{i,i}`, since if it zero then the stage update
+equation :eq:`ARKODE_MRISR_implicit` simplifies to a simple explicit Runge--Kutta-like stage
+update.
+
+The overall time step solution is given by the final internal stage solution,
+i.e., :math:`y_{n} = z_{s}`.  The embedded solution is computing using the above
+algorithm for stage index :math:`s+1`, under the definition that :math:`c_{s+1}^S=1`
+(and thus the fast IVP portion is evolved over the full time step,
+:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n}]`).
 
-These methods evolve the fast IVPs, :eq:`MRI_fast_IVP` and
-:eq:`MRI_embedding_fast_IVP`, on overlapping time intervals
-:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n-1}+c_i h^S]` and
-:math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n-1}, t_{n}]`, using initial
-conditions :math:`v_{0,i}=y_{n-1}` and :math:`\tilde{v}_0=y_{n-1}`.
 
 
 MERK Methods
 ------------
 
-MERK family of methods never include the stage updates, :eq:`MRI_implicit_solve`
-or :eq:`MRI_embedding_implicit_solve`.
+The MERK family of methods are only defined for multirate applications that
+are explicit at the slow time scale, i.e., :math:`f^I=0`, but otherwise they are
+nearly identical to IMEX-MRI-SR methods.  Specifically, like IMEX-MRI-SR methods,
+these evolve the fast IVPs
+:eq:`MRI_fast_IVP` and :eq:`MRI_embedding_fast_IVP` over the intervals
+:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n,i}^S]` and
+:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n}]`, respectively, and begin
+with the initial condition :math:`v_{0,i}=y_{n-1}`.  Furthermore, the fast IVP
+forcing functions are given by :eq:`IMEXMRISR_forcing` with :math:`f^I=0`.
+As MERK-based applications lack the implicit slow operator, they do not require
+the solution of implicit algebraic equations.
+
+However, unlike other MRI families, MERK methods were designed to admit a useful
+efficiency improvement.  Since each fast IVP begins with the same initial condition,
+:math:`v_{0,i}=y_{n-1}`, if multiple stages share the same forcing function
+:math:`r_i(t)`, then they may be "grouped" together.  This is achieved by sorting the
+final IVP solution time for each stage, :math:`t_{n,i}^S`, and then evolving the inner
+solver to each of these stage times in order, storing the corresponding inner solver
+solutions at these times as the stages :math:`z_i`.  For example, the
+:index:`ARKODE_MERK54` method involves 11 stages, that may be combined into 5 distinct
+groups.  The fourth group contains stages 7, 8, 9, and the embedding, corresponding to
+the :math:`c_i^S` values :math:`7/10`, :math:`1/2`, :math:`2/3`, and :math:`1`.
+Sorting these, a single fast IVP for this group must be evolved over the interval
+:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n}]`, first pausing at :math:`t_{n-1}+\frac12 h^S`
+to store :math:`z_8`, then pausing at :math:`t_{n-1}+\frac{2}{3} h^S` to store
+:math:`z_9`, then pausing at :math:`t_{n-1}+\frac{7}{10} h^S` to store :math:`z_7`,
+and finally finishing the IVP solve to :math:`t_{n-1}+h^S` to obtain :math:`\tilde{y}_n`.
+
+.. note::
+
+   Although all MERK methods were derived in :cite:p:`Luan:20` under an assumption that
+   the fast function :math:`f^F(t,y)` is linear in :math:`y`, in :cite:p:`Fish:24` it
+   was proven that MERK methods also satisfy all nonlinear order conditions up through
+   their linear order.  The lone exception is :index:`ARKODE_MERK54`, where it was only
+   proven to satisfy all nonlinear conditions up to order 4 (since :cite:p:`Fish:24` did
+   not establish the formulas for the order 5 conditions).  All our numerical tests to
+   date have shown :index:`ARKODE_MERK54` to achieve fifth order for nonlinear problems,
+   and so we conjecture that it also satisfies the nonlinear fifth order conditions.
 
-These methods evolve the fast IVPs, :eq:`MRI_fast_IVP` and
-:eq:`MRI_embedding_fast_IVP`, on overlapping time intervals
-:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n-1}+c_i h^S]` and
-:math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n-1}, t_{n}]`, using initial
-conditions :math:`v_{0,i}=y_{n-1}` and :math:`\tilde{v}_0=y_{n-1}`.
 
 
 .. _ARKODE.Mathematics.Error.Norm:
@@ -1063,17 +1126,25 @@ that it may have used a potentially different relative solution tolerance,
 estimates using either a "maximum accumulation" strategy,
 
 .. math::
-   \varepsilon^F_{max} = \text{RTOL}^F \max_{m\in M} \|\varepsilon^F_{n,m}\|_{WRMS},
+   \varepsilon^F_{max} = \text{RTOL}^F \max_{m\in \mathcal{S}} \|\varepsilon^F_{n,m}\|_{WRMS},
+   :label: maximum_accumulation
+
+an "additive accumulation" strategy,
+
+.. math::
+   \varepsilon^F_{sum} = \text{RTOL}^F \sum_{m\in \mathcal{S}} \|\varepsilon^F_{n,m}\|_{WRMS},
    :label: maximum_accumulation
 
 or using an "averaged accumulation" strategy,
 
 .. math::
-   \varepsilon^F_{avg} = \frac{\text{RTOL}^F}{|M|} \sum_{m\in M} \|\varepsilon^F_{n,m}\|_{WRMS},
+   \varepsilon^F_{avg} = \frac{\text{RTOL}^F}{\Delta t_{\mathcal{S}}} \sum_{m\in \mathcal{S}} h_{n,m} \|\varepsilon^F_{n,m}\|_{WRMS},
    :label: average_accumulation
 
-where the norm is taken using the tolerance-informed error-weight vector.  In
-both cases, the sum or the maximum is taken over the set of all steps :math:`M`
+where :math:`h_{n,m}` is the step size that gave rise to :math:`\varepsilon^F_{n,m}`,
+:math:`\Delta t_{\mathcal{S}}` denotes the elapsed time over which :math:`\mathcal{S}`
+is taken, and the norms are taken using the tolerance-informed error-weight vector.  In
+each case, the sum or the maximum is taken over the set of all steps :math:`\mathcal{S}`
 since the fast error accumulator has been reset.
 
 
@@ -1152,7 +1223,7 @@ initial state,
 where we again satisfy the inequality with a value of :math:`\frac12` to obtain
 
 .. math::
-   |h_0| = \frac{1}{2\left\| f(t_0,y_0)\right\|},
+   |h_0| = \frac{1}{2\left\| f(t_0,y_0)\right\|}.
    :label: H0_TSExp0
 
 
@@ -1565,7 +1636,7 @@ but may be modified by the user.
 
 When using the dense and band SUNMatrix objects for the linear systems
 :eq:`ARKODE_modified_Newton_system`, the Jacobian :math:`J` may be supplied
-by a user routine, or approximated internally by finite-differences.
+by a user routine, or approximated internally with finite-differences.
 In the case of differencing, we use the standard approximation
 
 .. math::

doc/arkode/guide/source/Usage/MRIStep/Custom_Inner_Stepper/Description.rst

@@ -237,7 +237,7 @@ Setting Member Functions
 
    :retval ARK_SUCCESS: if successful
    :retval ARK_ILL_INPUT: if the stepper is ``NULL``
-   
+
    .. versionadded: x.y.z
 
 
@@ -251,7 +251,7 @@ Setting Member Functions
 
    :retval ARK_SUCCESS: if successful
    :retval ARK_ILL_INPUT: if the stepper is ``NULL``
-      
+
    .. versionadded: x.y.z
 
 
@@ -265,7 +265,7 @@ Setting Member Functions
 
    :retval ARK_SUCCESS: if successful
    :retval ARK_ILL_INPUT: if the stepper is ``NULL``
-      
+
    .. versionadded: x.y.z
 
 
@@ -416,6 +416,7 @@ following member functions:
 
    **Example codes:**
       * ``examples/arkode/CXX_parallel/ark_diffusion_reaction_p.cpp``
+
    .. versionchanged:: v5.7.0
 
       Supplying a full right-hand side function was made optional.
@@ -424,7 +425,7 @@ following member functions:
 
    This function resets the inner (fast) stepper state to the provided
    independent variable value and dependent variable vector.
-   
+
    If provided, the :c:type:`MRIStepInnerResetFn` function will be called
    *before* a call to :c:type:`MRIStepInnerEvolveFn` when the state was
    updated at the slow timescale.
@@ -439,16 +440,17 @@ following member functions:
       value if a recoverable error occurred, or a negative value if it failed
       unrecoverably.
 
-   .. note::
-
-
    **Example codes:**
       * ``examples/arkode/CXX_parallel/ark_diffusion_reaction_p.cpp``
 
 
 .. c:type:: int (*MRIStepInnerGetAccumulatedError)(MRIStepInnerStepper stepper, sunrealtype* accum_error)
 
-   This function returns an estimate of the accumulated solution error arising from the inner stepper.
+   This function returns an estimate of the accumulated solution error arising from the
+   inner stepper.   Both the :c:type:`MRIStepInnerGetAccumulatedError` and
+   :c:type:`MRIStepInnerResetAccumulatedError` functions should be provided, or not; if only
+   one is provided then MRIStep will disable multirate temporal adaptivity and call neither.
+
 
    **Arguments:**
       * *stepper* -- the inner stepper object.
@@ -461,20 +463,19 @@ following member functions:
    .. note::
 
       This function is only called when multirate temporal adaptivity has been enabled,
-      using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL`.
+      using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL`.
 
       If provided, the :c:type:`MRIStepInnerGetAccumulatedError` function will always
       be called *after* a preceding call to the :c:type:`MRIStepInnerResetAccumulatedError`
       function.
 
-      Both the :c:type:`MRIStepInnerGetAccumulatedError` and
-      :c:type:`MRIStepInnerResetAccumulatedError` functions should be provided, or not; if only
-      one is provided then MRIStep will disable multirate temporal adaptivity and call neither.
-
 
 .. c:type:: int (*MRIStepInnerResetAccumulatedError)(MRIStepInnerStepper stepper)
 
    This function resets the inner stepper's accumulated solution error to zero.
+   This function performs a different role within MRIStep than the
+   :c:type:`MRIStepInnerResetFn`, and thus an implementation should make no
+   assumptions about the frequency/ordering of calls to either.
 
    **Arguments:**
       * *stepper* -- the inner stepper object.
@@ -486,7 +487,7 @@ following member functions:
    .. note::
 
       This function is only called when multirate temporal adaptivity has been enabled,
-      using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL`.
+      using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL`.
 
       The :c:type:`MRIStepInnerResetAccumulatedError` function will always be called
       *before* any calls to the :c:type:`MRIStepInnerGetAccumulatedError` function.
@@ -495,14 +496,13 @@ following member functions:
       :c:type:`MRIStepInnerResetAccumulatedError` functions should be provided, or not; if only
       one is provided then MRIStep will disable multirate temporal adaptivity and call neither.
 
-      This function performs a different role within MRIStep than the
-      :c:type:`MRIStepInnerResetFn`, and thus an implementation should make no
-      assumptions about the frequency/ordering of calls to either.
-
 
 .. c:type:: int (*MRIStepInnerSetRTol)(MRIStepInnerStepper stepper, sunrealtype rtol)
 
-   This function accepts a relative tolerance for the inner stepper to use in its upcoming adaptive solve.
+   This function accepts a relative tolerance for the inner stepper to use in its
+   upcoming adaptive solve.  It is assumed that if the inner stepper supports absolute
+   tolerances as well, then these have been set up directly by the user to indicate the
+   "noise" level for solution components.
 
    **Arguments:**
       * *stepper* -- the inner stepper object.
@@ -515,8 +515,4 @@ following member functions:
    .. note::
 
       This function is only called when multirate temporal adaptivity has been enabled
-      using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL`.
-
-      It is assumed that if the inner stepper supports absolute tolerances as well, then
-      these have been set up directly by the user to indicate the "noise" level for
-      solution components.
+      using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL`.

doc/arkode/guide/source/Usage/MRIStep/MRIStepCoupling.rst

@@ -273,8 +273,6 @@ are defined ``arkode/arkode_mristep.h``.
               A ``NULL`` pointer if ``stages`` was invalid, an allocation error occurred,
               or the input data arrays are inconsistent with the method type.
 
-   .. note::
-
 
 .. c:function:: MRIStepCoupling MRIStepCoupling_MIStoMRI(ARKodeButcherTable B, int q, int p)
 
@@ -425,19 +423,23 @@ with values specified for each method below (e.g., ``ARKODE_MIS_KW3``).
 
 Notes regarding the above table:
 
-#. The default method for each order when using fixed step sizes is marked with an asterisk (:math:`^*`).
+#. The default method for each order when using fixed step sizes is marked with an
+   asterisk (:math:`^*`).
 
-#. The default method for each order when using adaptive time stepping is marked with a circle (:math:`^\circ`).
+#. The default method for each order when using adaptive time stepping is marked
+   with a circle (:math:`^\circ`).
 
-#. The "Slow RHS Calls" column corresponds to the number of calls to the slow right-hand side function,
-   :math:`f^E`, per time step.
+#. The "Slow RHS Calls" column corresponds to the number of calls to the slow
+   right-hand side function, :math:`f^E`, per time step.
 
-#. Note A: although all MERK methods were derived in :cite:p:`Luan:20` under an assumption that the fast time
-   scale is linear in the solution, in :cite:p:`Fish:24` it was proven that they also satisfy all nonlinear
-   order conditions up through their linear order.  The lone exception is MERK54, where it was only proven to
-   satisfy all nonlinear conditions up to order 4, since :cite:p:`Fish:24` did not establish the formulas for
-   the order 5 conditions.  All our numerical tests to date have shown MERK54 to achieve fifth order for
-   nonlinear problems, and so we conjecture that it also satisfies the nonlinear fifth order conditions.
+#. Note A: although all MERK methods were derived in :cite:p:`Luan:20` under an assumption
+   that the fast function :math:`f^F(t,y)` is linear in :math:`y`, in :cite:p:`Fish:24` it
+   was proven that MERK methods also satisfy all nonlinear order conditions up through
+   their linear order.  The lone exception is :index:`ARKODE_MERK54`, where it was only
+   proven to satisfy all nonlinear conditions up to order 4 (since :cite:p:`Fish:24` did
+   not establish the formulas for the order 5 conditions).  All our numerical tests to
+   date have shown :index:`ARKODE_MERK54` to achieve fifth order for nonlinear problems,
+   and so we conjecture that it also satisfies the nonlinear fifth order conditions.
 
 
 .. table:: Diagonally-implicit, solve-decoupled MRI-GARK coupling tables. The default

doc/arkode/guide/source/Usage/MRIStep/Skeleton.rst

@@ -85,11 +85,16 @@ unchanged from the skeleton program presented in
    inputs to :c:func:`MRIStepCreate` is the ``MRIStepInnerStepper`` object for
    solving the fast (inner) IVP created in the previous step.
 
-#. Set the slow step size
-
-   Call :c:func:`ARKodeSetFixedStep()` on the MRIStep object to specify the
+#. If using fixed step sizes, then set the slow step size by calling
+   :c:func:`ARKodeSetFixedStep()` on the MRIStep object to specify the
    slow time step size.
 
+   If using adaptive slow steps, then specify the desired integration tolerances
+   as normal.  By default, MRIStep will use a "decoupled" (see
+   :numref:`ARKODE.Mathematics.MultirateControllers`) I controller (see
+   :numref:`SUNAdaptController.Soderlind`),  Alternately, create and attach a
+   multirate temporal controller (see :numref:`SUNAdaptController.MRIHTol`).
+
 #. Create and configure implicit solvers (*as appropriate*)
 
    Specifically, if MRIStep is configured with an implicit slow right-hand side