ARKODE: A Flexible IVP Solver Infrastructure for One-step Methods

2.1 ERKStep

The ERKStep time-stepping module is designed for IVPs of the form (3) \(\begin{equation} y^{\prime }(t) = f(t,y), \qquad y(t_0) = y_0, \end{equation}\) i.e., \(M(t)=I\) in (1). For such problems, ERKStep provides variable-step, embedded, explicit Runge–Kutta(ERK) methods that evolve a single time step \(t_{n-1}\rightarrow t_{n-1}+h_n = t_n\) using the algorithm (4a) \(\begin{align} z_i &= y_{n-1} + h_n \sum _{j=1}^{i-1} A_{i,j} f(t_{n,j}, z_j), \quad i=1,\ldots ,s, \end{align}\) (4b) \(\begin{align} y_n &= y_{n-1} + h_n \sum _{i=1}^{s} b_i f(t_{n,i}, z_i), \end{align}\) (4c) \(\begin{align} \tilde{y}_n &= y_{n-1} + h_n \sum _{i=1}^{s} \tilde{b}_i f(t_{n,i}, z_i), \end{align}\)

where the internal stage times are abbreviated using the notation \(t_{n,j}=t_{n-1}+c_jh_n\), and the coefficients \(A\in \mathbb {R}^{s\times s}\), \(b\in \mathbb {R}^s\), \(c\in \mathbb {R}^s\), and optional embedding coefficients \(\tilde{b}\in \mathbb {R}^s\) (a.k.a., the Butcher table) define the method. The embedding, \(\tilde{y}_n\), typically provides a slightly lower accuracy approximation than the computed solution, \(y_n\). ARKODE provides a number of ERK Butcher tables having orders of accuracy \(q=\lbrace 2,3,4,5,6,8\rbrace\); each is supplied with an embedding \(\tilde{b}\) that usually has order of accuracy \(p=q-1\). User-supplied Butcher tables are also supported and are not required to include embedding coefficients, although if these are not provided then automatic temporal adaptivity will be disabled. We have found this “fixed-step” mode to be particularly useful when existing simulation codes transition to ARKODE. By providing the ERK Butcher table that a code currently uses and ceding control over \(h_n\) to the application code, users may run verification tests to ensure reproducibility before switching to higher-order and/or adaptive ERK methods.

In order to allow users to “tune” ERKStep to better exploit the problem structure, numerous optional features may be enabled or modified by the user. The full set of features is described in [48], a few of which are highlighted in the example problems in Section 5.1.1.

2.2 ARKStep

The ARKStep time-stepping module is designed for IVPs where the right-hand side function may be additively split into two components: (5) \(\begin{equation} M(t)\, y^{\prime }(t) = f^E(t,y) + f^I(t,y), \qquad y(t_0) = y_0, \end{equation}\) where the user-provided right-hand side function \(f^E(t,y)\) contains the “nonstiff” components of the system, and \(f^I(t,y)\) contains the “stiff” components of the system. In solving Equation (5), we first convert to standard additive form: (6) \(\begin{equation} y^{\prime }(t) = \hat{f}^E(t,y) + \hat{f}^I(t,y), \qquad y(t_0) = y_0, \end{equation}\) where \(\hat{f}^E(t,y) = M(t)^{-1}f^E(t,y)\) and \(\hat{f}^I(t,y) = M(t)^{-1}f^I(t,y)\). ARKStep utilizes variable-step, embedded, ImEx additive Runge–Kutta(ARK) methods. Given a pair of s-stage Butcher tables \((A^E, b^E, \tilde{b}^E, c^E)\) and \((A^I, b^I, \tilde{b}^I, c^I)\), a single time step \(t_{n-1}\rightarrow t_{n}\) is advanced via the algorithm (7a) \(\begin{align} z_i &= y_{n-1} + h_n \sum _{j=1}^{i-1} A^E_{i,j} \hat{f}^E\left(t^E_{n,j}, z_j\right) + h_n \sum _{j=1}^{i} A^I_{i,j} \hat{f}^I\left(t^I_{n,j}, z_j\right)\!, \quad i=1,\ldots ,s, \end{align}\) (7b) \(\begin{align} y_n &= y_{n-1} + h_n \sum _{i=1}^{s} \left(b^E_i \hat{f}^E\left(t^E_{n,i}, z_i\right) + b^I_i \hat{f}^I\left(t^I_{n,i}, z_i\right)\right)\!, \end{align}\) (7c) \(\begin{align} \tilde{y}_n &= y_{n-1} + h_n \sum _{i=1}^{s} \left(\tilde{b}^E_i \hat{f}^E\left(t^E_{n,i}, z_i\right) + \tilde{b}^I_i \hat{f}^I\left(t^I_{n,i}, z_i\right)\right)\!, \end{align}\)

