Relaxation Runge-Kutta Methods (#298)

Add support for explicit, implicit, and IMEX relaxation Runge-Kutta methods in ERKStep and ARKStep --------- Co-authored-by: Daniel R. Reynolds <[email protected]> Co-authored-by: Cody Balos <[email protected]>
LLNL · Jul 17, 2023 · cd0d062 · cd0d062
1 parent df5f4b9
commit cd0d062
Show file tree

Hide file tree

Showing 61 changed files with 6,538 additions and 60 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,9 @@
 
 ## Changes to SUNDIALS in release 6.6.0
 
+Added support for relaxation Runge-Kutta methods to ERKStep and ARKStep in
+ARKODE.
+
 Added the second order IMEX method from Giraldo, Kelly, and Constantinescu 2013
 as the default second order IMEX method in ARKStep. The explicit table is given
 by `ARKODE_ARK2_ERK_3_1_2` and the implicit table by `ARKODE_ARK2_DIRK_3_1_2`.

diff --git a/doc/arkode/guide/source/Constants.rst b/doc/arkode/guide/source/Constants.rst
@@ -62,6 +62,16 @@ contains the ARKODE output constants.
    +---------------------------------------------+-----------------------------------------------------------+
    |                                                                                                         |
    +---------------------------------------------+-----------------------------------------------------------+
+   | **Relaxtion module input constants**                                                                    |
+   +---------------------------------------------+-----------------------------------------------------------+
+   | :index:`ARK_RELAX_BRENT`                    | Specifies Brent's method as the relaxation nonlinear      |
+   |                                             | solver.                                                   |
+   +---------------------------------------------+-----------------------------------------------------------+
+   | :index:`ARK_RELAX_NEWTON`                   | Specifies Newton's method as the relaxation nonlinear     |
+   |                                             | solver.                                                   |
+   +---------------------------------------------+-----------------------------------------------------------+
+   |                                                                                                         |
+   +---------------------------------------------+-----------------------------------------------------------+
    | **Explicit Butcher table specification**                                                                |
    +---------------------------------------------+-----------------------------------------------------------+
    | :index:`ARKODE_HEUN_EULER_2_1_2`            | Use the Heun-Euler-2-1-2 ERK method.                      |
@@ -377,6 +387,17 @@ contains the ARKODE output constants.
    +-------------------------------------+------+------------------------------------------------------------+
    | :index:`ARK_INVALID_TABLE`          | -41  | An invalid Butcher or MRI table was encountered.           |
    +-------------------------------------+------+------------------------------------------------------------+
+   | :index:`ARK_CONTEXT_ERR`            | -42  | An error occurred with the SUNDIALS context object         |
+   +-------------------------------------+------+------------------------------------------------------------+
+   | :index:`ARK_RELAX_FAIL`             | -43  | An error occurred in computing the relaxation parameter    |
+   +-------------------------------------+------+------------------------------------------------------------+
+   | :index:`ARK_RELAX_MEM_FAIL`         | -44  | The relaxation memory structure is ``NULL``                |
+   +-------------------------------------+------+------------------------------------------------------------+
+   | :index:`ARK_RELAX_FUNC_FAIL`        | -45  | The relaxation function returned an unrecoverable error    |
+   +-------------------------------------+------+------------------------------------------------------------+
+   | :index:`ARK_RELAX_JAC_FAIL`         | -46  | The relaxation Jacobian function returned an unrecoverable |
+   |                                     |      | error                                                      |
+   +-------------------------------------+------+------------------------------------------------------------+
    | :index:`ARK_UNRECOGNIZED_ERROR`     | -99  | An unknown error was encountered.                          |
    +-------------------------------------+------+------------------------------------------------------------+
    |                                                                                                         |

diff --git a/doc/arkode/guide/source/Introduction.rst b/doc/arkode/guide/source/Introduction.rst
@@ -121,6 +121,11 @@ Changes from previous versions
 Changes in v5.6.0
 -----------------
 
+
+Added support for relaxation Runge-Kutta methods in ERKStep and ARKStep, see
+:numref:`ARKODE.Mathematics.Relaxation`, :numref:`ARKODE.Usage.ERKStep.Relaxation`,
+and :numref:`ARKODE.Usage.ARKStep.Relaxation` for more information.
+
 Added the second order IMEX method from :cite:p:`giraldo2013implicit` as the
 default second order IMEX method in ARKStep. The explicit table is given by
 ``ARKODE_ARK2_ERK_3_1_2`` (see :numref:`Butcher.ARK2_ERK`) and the implicit

diff --git a/doc/arkode/guide/source/Mathematics.rst b/doc/arkode/guide/source/Mathematics.rst
@@ -2074,3 +2074,68 @@ step attempt, or fails with the minimum step size, then the integration is halte
 and an error is returned. In this case the user may need to employ other
 strategies as discussed in :numref:`ARKODE.Usage.ARKStep.Tolerances` and
 :numref:`ARKODE.Usage.ERKStep.Tolerances` to satisfy the inequality constraints.
+
+.. _ARKODE.Mathematics.Relaxation:
+
+Relaxation Methods
+==================
+
+When the solution of :eq:`ARKODE_IVP` is conservative or dissipative with
+respect to a smooth *convex* function :math:`\xi(y(t))`, it is desirable to have
+the numerical method preserve these properties. That is
+:math:`\xi(y_n) = \xi(y_{n-1}) = \ldots = \xi(y_{0})` for conservative systems
+and :math:`\xi(y_n) \leq \xi(y_{n-1})` for dissipative systems. For examples
+of such problems, see the references below and the citations there in.
+
+For such problems, ARKODE supports relaxation methods
+:cite:p:`ketcheson2019relaxation, kang2022entropy, ranocha2020relaxation, ranocha2020hamiltonian`
+applied to ERK, DIRK, or ARK methods to ensure dissipation or preservation of
+the global function. The relaxed solution is given by
+
+.. math::
+   y_r = y_{n-1} + r d = r y_n + (1 - r) y_{n - 1}
+   :label: ARKODE_RELAX_SOL
+
+where :math:`d` is the update to :math:`y_n` (i.e.,
+:math:`h_n \sum_{i=1}^{s}(b^E_i \hat{f}^E_i + b^I_i \hat{f}^I_i)` for ARKStep
+and :math:`h_n \sum_{i=1}^{s} b_i f_i` for ERKStep) and :math:`r` is the
+relaxation factor selected to ensure conservation or dissipation. Given an ERK,
+DIRK, or ARK method of at least second order with non-negative solution weights
+(i.e., :math:`b_i \geq 0` for ERKStep or :math:`b^E_i \geq 0` and
+:math:`b^I_i \geq 0` for ARKStep), the factor :math:`r` is computed by solving
+the auxiliary scalar nonlinear system
+
+.. math::
+   F(r) = \xi(y_{n-1} + r d) - \xi(y_{n-1}) - r e = 0
+   :label: ARKODE_RELAX_NLS
+
+at the end of each time step. The estimated change in :math:`\xi` is given by
+:math:`e \equiv h_n \sum_{i=1}^{s} \langle \xi'(z_i), b^E_i f^E_i + b^I_i f^I_i \rangle`
+where :math:`\xi'` is the Jacobian of :math:`\xi`.
+
+Two iterative methods are provided for solving :eq:`ARKODE_RELAX_NLS`, Newton's
+method and Brent's method. When using Newton's method (the default), the
+iteration is halted either when the residual tolerance is met,
+:math:`F(r^{(k)}) < \epsilon_{\mathrm{relax\_res}}`, or when the difference
+between successive iterates satisfies the relative and absolute tolerances,
+:math:`|\delta_r^{(k)}| = |r^{(k)} - r^{(k-1)}| < \epsilon_{\mathrm{relax\_rtol}} |r^{(k-1)}| + \epsilon_{\mathrm{relax\_atol}}`.
+Brent's method applies the same residual tolerance check and additionally halts
+when the bisection update satisfies the relative and absolute tolerances,
+:math:`|0.5 (r_c - r^{k})| < \epsilon_{\mathrm{relax\_rtol}} |r^{(k)}| + 0.5 \epsilon_{\mathrm{relax\_atol}}`
+where :math:`r_c` and :math:`r^{(k)}` bound the root.
+
+If the nonlinear solve fails to meet the specified tolerances within the maximum
+allowed number of iterations, the step size is reduced by the factor
+:math:`\eta_\mathrm{rf}` (default 0.25) and the step is repeated. Additionally,
+the solution of :eq:`ARKODE_RELAX_NLS` should be
+:math:`r = 1 + \mathcal{O}(h_n^{q - 1})` for a method of order :math:`q`
+:cite:p:`ranocha2020relaxation`. As such, limits are imposed on the range of
+relaxation values allowed (i.e., limiting the maximum change in step size due to
+relaxation). A relaxation value greater than :math:`r_\text{max}` (default 1.2)
+or less than :math:`r_\text{min}` (default 0.8), is considered as a failed
+relaxation application and the step will is repeated with the step size reduced
+by :math:`\eta_\text{rf}`.
+
+For more information on utilizing relaxation Runge--Kutta methods, see
+:numref:`ARKODE.Usage.ERKStep.Relaxation` and
+:numref:`ARKODE.Usage.ARKStep.Relaxation`.