Added explicit SSP Runge-Kutta time integrators

Closes idaholab#11855

framework/doc/content/bib/moose.bib

@@ -253,3 +253,30 @@ @article{macneal1985patch
   issue = {1},
   pages = {3--20}
 }
+
+@article{gottlieb2005,
+  title = {On High Order Strong Stability Preserving Runge–Kutta and Multi Step Time Discretizations},
+  author = {Sigal Gottlieb},
+  journal = {Journal of Scientific Computing},
+  volume = {25},
+  number = {1/2},
+  month = {November},
+  year = {2005},
+  doi = {10.1007/s10915-004-4635-5}
+}
+
+@article{zhao2019,
+       author = {{Zhao}, Weifeng and {Huang}, Juntao},
+        title = "{Boundary treatment of implicit-explicit Runge-Kutta method for hyperbolic systems with source terms}",
+      journal = {arXiv e-prints},
+     keywords = {Mathematics - Numerical Analysis},
+         year = "2019",
+        month = "Aug",
+          eid = {arXiv:1908.01027},
+        pages = {arXiv:1908.01027},
+archivePrefix = {arXiv},
+       eprint = {1908.01027},
+ primaryClass = {math.NA},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2019arXiv190801027Z},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

framework/doc/content/source/timeintegrators/ExplicitSSPRungeKutta.md