where we denote the internal stage times as \(t^E_{n,j} = t_{n-1} + c^E_j h_n\) and \(t^I_{n,j} = t_{n-1} + c^I_j h_n\). We note that ARKStep enforces the constraint that the explicit and implicit Butcher tables share the same number of stages, s, but it allows for the possibility of different explicit and implicit abscissae, i.e., \(c^E \ne c^I\) (as required for various SSP ImEx-ARK methods [29, 30, 31, 42]).

Through appropriate selection of the Butcher tables and right-hand side functions, ARKStep supports three classes of Runge–Kutta methods: ImEx, explicit, and diagonally implicit. Several Butcher tables are provided for each class, all of which include the embedding coefficients \(\tilde{b}^E\) and \(\tilde{b}^I\), although users may provide their own tables without embeddings (in which case adaptivity will be disabled). For mixed stiff/nonstiff problems, a user should provide both of the functions \(f^E\) and \(f^I\) that define the IVP system, in which case ARKStep provides the ImEx-ARK Butcher tables proposed in [35, 37] with orders of accuracy \(q = \lbrace 3,4,5\rbrace\) and embeddings of orders \(p = \lbrace 2,3,4\rbrace\), respectively.

Either of the function pointers for \(f^I\) or \(f^E\) may be NULL to disable that component. For nonstiff problems one can disable \(f^I\), in which case Equation (5) reduces to the non-split IVP: (8) \(\begin{equation} M(t)\, y^{\prime }(t) = f^E(t,y), \qquad y(t_0) = y_0. \end{equation}\) In this scenario, the implicit Butcher table (\(A^I\), \(b^I\), \(\tilde{b}^I\), \(c^I\)) in Equation (7) is ignored, and the ARK methods reduce to classical ERK methods. Here, all methods from ERKStep are again available in ARKStep. We note that in this mode, both the problem in Equation (8) and the method in Equation (7) fully encapsulate the ERKStep problem in Equation (3) and method in Equation (4). While it therefore follows that ARKStep can be used to solve every problem solvable by ERKStep (even with the same Butcher tables), we retain ERKStep as a distinct time-stepping module since its simplified form admits a more efficient and memory-friendly implementation than that provided by ARKStep.

Alternately, for stiff problems \(f^E\) may be disabled so that Equation (5) reduces to (9) \(\begin{equation} M(t)\, y^{\prime }(t) = f^I(t,y), \qquad y(t_0) = y_0. \end{equation}\) Here the explicit Butcher table (\(A^E\), \(b^E\), \(\tilde{b}^E\), \(c^E\)) in Equation (7) is ignored, and the ARK methods reduce to standard diagonally implicit Runge–Kutta(DIRK) methods, for which ARKODE provides tables with orders of accuracy \(q = \lbrace 2,3,4,5\rbrace\) and embeddings of orders \(p = \lbrace 1,2,3,4\rbrace\), respectively.

As with ERKStep, the ARKStep module includes a number of features that allow users to optimize its usage on their problems. The full list of these is included in [48]. One such feature is the routine ARKStepSetLinear, which allows the user to specify that \(f^I\) depends linearly on y, and thus each implicit stage equation for \(z_i\) from Equation (7) corresponds to a linear system of equations. In this case, only one iteration of the Newton nonlinear solver is performed, avoiding extraneous nonlinear iterations and residual norm computations to assess convergence.

A second notable feature in ARKStep is its native interface to the scalable parallel-in-time (PinT) library XBraid [12]. Instead of performing a standard time-marching approach from one step to the next \((t_{n-1}\rightarrow t_{n})\), XBraid solves for all time values simultaneously using multigrid reduction in time(MGRIT) [13], a highly parallel iterative method that exposes parallelism in the time domain in addition to spatial parallelization. We have designed the ARKStep + XBraid interface so that simulation codes using ARKStep should require minimal modifications to explore using PinT.

2.3 MRIStep

MRIStep is the newest time-stepping module in ARKODE and targets multirate IVPs of the form (10) \(\begin{equation} y^{\prime }(t) = f^{E}(t,y) + f^{I}(t,y) + f^F(t,y), \qquad y(t_0) = y_0. \end{equation}\) The right-hand side function is split into slow components, \(f^E(t,y) + f^I(t,y)\), that should be integrated with a large step size H, and where \(f^E(t,y)\) contains nonstiff terms and \(f^I(t,y)\) contains stiff terms, and fast components, \(f^F(t,y)\), that should be integrated with a small step size \(h \ll H\).

For approximating solutions to Equation (10), MRIStep utilizes multirate infinitesimal (MRI) methods, which are characterized by an alternation between an outer ARK-based method for the slow time scale and a sequence of fast time scale IVP systems that are evolved using an inner solver.

When both slow components \(f^E\) and \(f^I\) are present, MRIStep uses implicit-explicit multirate infinitesimal generalized-structure additive Runge–Kutta (ImEx-MRI-GARK) methods [8]. When either \(f^E\) or \(f^I\) are NULL-valued, ImEx-MRI-GARK methods simplify to either implicit or explicit multirate infinitesimal GARK (MRI-GARK) methods [51]. The following algorithm outlines a single step \(t_{n-1}\rightarrow t_{n-1}+H=t_n\) of an s-stage MRI method from either family:

Here the coefficients, \(0=c_1 \le c_s= 1\), define the abscissae for the method. The coefficient functions, \(\omega _{i,j}\) and \(\gamma _{i,j}\), are polynomials in time that dictate the couplings from the slow to the fast time scale and can be expressed as in [8] as \(\begin{equation*} \omega _{i,j}(\theta) = \sum _{k=0}^K \omega _{i,j}^{\lbrace k\rbrace }\theta ^k, \quad \gamma _{i,j}(\theta) = \sum _{k=0}^K \gamma _{i,j}^{\lbrace k\rbrace } \theta ^k, \end{equation*}\) where typically \(0\le K \le 2\). The coupling tables, \(\Omega ^{\lbrace k\rbrace }\in \mathbb {R}^{s\times s}\) and \(\Gamma ^{\lbrace k\rbrace }\in \mathbb {R}^{s\times s}\), contain the explicit, \(\omega _{i,j}^{\lbrace k\rbrace }\), and implicit, \(\gamma _{i,j}^{\lbrace k\rbrace }\), coefficients, respectively.

Like Runge–Kutta methods, implicitness at the slow time scale is characterized by nonzero values on or above the diagonal. MRIStep currently supports diagonally implicit, “solve-decoupled” methods, i.e., \(\omega _{i,j}^{\lbrace k\rbrace } = 0\) for \(j \ge i\), \(\gamma _{i,j}^{\lbrace k\rbrace }=0\) for \(j\gt i\), and \(\Delta c_i=0\) if \(\gamma _{i,i}^{\lbrace k\rbrace } \ne 0\). When \(\Delta c_i=0\), the “fast” IVP solve in step 2b reduces to a standard ImEx-ARK stage in Equation (7a), with step size H and coefficients \(\begin{equation*} A^E_{i,j} = \sum _{k=0}^K \frac{\omega _{i,j}^{\lbrace k\rbrace }}{k+1}, \quad A^I_{i,j} = \sum _{k=0}^K \frac{\gamma _{i,j}^{\lbrace k\rbrace }}{k+1}. \end{equation*}\) Thus, an MRI stage only requires a nonlinear implicit solve if \(A^I_{i,i} \ne 0\).

MRIStep also supports multirate infinitesimal step (MIS) methods [53, 54, 55]. These methods are a subset of MRI-GARK methods where the coupling coefficients are uniquely defined based on a slow explicit Butcher table \((A^E,b^E,c^E)\) with \(s-1\) stages and sorted abscissae \(c^E\). The MIS method then has abscissae \(c = [c^E_1 \cdots c^E_{s-1} 1]^T\) and coupling coefficients \(\gamma _{i,j}^{\lbrace 0\rbrace } = 0\) and (11) \(\begin{equation} \omega _{i,j}^{\lbrace 0\rbrace } = {\left\lbrace \begin{array}{ll} 0, & \text{if} i=1,\\ A^E_{i,j} - A^E_{i-1,j}, & \text{if} 2\le i\le s-1,\\ b^E_{j} - A^E_{s-1,j}, & \text{if} i=s. \end{array}\right.} \end{equation}\)

At present, the MRIStep module implements all third- and fourth-order ImEx-MRI-GARK methods from [8], the third-order MIS method defined by the slow Butcher table from [38], and all of the third- and fourth-order MRI-GARK methods from [51]. Alternately, users may provide a custom set of coupling tables, \(\Gamma ^{\lbrace k\rbrace }\) and/or \(\Omega ^{\lbrace k\rbrace }\), to define their own MRI method. MRIStep also includes a utility routine that will convert a slow Butcher table \((A^E,b^E,c^E)\) with sorted abscissae \(c^E\) into a set of MRI-GARK coefficients using Equation (11). If the slow Butcher table has order of accuracy at least two, then the MIS method will also have order two. Third-order accuracy can be obtained if the slow table itself is at least third order and satisfies the condition [38] \(\begin{equation*} \sum _{i=2}^{\hat{s}} \left(c_i^E-c_{i-1}^E\right) \left(e_i+e_{i-1}\right)^T A^E c^E + \left(1-c_{\hat{s}}^E\right) \left(\frac{1}{2}+e_{\hat{s}}^T A^E c^E\right) = \frac{1}{3}, \end{equation*}\) where \(A^{E}\) has \(\hat{s}\) stages and \(e_k\) is the kth column of an \(\hat{s}\times \hat{s}\) identity matrix.