@@ -0,0 +1,222 @@
+# ExplicitSSPRungeKutta
+
+## Introduction
+
+This is a base class for Strong Stability Preserving (SSP) Runge-Kutta time
+integration methods, which are fully explicit methods of varying orders of
+accuracy and number of stages. The key feature of these methods is that they
+preserve the strong stability properties (in any norm or seminorm) of the
+explicit/forward Euler method [!cite](gottlieb2005).
+
+These methods derive from [/ExplicitTimeIntegrator.md].
+
+## Formulation
+
+For the ODE
+\begin{equation}
+  \frac{d u}{d t} = f(u(t), t) \,,
+\end{equation}
+the SSP Runge-Kutta methods up to order 3 can be expressed in the following form
+for a time step:
+\begin{equation}
+  u^{(0)} = u^n
+\end{equation}
+\begin{equation}
+  u^{(s)} = \sum\limits_{k=0}^{s-1} a_{s,k} u^{(k)} + b_s \Delta t f\left(u^{(s-1)}, t^n + c_s\Delta t\right) \,,
+    \qquad i = 1,\ldots,N_s
+\end{equation}
+\begin{equation}
+  u^{n+1} = u^{(N_s)}
+\end{equation}
+where $N_s$ is the number of stages and for methods up to order 3, is also the
+order of accuracy. The coefficients $a_{i,j}$, $b_i$, and $c_i$ can be
+conveniently expressed in the following tabular form:
+\begin{equation}
+  \begin{array}{c|c}
+    \mathbf{c} & \mathbf{A}\\\hline
+    & \mathbf{b}^T\\
+  \end{array}
+\end{equation}
+These methods have the following time step size requirement for stability:
+\begin{equation}
+  \Delta t < \gamma \Delta t_{FE} \,, \quad \gamma = \min\limits_{i,k} \frac{a_{i,k}}{b_i}
+\end{equation}
+where $\Delta t_{FE}$ is the maximum time step size for stability of the forward
+Euler method.
+
+In MOOSE, generally the system of ODEs to be solved result from discretization
+using the finite element method, and thus a mass matrix exists:
+\begin{equation}
+  \mathbf{M}\frac{d \mathbf{u}}{d t} = \mathbf{f}\left(\mathbf{u}(t), t\right)
+\end{equation}
+In this case, the stage $s$ solution is actually the following:
+\begin{equation}
+  \mathbf{u}^{(s)} = \sum\limits_{k=0}^{s-1} a_{s,k} \mathbf{u}^{(k)}
+    + \mathbf{M}^{-1} b_s \Delta t \mathbf{f}\left(\mathbf{u}^{(s-1)}, t^n + c_s\Delta t\right)
+\end{equation}
+As an implementation note, the usual mass matrix entry is
+\begin{equation}
+  m_{i,j} = (\phi_i, \phi_j)
+\end{equation}
+However, in MOOSE, the mass matrix includes the time step size:
+\begin{equation}
+  m_{i,j} = \frac{(\phi_i, \phi_j)}{b_i \Delta t}
+\end{equation}
+\begin{equation}\label{eq:ssprk_no_bc}
+  \mathbf{u}^{(s)} = \sum\limits_{k=0}^{s-1} a_{s,k} \mathbf{u}^{(k)}
+    + \mathbf{M}^{-1} \mathbf{f}\left(\mathbf{u}^{(s-1)}, t^n + c_s\Delta t\right)
+\end{equation}
+
+## Dirichlet Boundary Conditions Treatment
+
+Now consider the case where one or more degrees of freedom are subject to strong
+(Dirichlet) boundary conditions:
+\begin{equation}
+  u_i(t) = g(\mathbf{x}_i, t) \,, \qquad i \in \mathcal{I}_\text{Dirichlet}
+\end{equation}
+For a nonlinear solve with Newton's method, each iteration consists of the
+solution of a linear system:
+\begin{equation}\label{eq:ssprk_newton_update}
+  \mathbf{J}^{(\ell-1)}\mathbf{\delta u}^{(\ell)} = -\mathbf{r}^{(\ell-1)}
+\end{equation}
+and then updating the solution:
+\begin{equation}
+  \mathbf{u}^{(\ell)} = \mathbf{u}^{(\ell-1)} + \mathbf{\delta u}^{(\ell)} \,.
+\end{equation}
+In MOOSE, Dirichlet boundary conditions are implemented by modifying the residual
+vector $\mathbf{r}$ to replace entries for the affected degrees of freedom:
+\begin{equation}
+  \tilde{r}_i = \left\{\begin{array}{l l}
+    r_i \,,                      & i \notin \mathcal{I}_\text{Dirichlet}\\
+    u_i - g(\mathbf{x}_i, t) \,, & i \in \mathcal{I}_\text{Dirichlet}
+    \end{array}\right.
+\end{equation}
+By modifying the Jacobian matrix as follows, one can guarantee that the
+boundary conditions are enforced, i.e., $u_i^{(\ell)} = g(\mathbf{x}_i, t)$ for
+$i \in \mathcal{I}_\text{Dirichlet}$:
+\begin{equation}
+  \tilde{J}_{i,j} = \left\{\begin{array}{l l}
+    J_{i,j} \,, & i \notin \mathcal{I}_\text{Dirichlet}\\
+    1 \,,       & i \in \mathcal{I}_\text{Dirichlet}, \quad j = i\\
+    0 \,,       & i \in \mathcal{I}_\text{Dirichlet}, \quad j \ne i
+    \end{array}\right.
+\end{equation}
+To work with MOOSE's Dirichlet boundary condition implementation, [eq:ssprk_no_bc]
+must be put in an update form, similar to [eq:ssprk_newton_update]:
+\begin{equation}\label{eq:ssprk_stage_update}
+  \mathbf{M}\mathbf{\delta u}^{(s)} = \mathbf{b}
+\end{equation}
+\begin{equation}\label{eq:ssprk_b}
+  \mathbf{b} = \mathbf{M}\left(\sum\limits_{k=0}^{s-1} a_{s,k} \mathbf{u}^{(k)}
+    - \mathbf{u}^{(s-1)} \right)
+    + \mathbf{f}\left(\mathbf{u}^{(s-1)}, t^n + c_s\Delta t\right)
+\end{equation}
+\begin{equation}
+  \mathbf{\delta u}^{(s)} \equiv \mathbf{u}^{(s)} - \mathbf{u}^{(s-1)}
+\end{equation}
+To impose the Dirichlet boundary conditions, the mass matrix and right-hand
+side vector are modified as for the Newton case:
+\begin{equation}
+  \tilde{m}_{i,j} = \left\{\begin{array}{l l}
+    m_{i,j} \,, & i \notin \mathcal{I}_\text{Dirichlet}\\
+    1 \,,       & i \in \mathcal{I}_\text{Dirichlet}, \quad j = i\\
+    0 \,,       & i \in \mathcal{I}_\text{Dirichlet}, \quad j \ne i
+    \end{array}\right.
+\end{equation}
+\begin{equation}\label{eq:ssprk_b_mod}
+  \tilde{b}_i = \left\{\begin{array}{l l}
+    b_i \,,                              & i \notin \mathcal{I}_\text{Dirichlet}\\
+    g_i^{(s)} - u_i^{(s-1)} \,, & i \in \mathcal{I}_\text{Dirichlet}
+    \end{array}\right.
+\end{equation}
+where $g_i^{(s)}$ is an appropriate value to impose for degree of freedom $i$
+in stage $s$. For most cases, this is simply
+\begin{equation}
+  g_i^{(s)} = g(\mathbf{x}_i, t^{n+1})
+\end{equation}
+However, in general, certain conditions must be enforced on the imposed boundary values
+for intermediate stages to preserve the formal order of accuracy of the method.
+For methods up to order 2, it is safe to impose each stage as shown above.
+For methods greater than order 2, the conditions to preserve order of accuracy
+will be presented on a case-by-case basis.
+
+## Usage
+
+Classes that derive from this class are responsible for the following:
+
+- overriding the `order()` method
+- setting the Runge-Kutta coefficient data `_a`, `_b`, `_c`, and `_s` in the
+  constructor
+- resizing `_solution_stage` to size `_s + 1` in the constructor
+- overriding the `solve()` method (see below)
+
+### Implementing the `solve()` method
+
+Firstly, the counts of nonlinear and linear iterations should be reset to zero:
+
+```
+_n_nonlinear_iterations = 0;
+_n_linear_iterations = 0;
+```
+
+Then, each stage needs to be solved. For each stage, the first step is to set
+the stage index (first is 0):
+
+```
+_stage = 0;
+```
+
+Then, for stage $i$, pointers to solutions of stages up to $i-1$ need to be
+stored. For example, for the first stage of any method,
+
+```
+_solution_stage[0] = &_solution_old;
+```
+
+This becomes trickier with increasing number of stages because storage for
+intermediate stages may need to be created. The next step is to set the time
+to the stage time and then update the nonlinear implicit system:
+
+```
+_fe_problem.time() = time_old + _c[_stage] * dt;
+_nonlinear_implicit_system->update();
+```
+
+Finally, the stage solve is performed with `solveStage()`, and a boolean is
+returned that is false if the solved failed. If more stages are to be performed,
+and the stage failed, then one should just return from `solve()` prematurely.
+
+## Implementation
+
+[eq:ssprk_stage_update] is implemented as described in the following sections:
+
+### `computeTimeDerivatives()`
+
+Only the Jacobian `_du_dot_du` is implemented, which is needed by the mass matrix.
+The time derivative itself is not needed because only part of it appears in the
+residual vector.
+
+### `solveStage()`
+
+First the mass matrix is computed by calling `computeJacobianTag()` with the
+time tag. Because the mass matrix is computed before the call
+`computeResidual()`, the call to `computeTimeDerivatives()` must be made before
+`computeJacobianTag()`, even though it will be called again in
+`computeResidual()`. The Jacobian must be computed before the call to
+`computeResidual()` because the mass matrix will be used in `computeResidual()`
+via the call to `postResidual()`. In `computeResidual()`, the following steps
+occur:
+
+- The steady-state residual $\mathbf{f}\left(\mathbf{u}^{(s-1)}, t^n + c_s\Delta t\right)$
+  is computed.
+- `postResidual()` is called, which adds the mass matrix product shown in [eq:ssprk_b].
+- Dirichlet boundary condition modifications are made to the residual as shown
+  in [eq:ssprk_b_mod].
+
+### `postResidual()`
+
+Here $\mathbf{b}$ is assembled as shown in [eq:ssprk_b]. The mass matrix product
+here is responsible for the need to call `computeJacobianTag()` before `computeResidual()`
+in `solveStage()`.
+
+!bibtex bibliography

framework/doc/content/source/timeintegrators/ExplicitSSPRungeKutta1.md

@@ -0,0 +1,20 @@
+# ExplicitSSPRungeKutta1
+
+This time integrator is the first order method of the SSP Runge-Kutta family,
+which has the base class [/ExplicitSSPRungeKutta.md]. It has the following tableaux:
+
+\begin{equation}
+  \begin{array}{c|c}
+    0 & 1\\\hline
+    & 1\\
+  \end{array}
+\end{equation}
+
+The first order method is simply the explicit/forward Euler method; thus it
+should be equivalent to [/ActuallyExplicitEuler.md].
+
+!syntax parameters /Executioner/TimeIntegrator/ExplicitSSPRungeKutta1
+
+!syntax inputs /Executioner/TimeIntegrator/ExplicitSSPRungeKutta1
+
+!syntax children /Executioner/TimeIntegrator/ExplicitSSPRungeKutta1

framework/doc/content/source/timeintegrators/ExplicitSSPRungeKutta2.md

@@ -0,0 +1,18 @@
+# ExplicitSSPRungeKutta2
+
+This time integrator is the second order method of the SSP Runge-Kutta family,
+which has the base class [/ExplicitSSPRungeKutta.md]. It has the following tableaux:
+
+\begin{equation}
+  \begin{array}{c|c c}
+    0 & 1           & \\
+    1 & \frac{1}{2} & \frac{1}{2}\\\hline
+      & 1           & \frac{1}{2}\\
+  \end{array}
+\end{equation}
+
+!syntax parameters /Executioner/TimeIntegrator/ExplicitSSPRungeKutta2
+
+!syntax inputs /Executioner/TimeIntegrator/ExplicitSSPRungeKutta2
+
+!syntax children /Executioner/TimeIntegrator/ExplicitSSPRungeKutta2

framework/doc/content/source/timeintegrators/ExplicitSSPRungeKutta3.md

@@ -0,0 +1,36 @@
+# ExplicitSSPRungeKutta3
+
+This time integrator is the third order method of the SSP Runge-Kutta family,
+which has the base class [/ExplicitSSPRungeKutta.md]. It has the following tableaux:
+
+\begin{equation}
+  \begin{array}{c|c c}
+    0           & 1           &             &\\
+    1           & \frac{3}{4} & \frac{1}{4} &\\
+    \frac{1}{2} & \frac{1}{3} & 0           & \frac{2}{3}\\\hline
+                & 1           & \frac{1}{4} & \frac{2}{3}\\
+  \end{array}
+\end{equation}
+
+## Dirichlet Boundary Conditions Treatment
+
+As discussed [here](/ExplicitSSPRungeKutta.md), the values imposed for Dirichlet
+boundary conditions in intermediate stages must satisfy certain conditions
+to preserve the formal order of accuracy. For this method, the boundary
+values imposed in each stage should be as follows, according to [!cite](zhao2019):
+\begin{equation}
+  u_i^{(1)} = g(\mathbf{x}_i, t^n) + \Delta t g'(\mathbf{x}_i, t^n)
+\end{equation}
+\begin{equation}
+  u_i^{(2)} = g(\mathbf{x}_i, t^n) + \frac{1}{2}\Delta t g'(\mathbf{x}_i, t^n)
+    + \frac{1}{4} \Delta t^2 g''(\mathbf{x}_i, t^n)
+\end{equation}
+\begin{equation}
+  u_i^{(3)} = g(\mathbf{x}_i, t^{n+1})
+\end{equation}
+
+!syntax parameters /Executioner/TimeIntegrator/ExplicitSSPRungeKutta3
+
+!syntax inputs /Executioner/TimeIntegrator/ExplicitSSPRungeKutta3
+
+!syntax children /Executioner/TimeIntegrator/ExplicitSSPRungeKutta3

framework/include/timeintegrators/ExplicitSSPRungeKutta.h

@@ -0,0 +1,57 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "ExplicitTimeIntegrator.h"
+
+class ExplicitSSPRungeKutta;
+
+template <>
+InputParameters validParams<ExplicitSSPRungeKutta>();
+
+/**
+ * Explicit strong stability preserving Runge-Kutta methods
+ */
+class ExplicitSSPRungeKutta : public ExplicitTimeIntegrator
+{
+public:
+  ExplicitSSPRungeKutta(const InputParameters & parameters);
+
+  virtual void computeTimeDerivatives() override;
+  virtual void computeADTimeDerivatives(DualReal & ad_u_dot,
+                                        const dof_id_type & dof) const override;
+  virtual void postResidual(NumericVector<Number> & residual) override;
+
+protected:
+  /**
+   * Solves a stage of the time integrator
+   *
+   * @returns true if converged, false if not
+   */
+  bool solveStage();
+
+  /// Number of stages
+  unsigned int _s;
+  /// Runge-Kutta "a" coefficient matrix
+  std::vector<std::vector<Real>> _a;
+  /// Runge-Kutta "b" coefficient vector
+  std::vector<Real> _b;
+  /// Runge-Kutta "c" coefficient vector
+  std::vector<Real> _c;
+  /// Current stage
+  unsigned int _stage;
+  /// Pointer to solution vector for each stage
+  std::vector<const NumericVector<Number> *> _solution_stage;
+
+  /// Temporary solution vector
+  NumericVector<Number> & _tmp_solution;
+  /// Temporary mass-matrix/solution vector product
+  NumericVector<Number> & _tmp_mass_solution_product;
+};