The IVPs in step 2b may be solved using any sufficiently accurate method by supplying a “fast IVP” stepper to the MRIStep module. For convenience, ARKODE includes a constructor that will attach an ARKStep instance as the inner integrator, thereby enabling an adaptive explicit, implicit, or ImEx treatment of the fast time scale. Alternatively, applications can supply a user-defined integrator for the fast time scale as demonstrated in the simple example problem in Section 5.1.3 and in the large-scale demonstration problem in Section 5.2. Currently MRIStep supports temporal adaptivity only at the fast time scale. Efficient multirate temporal adaptivity is an open research question (see, e.g., [14]), but future support is anticipated to follow new developments in the field. Finally, like ERKStep and ARKStep, MRIStep provides numerous configuration options for users to control integrator behavior. The full list of options is included in [48].

3.1 User-supplied Tolerances and Error Control

To control errors at various levels (time integration and algebraic solvers), ARKODE, like other SUNDIALS solvers [32, Equation (4)], uses a weighted root-mean-square norm for all error-like quantities: (12) \(\begin{equation} \Vert v\Vert _{\text{WRMS}} = \left(\frac{1}{N} \sum _{i=1}^N \left(v_i\,w_i\right)^2\right)^{1/2}. \end{equation}\) This norm allows problem-specific error measurement via the weighting vector, w, that combines the units of the problem with user-supplied values that specify “acceptable” error levels. Thus, at the start of each time step \(t_{n-1} \rightarrow t_n\), we construct two weight vectors. The first is an error weight vector that, as in CVODE [32, Equation (5)], uses the most recent step solution and user-supplied relative (rtol) and absolute (atol, scalar- or vector-valued) tolerances: (13) \(\begin{equation} w_i = \left(\text{rtol}\, |y_{n-1,i}| + \text{atol}_{i}\right)^{-1}. \end{equation}\) The second focuses on problems with a non-identity operator \(M(t)\), since the units of Equation (1) may differ from the units of the solution y. Here, we construct a residual weight vector: (14) \(\begin{equation} \sigma _i = \left(\text{rtol}\, |r_i| + \text{ratol}_i\right)^{-1}, \quad r = M(t_{n-1}) y_{n-1}, \end{equation}\) where the user may specify a separate residual absolute tolerance value or vector, \(\text{ratol}\).

ARKODE uses w or \(\sigma\) based on the quantity being measured: errors having “solution” units use w, whereas those having “equation” units use \(\sigma\). Obviously, when \(M=I\), the solution and equation units match, and ARKODE uses w for all error norms. In both cases, since \(w_i^{-1}\) (or \(\sigma _i^{-1}\)) represents a tolerance in the ith component of the solution (or equation) vector, then an error whose WRMS norm is \(\le \!\! 1\) is regarded as being within an acceptable error tolerance.

3.2 Stepsize Adaptivity

A critical utility provided by ARKODE is its adaptive control of local truncation error (LTE), \(T_n\), at each step, which corresponds to the error induced within the time step \(t_{n-1}\rightarrow t_n\), under an assumption that the initial condition for that step was exact. At every internal step, ARKODE requests both an approximate solution \(y_n\) and an estimate of \(T_n\) from the time-stepping module. These approximations have global orders of accuracy q and p, respectively, where generally \(p=q-1\). Since the global and local errors for one-step methods usually differ by one factor of \(h_n\) [25], then (15) \(\begin{align} \Vert y_n - y(t_n) \Vert = C h_n^{q+1} + \mathop {}\mathopen {}\mathcal {O}\mathopen {}\left(h_n^{q+2}\right)\!, \end{align}\) (16) \(\begin{align} \Vert T_n \Vert = D h_n^{p+1} + \mathop {}\mathopen {}\mathcal {O}\mathopen {}\left(h_n^{p+2}\right)\!, \end{align}\) where C and D are constants independent of \(h_n\), and where we have assumed exact initial conditions for the step, i.e., \(y_{n-1} = y(t_{n-1})\). Given this time-stepper-estimated value of \(T_n\), ARKODE adopts the same local temporal error test as other SUNDIALS integrators, simply (17) \(\begin{equation} \Vert T_n\Vert _{\text{WRMS}} \le 1. \end{equation}\) If this test passes, the step is accepted and \(T_n\) is used to estimate a prospective next step size, \(h^{\prime }\), using an adaptive error controller. If Equation (17) fails, then the step is rejected, and \(T_n\) is used to estimate a reduced step size, \(h^{\prime }\). In either case, temporal adaptivity in ARKODE focuses on choosing the maximum \(h^{\prime }\) such that Equation (17) should pass on the next attempted step. We thus define the biased local error estimate: \(\begin{equation*} \varepsilon _k \ \equiv \ \beta \, \Vert T_k\Vert _{\text{WRMS}}, \end{equation*}\) where the error bias \(\beta \gt 1\) helps to account for the error constant D in Equation (16). ARKODE stores the biased error estimates corresponding to the three most recent steps, \(\varepsilon _n\), \(\varepsilon _{n-1}\), and \(\varepsilon _{n-2}\). Of these, both \(\varepsilon _{n-1}\) and \(\varepsilon _{n-2}\) correspond to the last two successful steps, whereas \(\varepsilon _n\) corresponds to the most recently attempted step. ARKODE provides a variety of error controllers that utilize these estimates, including PID (the ARKODE default) [35, 57, 58, 59], PI [35, 57, 58, 59], I (the “standard” method used in most publicly available IVP solvers) [33], the explicit Gustafsson controller [23], the implicit Gustafsson controller [24], and an “ImEx” Gustaffson controller that sets \(h^{\prime }\) as the minimum of the explicit and implicit Gustafsson controller values. Alternately, users may provide their own functions to produce \(h^{\prime }\) using the inputs \(y_n, t_n, h_n, h_{n-1}, h_{n-2}, \varepsilon _n, \varepsilon _{n-1}\), and \(\varepsilon _{n-2}\). Each controller bases its adaptivity off the order of accuracy for the error estimate p, although the user may request to use the solution order q instead. All of the provided controllers include tuning parameters, with default values determined using an exhaustive genetic optimization process based off of dozens of challenging IVP test problems, including nearly all from the “Test Set For IVP Solvers” [62].

In addition to error-based temporal adaptivity, ARKODE adopts many of the heuristics from CVODE [32] for bounding \(h^{\prime }\), including approaches for handling repeated temporal error test failures, overly aggressive growth or decay in prospective steps (\(h^{\prime } \gg h_n\) or \(h^{\prime } \ll h_n\)), step size reduction following a failed implicit solve, and user-supplied step size bounds \(h_\text{min} \le h^{\prime } \le h_\text{max}\). Complete details on these heuristics, as well as information on how each parameter may be adjusted by users, are provided in the ARKODE documentation [48].

3.3 Interpolation Modules

ARKODE supports interpolation of approximate solutions \(y(t_{out})\) and derivatives \(y^{(d)}(t_{out})\), where \(t_{out}\) occurs within the most recently completed step \(t_{n-1}\rightarrow t_n\). This utility also supports extrapolation of y and its derivatives outside this time step to construct predictors for iterative algebraic solvers. These approximations are based on polynomial interpolants \(\pi _{\xi }(t)\) of degree up to \(\xi =5\) that are constructed using one of two complementary approaches: “Hermite” and “Lagrange.”

3.4 Implicit Solvers

For both ARKStep and MRIStep, if \(f^{I}\) is nonzero, then each implicit stage, \(z_i\in \mathbb {R}^N\), may require the solution of an algebraic system of equations: (18) \(\begin{equation} M(t_{n,i})\left(z_i - a_i\right) - \gamma _i f^I(t_{n,i},z_i) = 0, \end{equation}\) where \(t_{n,i}\) is the corresponding implicit stage time, \(\gamma _i\in \mathbb {R}\) is proportional to \(h_n\), \(f^I(t,y)\) is the implicit portion of the IVP right-hand side, and \(a_i\) is “known data” that may include previous stages or IVP right-hand side evaluations. ARKODE provides utilities to map between \(f^I\) and the nonlinear or linear algebraic system for each stage (18), interfaces to a range of algebraic solvers from SUNDIALS, and optimizations to enhance the efficiency of these solvers for the IVP context.

3.5 Temporal Root-finding

A particularly useful feature of SUNDIALS integrators that is also provided by the ARKODE infrastructure is event detection. While integrating an IVP, ARKODE can find roots of a set of user-defined functions, \(g_i(t,y(t)) = \tilde{g}_i(t)\). The number of these functions is arbitrary, and if more than one \(\tilde{g}_i\) is found to have a root in any given interval, the various root locations are found and reported in the order that they occur in the direction of integration. The basic scheme and heuristics match other SUNDIALS integrators [32]. During integration of Equation (1), following each successful step, ARKODE checks for sign changes of any \(\tilde{g}_i(t)\); if one is found, it then will home in on the root (or roots) with a modified secant method [28].

3.6 Inequality Constraints

A final notable feature of the ARKODE infrastructure, also shared with other SUNDIALS integrators, is support for user-imposed inequality constraints on components of the solution y: \(y_i \gt 0\), \(y_i \lt 0\), \(y_i \ge 0\), or \(y_i \le 0\). Constraint satisfaction is tested after a successful step and before the temporal error test. If any constraint fails, ARKODE estimates a reduced step size based on a linear approximation in time for each violated component (including a safety factor to cover the strict inequality case). If applicable, a flag is set to update the Jacobian or preconditioner in the next step attempt.

5.1 Simple Example Problem

For each stepper module example, we consider a one-dimensional advection-diffusion-reaction equation using a variation of the Brusselator chemical kinetics model [26]. The system is given by (22) \(\begin{equation} \begin{aligned}u_t &= -c\, u_x + d\, u_{xx} + a - (w + 1)\, u + v\, u^2, \\ v_t &= -c\, v_x + d\, u_{xx} + w\, u - v\, u^2, \\ w_t &= -c\, w_x + d\, u_{xx} + (b - w) / \epsilon - w\, u, \end{aligned} \end{equation}\) for \(t \in [0, 10]\) with \(x \in [0, 1]\), where u, v, and w are the chemical species concentrations; \(c = 0.001\) is the advection speed; d is the diffusion rate (typically 0.01); \(a = 0.6\) and \(b = 2\) are the concentrations of species that remain constant over time; and \(\epsilon = 0.01\) is a parameter that determines the stiffness of the system. The initial conditions are \(\begin{equation*} \begin{gathered}u(0,x) = a + 0.1 \sin (\pi x), \qquad v(0,x) = b/a + 0.1 \sin (\pi x), \qquad w(0,x) = b + 0.1 \sin (\pi x). \end{gathered} \end{equation*}\) Spatial derivatives are computed using second-order centered differences with the data distributed over \(N = 512\) points on a uniform spatial grid with stationary boundary conditions, i.e., \(\begin{equation*} \begin{gathered}u_t(t,0) = u_t(t,1) = 0, \qquad v_t(t,0) = v_t(t,1) = 0, \qquad w_t(t,0) = w_t(t,1) = 0. \end{gathered} \end{equation*}\) In the subsections that follow, we explore various time integration methods and ImEx splittings of Equation (22), so we first summarize the relevant time scales present in this problem for each of the advection, diffusion, and reaction terms. For this problem setup, if explicit time-stepping methods were used for all operators and linear stability were the only consideration, then diffusion would require the smallest time step, followed by reaction, followed by advection. However, if linear stability were not an issue (e.g., when using A-stable implicit methods) and temporal accuracy were the only consideration, then reaction would require the smallest time steps, followed by advection and diffusion, which would evolve on comparable time scales. Finally, we note that the \((b-w)/\epsilon\) term in Equation (22) induces a mild amount of additional stiffness into the reaction processes, thus requiring slightly smaller “stable” than “accurate” steps.

All of the examples use the serial N_Vector implementation and results were obtained on a single node of the Quartz cluster at Lawrence Livermore National Laboratory. Each Quartz node contains two Intel Xeon 5-2695 v4 CPUs. The source files for this example problem are ark_advection_diffusion_reaction.hpp and ark_advection_diffusion_reaction.cpp. In each test case the maximum relative error at the final time is computed using a reference solution generated with the default fifth-order DIRK method (ARK5(4)8L[2]SA-ESDIRK in [35]) with relative and absolute tolerances of \(10^{-8}\) and \(10^{-14}\), respectively.

5.1.1 ERKStep.

We first consider an advection-reaction setup (i.e., \(d = 0\) in Equation (22)) and evolve the system in time using the default second- (Heun-Euler), third- [5], fourth- [65], and fifth-order [6] ERK methods in ARKODE. The optimal method choice will depend on the user’s desired accuracy, the cost per step of the method, and the step sizes selected by the adaptivity controller. As such, we show how the options in ARKODE can help users determine more effective methods in their accuracy regimes of interest. We compare the performance of each method using “loose,” “medium,” and “tight” sets of relative and absolute tolerances pairs (\(10^{-4}\)/\(10^{-9}\), \(10^{-5}\)/\(10^{-10}\), and \(10^{-6}\)/\(10^{-11}\), respectively) and the PID, PI, I, and explicit Gustafsson controllers.

Figure 2 shows a work-precision plot for the various configurations where the total number of right-hand side function evaluations for successful and failed steps is used as the work metric. Given its lower cost per step, the second-order method gives the best performance at coarser levels of accuracy (three left-most solid red markers) despite requiring more than double the number of step attempts as the third-order method. For errors between \(10^{-5}\) and \(10^{-7}\) the third-order method is the most efficient (lower left orange markers) since the second-order method requires at least three times as many step attempts, while the fourth- and fifth-order methods make fewer attempts but the reduced number of attempts is insufficient to offset their higher cost per step. The fifth-order method becomes most efficient at accuracies below \(10^{-7}\) (blue markers near the bottom left).

Fig. 2.

View Figure

Across the methods and tolerances, the I controller is frequently the least efficient for this example as it only considers the error in the current step leading to a high step rejection rate (between 13% and 50%). The additional error history retained by the PI and Gustafsson controllers leads to better step size predictions and often the greatest efficiency (rejection rates <7%). The longer history of the PID controller does not lead to further improvement and it generally produces results that fall between the PI and I controllers. As explored in [45], the choice of controller parameters (tunable in ARKODE but not explored here) also plays an important role in performance.

5.2 Large-scale 3D Demonstration Problem

As a final demonstration of ARKODE’s flexibility and performance, we consider the three-dimensional nonlinear inviscid compressible Euler equations, with advection and reaction of chemical species: (23) \(\begin{equation} \mathbf {w}_t = -\nabla \cdot F(\mathbf {w}) + R(\mathbf {w}), \end{equation}\) with the independent variables \((X,t) = (x,y,z,t) \in \Omega \times [0, t_f]\), and where the spatial domain is a three-dimensional cube, \(\Omega = [0, 1] \times [0, 1] \times [0, 1]\). The solution is given by \(w = [ \rho \rho v_x \rho v_y \rho v_z e_t \mathbf {c} ]^T\), corresponding to the density; momentum in the x, y, and z directions; total energy per unit volume; and some number of chemical densities \(\mathbf {c}\in \mathbb {R}^{nchem}\) that are advected along with the fluid. The fluxes are given by (24) \(\begin{align} F_x(w) &= \begin{bmatrix} \rho v_x & \rho v_x^2 + p & \rho v_x v_y & \rho v_x v_z & v_x (e_t+p) & \mathbf {c} v_x \end{bmatrix}^T, \end{align}\) (25) \(\begin{align} F_y(w) &= \begin{bmatrix} \rho v_y & \rho v_x v_y & \rho v_y^2 + p & \rho v_y v_z & v_y (e_t+p) & \mathbf {c} v_y \end{bmatrix}^T, \end{align}\) (26) \(\begin{align} F_z(w) &= \begin{bmatrix} \rho v_z & \rho v_x v_z & \rho v_y v_z & \rho v_z^2 + p & v_z (e_t+p) & \mathbf {c} v_z \end{bmatrix}^T. \end{align}\) The reaction mechanism is encoded in \(R(w)\), and the ideal gas equation of state gives \(p = (\gamma -1) (e_t - \frac{\rho }{2} (v_x^2 + v_y^2 + v_z^2))\), where \(\gamma = c_p/c_v\) is the constant ratio of specific heats for the gas.

We discretize this problem using the method of lines, where we first semi-discretize in space using a regular finite volume grid with dimensions \(n_x \times n_y \times n_z\), with fluxes at cell faces calculated using a fifth-order FD-WENO reconstruction [56]. MPI parallelization is achieved using a standard 3D domain decomposition of \(\Omega\). We organize the corresponding spatially discretized solution vector using the SUNDIALS MPIManyVector [18] that wraps node-local vectors for each of the fields, \(\rho , \ldots , \mathbf {c}\), to create the overall state vector, \(\mathbf {w}\), and provides the requisite MPI functionality for coordinating vector operations among the subvectors. The hydrodynamic fields (\(\rho\), \(\rho \mathbf {v}\), and e) are stored in CPU memory, and the chemical densities, \(\mathbf {c}\), are stored in GPU memory. After spatial semi-discretization, we are faced with a large IVP system: (27) \(\begin{equation} {w}^{\prime }(t) = f_1({w}) + f_2({w}), \quad {w}(t_0)={w}_0, \end{equation}\) where \(f_1({w})\) and \(f_2({w})\) contain the spatially discretized forms of \(-\nabla \cdot F(\mathbf {w})\) and \(R(\mathbf {w})\), respectively. For additional information on this test problem or our computational implementation, please see the technical report [47] and the GitHub repository [46, v4.0].

In the results that follow, we simulate the advection and reaction of a low-density primordial gas, present in models of the early universe [1, 21, 39, 61]. This gas consists of 10 advected species (i.e., \(nchem=10\)), corresponding to various ionization states of atomic and molecular Hydrogen and Helium, free electrons, and the internal gas energy. For this setup the reaction mechanism is considerably stiffer than the advection, so we consider two forms for temporal evolution:

For (a) we use SUNDIALS’s modified Newton solver to handle the global nonlinear algebraic systems arising at each implicit stage of each time step. Since only \(f_2\) is treated implicitly and the reactions are purely local in space, the Newton linear systems are block-diagonal. As such, we provide a custom SUNLinearSolver implementation that solves each MPI rank-local linear system independently. The portion of the Jacobian matrix on each rank, \(J_p\), is itself block-diagonal, i.e., \(\begin{equation*} J_p = \begin{bmatrix} J_{p,1,1,1} & & & \\ & J_{p,2,1,1} & & \\ & & \ddots & \\ & & & J_{p,n_{xloc},n_{yloc},n_{zloc}} \end{bmatrix}, \end{equation*}\) with each cell-local block \(J_{p,i,j,k} \in \mathbb {R}^{10\times 10}\). We further leverage this structure by solving each \(J_p x_p = b_p\) linear system using the SUNDIALS SUNLinearSolver implementation that interfaces to the GPU-enabled, batched, dense LU solver routines from MAGMA [60].

The multirate approach (b) can leverage the structure of \(f_2\) at a higher level. Since the MRI method applied to this problem evolves “fast” sub-problems of the form (28) \(\begin{equation} v^{\prime }(t) = f_2(t,v) + r_i(t), \quad i=2,\ldots ,s, \end{equation}\) and all MPI communication necessary to construct the forcing functions, \(r_i(t)\), has already been performed, each sub-problem in Equation (28) consists of \(n_x\,n_y\,n_z\) spatially decoupled fast IVPs. We construct a custom fast integrator that groups all the independent fast IVPs on an MPI rank together as a single system evolved using a rank-local ARKStep instance. The code for this custom integrator itself is minimal, primarily consisting of steps to access the local subvectors in w on a given MPI rank and wrapping them in MPI-unaware ManyVectors provided to the local ARKStep instance. The collection of independent local IVPs also leads to a block-diagonal Jacobian, and we again utilize the SUNDIALS interface to MAGMA for linear solves within the modified Newton iteration.

These tests use periodic boundary conditions and a “clumpy” initial background density field: \(\begin{equation*} \rho _b(X) = 100 \left(1 + \sum _{i=1}^{N_c} s_i \exp \left(-2 \left(\frac{\Vert X-X_i\Vert }{r_i}\right)^2\right)\right), \end{equation*}\) where \(s_i \in [0,5]\), \(r_i \in [3\Delta x,6\Delta x]\), \(X_i\in \Omega _p\) are uniformly distributed random values, and there are 10 “clumps” per MPI rank. The initial background temperature is \(T_b(X)= 10\), and the initial velocity field is identically zero. The initial conditions for density and temperature are obtained by adding a Gaussian bump to the background states: \(\begin{equation*} \rho _0(X) = \rho _b(X) + 500 \exp \left(-2 \left(\frac{\Vert X-0.5\Vert }{0.1}\right)^2 \right), \quad T_0(X) = T_b(X) + 50 \exp \left(-2\left(\frac{\Vert X-0.5\Vert }{0.1}\right)^2 \right), \end{equation*}\) causing a high-pressure front that emanates from the center of the domain. We initialize the chemical species to consist almost exclusively of neutral Hydrogen (76%) and Helium (24%), with trace amounts of ionized Hydrogen and Helium (0.0000000001%) outside of this high-pressure region. As the region expands and increases the surrounding temperature, the chemical concentrations in those cells rapidly transition to a different equilibrium, resulting in a chemical time scale that is multiple orders of magnitude faster than the spatial propagation of the high-pressure front.

In Figure 3 we present timing results for this problem run on Summit at the Oak Ridge Leadership Computing Facility. Each Summit node contains two IBM POWER9 processors and six NVIDIA Tesla V100 GPUs. The nodes are connected to a dual-rail EDR InfiniBand network in a non-block fat tree topology. We performed a weak-scaling simulation with:

Fig. 3.

View Figure

• a \(25 \times 25 \times 25\) spatial grid per MPI rank,

• multirate integrators using a slow time step size satisfying the CFL constraint \(H \le 10\Delta x\),

• fast multirate time step sizes and ImEx time step sizes chosen adaptively by the ARKStep integrator with relative and absolute tolerances of \(10^{-5}\) and \(10^{-9}\), respectively,

• temporal evolution over \([0,5H]\), and

• six MPI ranks per Summit node, with one NVIDIA Volta GPU tied to each MPI rank.

We note that both the ImEx and MRI approaches select time steps to accurately track the chemical evolution. In the ImEx case, the entire simulation evolves with a shared step that is therefore set by the chemistry in the “hardest” spatial cell. However, in the multirate case, only the fast time scale must track the chemistry, and that occurs on a rank-by-rank basis, allowing other MPI ranks and other physical processes to evolve with larger time steps.

The top row of Figure 3 shows the total simulation times for both configurations of the problem, while the second row shows the corresponding algorithmic scalability. Examining the second row, we see that for all solver statistics both configurations demonstrate ideal algorithmic scalability, implying that our algorithm choices map optimally to the underlying problem structure.

From the first row of Figure 3 two observations are immediately clear. First, the high cost associated with each evaluation of the advection operator (shown as fAdv) causes the MRI approach to easily outperform the ImEx method, due to its much smaller number of evaluations. Second, the primary runtime slowdown experienced within both approaches arises in the green MPI component. For the ImEx method, this component corresponds solely to a single MPI_Allreduce after each rank-local linear solve to determine the overall “success” or “failure” of the global linear solve. Thus, the MPI cost reflects the cost of global synchronization. Similarly, in the MRI case the MPI component measures an MPI_Allreduce that determines success of the fast local IVP solves. Hence, the fundamental scalability difference between the ImEx and MRI approaches results from the frequency of these synchronizations. For the ImEx method synchronizations occur following each linear solve, within each Newton iteration, within each implicit stage, within each time step. Since ARKStep’s temporal adaptivity tracked the (fast) reaction time scale, the simulation required a large number of synchronizations. On the other hand, with the MRI method these synchronizations occur at the slow time scale, i.e., considerably less often than in the ImEx case.

Taking a step back to focus on ARKODE as an infrastructure, we conclude that its algorithmic flexibility, and in particular its support for custom solver components and the new MRIStep multirate module, enables a rich exploration of the “solver space” for this and many problems.