Performance of explicit and IMEX MRI multirate methods on complex reactive flow problems within modern parallel adaptive structured grid frameworks

The first class of multirate schemes that we consider is the implicit-explicit multirate infinitesimal GARK (IMEX-MRI-GARK) methods from Chinomona and Reynolds (2021). These schemes are a recent addition to the broader class of infinitesimal methods, that generalize the previously-mentioned MIS and MRI-GARK methods constructed for problems of the form (1). Methods within this family share two fundamental traits. First, they utilize a Runge–Kutta like structure for each time step $t_{n}\to t_{n}+H$, wherein the overall time step is achieved through computation of $s$ internal stages that are eventually combined to form the overall time step solution $y_{n+1}\approx y(t_{n}+H)$. These internal stages roughly correspond to approximations $Y_{i}\approx y(t_{n}+c_{i}H)$, $i=1,\ldots,s$, where a specific scheme includes an array of coefficients $c\in\mathbb{R}^{s}$ that define these stage times. Second, these methods couple the fast and slow time scales by solving a set of modified fast IVPs to obtain the internal stages. Specifically, to compute the updated solution $y_{n+1}$ from $y_{n}$, these methods perform the algorithm:

1. 1.

   Set $Y_{1}=y_{n}$ and $T_{1}=t_{n}$.

2. 2.

   For $i=2,\ldots,s$:

   • a.

     Let $v(T_{i-1})=Y_{i-1}$, $T_{i}=t_{n}+c_{i}H$ and $\Delta c_{i}=c_{i}-c_{i-1}$, and define a “forcing” polynomial $r_{i}(\theta)$ constructed from evaluations of $f^{S}(T_{j},Y_{j})$ where $j=1,\ldots,i$.

   • b.

     Solve the “fast” IVP over $\theta\in[T_{i-1},T_{i}]$:

     $\displaystyle v^{\prime}(\theta)$

     $\displaystyle=f^{F}(\theta,v(\theta))+r_{i}(\theta).$

     (3)

   • c.

     Set $Y_{i}=v(T_{i})$.

3. 3.

   Set $y_{n+1}=Y_{s}$.

When $r_{i}(\theta)$ depends on $f^{S}(T_{i},Y_{i})$ the step 2b is implicit in $Y_{i}$; however when $\Delta c_{i}=0$, that step simplifies to a standard Runge–Kutta update. Thus, infinitesimal methods that allow implicit treatment of $f^{S}$ typically only allow implicit dependence on $Y_{i}$ for stages where $\Delta c_{i}=0$, resulting in a “solve-decoupled” approach. Here, the stages alternate between a standard diagonally-implicit slow nonlinear system solve and evolving the modified fast IVP. We note that since the step 2b may be solved using any sufficiently-accurate algorithm, these methods offer extreme flexibility for practitioners to employ problem-specific integrators.

Within this infinitesimal family, the MIS methods first introduced by Schlegel et al. (2009) and Schlegel (2011) treat the slow time scale operator $f^{S}$ explicitly, and result in up to $\mathcal{O}(H^{3})$ accuracy. The MRI-GARK methods introduced by Sandu (2019) treat $f^{S}$ using either an explicit or a solve-decoupled implicit scheme, and result in up to $\mathcal{O}(H^{4})$ accuracy. The IMEX-MRI-GARK methods proposed by Chinomona and Reynolds (2021) treat the slow time scale using ImEx methods, $f^{S}=f^{I}+f^{E}$, and result in up to $\mathcal{O}(H^{4})$ accuracy. Additionally, within IMEX-MRI-GARK methods, either of the components, $f^{I}$ or $f^{E}$, may be zero, resulting in an MRI-GARK method.

We note that these flexible infinitesimal formulations are not new, as the longstanding Lie-Trotter and Strang-Marchuk operator splitting methods similarly allow arbitrary algorithms at the fast time scale, as well as explicit, implicit, or IMEX treatment of the slow time scale. However, both of these approaches use the forcing function $r_{i}(\theta)=0$, resulting in $\mathcal{O}(H)$ and $\mathcal{O}(H^{2})$ accuracy, respectively.

Spectral deferred correction algorithms are a class of numerical methods that represent the solution as an integral in time and iteratively solve a series of correction equations designed to reduce the error. The SDC approach is introduced in Dutt et al. (2000) for the solution of ordinary differential equations. The correction equations are typically formed using a first-order time-integration scheme, but are applied iteratively to construct schemes of arbitrarily high accuracy. In practice, the time step is divided into subintervals, and the solution at the end of each subinterval is iteratively updated using (typically) forward or backward Euler temporal discretizations. Each set of correction sweeps over all the subintervals increases the overall order of the method by one, up to the underlying order of the numerical quadrature rule over the subintervals. Thus, it is advantageous to define the subintervals using high-order numerical quadrature rules, such as Gauss-Lobatto, Gauss-Radau, Gauss-Legendre, or Clenshaw-Curtis.

There have been a number of works that have extended the SDC approach to handle physical systems with multiple disparate time scales. Minion et al. (2003) introduces a semi-implicit version for ODEs with stiff and non-stiff processes, such as advection–diffusion systems. The correction equations for the non-stiff terms are discretized explicitly, whereas the stiff term corrections are treated implicitly. Bourlioux et al. (2003) introduces a multirate SDC (MRSDC) approach for PDEs with advection-diffusion-reaction processes where advection terms are evaluated explicitly, reaction and diffusion terms treated implicitly, and different time steps are used for each process. For example, the “slow” advection process can be discretized using Gauss-Lobatto quadrature over the entire time step while the stiffer diffusion/reaction processes can be sub-divided into Gauss-Lobatto quadrature(s) within each of the advection nodes. This approach was also used successfully by Layton and Minion (2004) for a conservative formulation of reacting compressible gas dynamics. More recently, approaches have been developed for fourth-order, finite volume, adaptive mesh simulations (Emmett et al., 2019), as well as eighth-order finite differences (Emmett et al., 2014). The latter paper, known as the “SMC” algorithm, forms the basis of the comparisons in this paper to the aforementioned MRI integrators.

Multirate SDC offers flexibility in (i) the number of physical processes and how they are grouped together, (ii) the choices of intermediate quadrature points each process uses, and (iii) the form of the correction equation for each process. A long review can be found in the references above and mentioned within. Here, we will describe one particular approach taken in the SMC algorithm that forms the basis of our comparisons that follow.

For the MRSDC methods in this paper we consider Gauss-Lobatto quadrature, noting that both endpoints of any given interval are included in the quadrature. We divide the overall time step into a series of coarse substeps using nodes at time $t_{p}$, and each coarse substep is further divided into a series of finer substeps using nodes at time $t_{q}$. With this choice of quadrature each coarse node at time $t_{p}$ corresponds to a location where a fine node also exists. For comparison against the MRI schemes in this paper, we group advection and diffusion together as a “slow” explicit process, $f^{S}$, and evaluate these on the coarse nodes. We treat reactions as the “fast” explicit process, $f^{F}$, and evaluate their terms on fine nodes. Other configurations are possible e.g., mixed implicit-explicit methods for the slow time scale with diffusion treated implicitly in $f^{I}$ (Bourlioux et al., 2003). To describe the temporal scheme, we use the terminology “MRSDC-XYZ”, where X is the number of coarse nodes, and we apply a Y-node quadrature Z times in each coarse subinterval. As an example, Figure 1 shows a graphical illustration of an MRSDC-332 scheme with the 3-point Gauss-Lobatto quadrature, which consists of 3 equally-spaced points. Note that in general, there are $n_{q}=X+(X-1)(Z-1)+(X-1)(Y-2)Z$ total fine nodes.

The MRSDC scheme begins by supplying an initial guess for the solution at all nodes at time $t_{q}$ (typically the solution at $t_{n}$ is used), denoted $y_{q}^{k}$, where $k$ is the iteration index. Each MRSDC iteration consists of a sweep over all the fine nodes at time $t_{q}$ in order, updating the solution from $y_{q}^{k}$ to $y_{q}^{k+1}$. Each iteration in $k$ increases the overall order of accuracy by one (when using forward or backward Euler), up to the order of accuracy of the underlying quadrature. When all processes are treated explicitly, the forward-Euler correction equation takes the form

where $H_{q}=t_{q+1}-t_{q}$, and $p$ is the index of the coarse node closest to the left of time $t_{q}$, $S_{F}^{k,q}$ is the numerical integral over $[t_{q},t_{q+1}]$ of $f^{F}(t,y^{k})$ using values from the fine nodes, and $S_{C}^{k,q}$ is the numerical integral over $[t_{q},t_{q+1}]$ of $f^{S}(t,y^{k})$ using values from the coarse nodes. In different formulations, some or all of the processes can be treated implicitly; for example, replacing $f^{F}(t_{q},y_{q}^{k+1})-f^{F}(t_{q}.y_{q}^{k})$ with $f^{F}(t_{q+1},y_{q+1}^{k+1})-f^{F}(t_{q+1}.y_{q+1}^{k})$ in (4) represents a backward-Euler discretization of $f^{F}$, forming an implicit update. Equation (4) is applied over all $n_{q}$ nodes (for $q>0$) a total of $N_{\rm sweep}$ times; for this explicit approach the computational expense is largely determined by the cost of the evaluations of $f^{S}$ and $f^{F}$. The integration matrices $S_{F}$ and $S_{S}$ are (generally inexpensive) numerical integral quadratures over a small number of points. For implicit or semi-implicit updates, the cost of each application of (4) can be much more expensive than the explicit case.

We thus see that both MRI and SDC approaches achieve high-order accuracy by coupling different (possibly multirate) physical processes via forcing terms. Due to the iterative nature of SDC approaches, the computational expense tends to be larger than MRI methods of the same order, since MRI approaches advance the solution through a series of stages in “one shot”. We demonstrate the increased computational efficiency of MRI approaches in our results. Also, it is worth noting that this difference in coupling strategies currently limits MRI approaches to at-best fourth order in time, whereas there are ostensibly no upper bounds to the temporal order of accuracy for SDC.

The spatial discretization for (5) uses a finite difference approach with a uniform, single-level grid structure as described in Emmett et al. (2014). Centered, compact eighth-order spatial stencils are used in the interior of the domain, which require four cells on each side of a point where a derivative is computed. As all calculations in our tests were done in parallel using spatial domain decomposition, four ghost cells are exchanged between each spatial partition when computing derivatives (see Section 4.1 below). Stencils near non-periodic physical domain boundaries are handled as follows. Cells four away from a physical boundary use sixth-order centered stencils. Cells three away use fourth-ordered centered stencils. At cells two away, diffusion stencils use a centered second-order stencil while advection stencils use a biased third-order stencil. For the boundary cell, diffusion uses a completely biased second-order stencil and advection a wholly biased third-order stencil. While this mixture of discretizations theoretically results in overall second-order spatial accuracy, we replicate it here to allow for direct comparisons against the results from Emmett et al. (2014).

The advection and diffusion operators of the second model, (9), are discretized using standard second-order finite-volume stencils. These stencils require that one ghost cell be exchanged between spatial partitions during parallel computation of derivatives. Test problems using this model all used periodic boundary conditions so the stencils remain second-order at physical boundaries.

In this work, we explore a third-order explicit MRI method, as well as third- and fourth-order solve-decoupled IMEX-MRI-GARK methods. The overall cost at the fast time scale for each method is commensurate with a single sweep over the time interval, while the overall cost at the slow time scale is commensurate with a standard additive Runge–Kutta IMEX (ARK-IMEX) method (Ascher et al., 1995). We compare their performance on the two model problems against single- and multirate SDC methods, as well as a standard ARK-IMEX method. The configuration and rough operation counts for a single step of each of the methods is as follows.

• •

  SDC: Explicit fourth-order single-rate SDC method that uses a 3-point quadrature and four iterations per step.

• •

  MRSDC-338: Explicit fourth-order multirate SDC method reported in Emmett et al. (2014). Uses a 3-point quadrature for $f^{S}$, a 3-point quadrature repeated 8 times per iteration for $f^{F}$, and four iterations per step.

• •

  MRSDC-352: Explicit fourth-order multirate SDC method reported in Emmett et al. (2014). Uses a 3-point quadrature for $f^{S}$, a high-order 5-point quadrature repeated only twice per iteration for $f^{F}$, and four iterations per step.

• •

  MRI3: The slow method is the third-order explicit method from Knoth and Wolke (1998) that requires four evaluations of $f^{S}$ and evolution of four modified fast IVPs (3). At the fast time scale we use explicit Runge–Kutta methods of orders two through five, corresponding to Heun, Bogacki and Shampine (1989) without use of the second-order embedding, Zonneveld (1963), and Cash and Karp (1990) respectively.

• •

  ARK-IMEX: A fourth-order single-rate ARK-IMEX Runge-Kutta method from Kennedy and Carpenter (2003). Each step requires six evaluations of $f^{E}$ and the solution of five diagonally-implicit nonlinear systems using $f^{I}$.

• •

  IMEX-MRI3: The slow method is the IMEX-MRI-GARK3b method from Chinomona and Reynolds (2021) and the fast method is Kutta’s third order method Kutta (1901). Within each slow time step, this approach requires four evaluations of $f^{E}$, the solution of three diagonally-implicit nonlinear systems using $f^{I}$, and the evolution of three modified fast IVPs (3).

• •

  IMEX-MRI4: The slow method is the IMEX-MRI-GARK4 method from Chinomona and Reynolds (2021) and the fast method is the classic fourth order Runge-Kutta method. Within each slow time step, this approach requires six evaluations of $f^{E}$, the solution of five diagonally-implicit nonlinear systems using $f^{I}$, and the evolution of five modified fast IVPs (3).

The model problem of (5) is simulated using the multirate MRSDC-338, MRSDC-352, and MRI3 methods, plus the single-rate SDC. The right hand side is kept unpartitioned for SDC but is divided into fast and slow partitions for the multirate methods (both MRSDC and MRI). Since the simulation accuracy for all approaches hinges on resolution of the chemistry, the multirate methods put the chemical reaction term in Equation (5b) into the fast partition, $f^{F}=\rho\dot{\omega}_{k}$, and the remaining transport terms are put into the slow, $f^{S}$. SDC must use the same step size for the transport as needed to accurately resolve the chemistry, whether that is necessary for them or not.

Problem (9) is computed with the multirate IMEX-MRI3 and IMEX-MRI4 methods, as well as the standard ARK-IMEX method. For the same reason as the previous problem, the multirate integrators place the chemical reaction term into the fast partition, $f^{F}=\tfrac{\dot{\omega}_{k}W_{k}}{\rho}$. The remaining slow terms are split into an explicitly stepped partition containing the advection term, $f^{E}=\nabla\cdot(Y_{k}{\bf u})$, and an implicitly stepped one with the diffusion term, $f^{I}=\nabla\cdot D_{k}\nabla Y_{k}$. For ARK-IMEX, the chemistry and advection are assigned to the explicit partition, $f^{E}$, and the diffusion to the implicit one, $f^{I}$. In that method, all partitions are stepped at the step size needed to resolve the chemistry, even if it is smaller than needed for the transport.

AMReX is a software framework that supports the development of block-structured AMR algorithms for solving systems of partial differential equations on current and emerging architectures. AMReX supplies data containers and iterators for mesh-based fields, supporting data fields of any nodality. There is no parent-child relationship between coarser grids and finer grids. See Figure 2 for a sample grid configuration. In order to create refined grids, users define a “tagging” criteria to mark regions where refined grids are necessary. The tagged cells are grouped into rectangular grids using the clustering algorithm in Berger and Rigoutsos (1991). There are user-defined parameters controlling the maximum length a grid can have, as well as forcing each side of the grid to be divisible by a particular factor. Furthermore, grids can be subdivided if there is an excess of MPI processes with no load. Proper nesting is enforced, such that a fine grid is strictly contained within the grids at the next-coarser level.

In this work we leverage (a) performance enhancements in the context of single-core and MPI performance, (b) built in grid generation, load balancing, multilevel communication patterns, (c) linear solvers via geometric multigrid, (d) efficient parallel I/O, and (d) visualization support through VisIt, ParaView, and yt. The functionality in AMReX makes it easy to efficiently solve for the implicit terms in IMEX methods and is straightforward to use with both single-level and multi-level grids. The code for the first model is implemented in the pure Fortran90 BoxLib framework. For testing IMEX-MRI-GARK methods, the second model retains the complex chemistry from the first problem while taking advantage of the most updated version of AMReX.

AMReX uses an MPI+X strategy where grids are distributed over MPI ranks, and X represents a fine-grained parallelization strategy. In recent years, the most common X in the MPI+X strategy has been OpenMP but AMReX now runs effectively on NVIDIA, AMD, and Intel GPUs. The focus in this paper is on pure MPI; future efforts will be able to easily adopt the multicore/GPU strategies. Load balancing over MPI ranks is done on a level-by-level basis using either options of a Morton-order space filling curve or knapsack algorithm. AMReX efficiently manages and caches communication patterns between grids at the same level, grids at different levels, as well as data structures at the same level with different grid layouts. The latter case is useful for when certain physical processes have widely varying computational requirements in different spatial regions, and grid generation with load balancing in mind may be weighted by this metric, rather than cell count.

The MRI methods evaluated in this work are provided by the SUNDIALS suite of time integrators and nonlinear solvers (Hindmarsh et al., 2005; Gardner et al., 2022). The suite is comprised of the ARKODE package of one-step methods for ODEs, the multistep integrators CVODE and IDA for ODEs and DAEs, respectively, the forward and adjoint sensitivity analysis enabled variants, CVODES and IDAS, and the nonlinear algebraic equation solver KINSOL. The ARKODE package provides an infrastructure for developing one-step methods for ODEs and includes time integration modules for stiff, non-stiff, mixed stiff/non-stiff, and multirate systems (Reynolds et al., 2022). ARKODE’s MRIStep module implements explicit MIS methods and high-order explicit and solve-decoupled implicit and IMEX MRI-GARK methods. In this work, we solve each fast time scale ODE (3) using ARKODE’s ARKStep module which provides explicit, implicit, and ARK-IMEX Runge-Kutta methods, though a user may also provide their own fast scale integrator. The specific MRI methods evaluated in this work are built into ARKODE but users may alternatively supply their own coupling tables for MRI methods or Butcher tables for Runge-Kutta methods. At a minimum, to utilize the methods in ARKODE, an application must provide routines to evaluate the ODE right-hand side functions, and ARKODE must be able to perform vector operations on application data.

SUNDIALS packages utilize an object-oriented design and are written in terms of generic mathematical operations defined by abstract vector, matrix, and solver classes. Thus, the packages are agnostic to an application’s particular choice of data structures/layout, parallelization, algebraic solvers, etc. as the specific details for carrying out the required operations are encapsulated within the corresponding classes. Several vector, matrix, and solver modules are provided with SUNDIALS, but users can also provide their own class implementations tailored to their use case. In this work, we supplied custom vector implementations to interface with the native BoxLib and AMReX data structures (one for BoxLib in Fortran, one each for the single and multi-level grid cases of AMReX in C++). Often in-built functionality from BoxLib or AMReX could be used to aid these routines, but sometimes the cells in the grid needed to be iterated over and the operation computed per-cell explicitly. The following is the listing for the vector multiplication operation (multiplies per-element like Matlab’s .* operation) that leverages in-built AMReX operations in the multi-level case:

With explicit MRI methods, as used in the first problem (5), a vector implementing the necessary operations and routines for evaluating the right-hand side functions are all that is required for interfacing with SUNDIALS. For the IMEX-MRI methods, as used in the second problem (9), a nonlinear solver and a linear solver are needed to solve the nonlinear systems that arise at each implicit slow stage. Here we use the native SUNDIALS Newton solver and matrix-free GMRES implementation to solve these algebraic systems. As both solvers are written entirely in terms of vector operations no additional interfacing is required. However, to accelerate convergence of the GMRES iteration an additional wrapper function was created to leverage the AMReX multigrid linear solver as a preconditioner. This solver is based on the “MLABecLaplacian” linear operator,

where $\alpha$ and $\beta$ are spatially-varying coefficients.

The experiment in this section is close in configuration to the three-dimensional flame problem from Section 5.2.2 in Emmett et al. (2014). The problem domain is a cubic box with 1.6 cm on each side. The boundaries are periodic in the y and z directions and inlet/outlet in the x-direction. The chemistry mechanism for the flame is the GRI-MECH 3.0 network (Frenklach et al., 1995) that contains 53 species and 325 reactions. Combined with the other state variables such as density and velocity, the resulting system has 58 unknowns at each spatial location. The initial conditions are taken from the PREMIX code from the CHEMKIN library (Kee et al., 1998), where we use a precomputed realistic distribution of species in a premixed methane flame in 1D along the x-axis with offsets in the pattern varying in the x-y plane by perturbing spherical harmonics, resulting in a bumpy pattern at the front. The problem is integrated over a total time of $t=2.0\times 10^{-6}$ s using coarse time step sizes that range from $H=2.0\times 10^{-6}$ to $H=1.5625\times 10^{-8}$ s.

The domain interior is discretized with $32^{3}$ points. Since there are $58$ unknowns per spatial location, this discretization gives a total of 1,900,544 unknowns over the interior of the domain. The problem is parallelized by splitting each dimension in half, resulting in 8 boxes distributed over 8 MPI ranks. Since the discretization also requires an additional four layers of ghost cells at each face, the total number of data in the problem becomes 6,414,336 values when including ghost cells.

As with Emmett et al. (2014), we test this problem on the single-rate SDC and the two multirate SDC methods MRSDC-338 and MRSDC-352. Here we also compare against MRI3. As noted in Section 3.2, the MRI3 scheme is paired with fast methods of orders two through five. For each method we perform k fast substeps per slow step (where that ratio is held constant over the problem), and denote the resulting combinations as MRI-P(k), where is P the order of the fast method. Through exhaustive testing, on this application we found the following options to be the most efficient for each fast method order: MRI-2(8), MRI-3(4), MRI-4(8), and MRI-5(4).

The results were gathered on a node of the Quartz cluster at Lawrence Livermore National Laboratory. Each node is composed of two Intel Xeon E5-2695 v4 chips. The tests used 8 ranks of MPI parallelism distributed over the cores of the node.

The range of magnitudes of the different solution components (e.g. density, velocity, mass fractions of different species) is quite large in this problem, so a norm of error over the solution would not properly capture the total error over all components without using an appropriate weighting vector in the norm. As in Emmett et al. (2014), we find the order of convergence of the methods to not depend on which solution component is chosen. Therefore, rather than use a norm over all components, we measure the error as the max norm in density. The results are similar for other components.

The accuracy of the methods is tested against a reference solution calculated using a six-stage fourth-order explicit Runge-Kutta method (Kennedy et al., 2000) with a step size of $5.0\times 10^{-10}$ s. The integration of each method is done using constant time stepping, starting each with the coarsest time step that exhibited numerical stability and halving the step for each successive run. The coarsest allowable time step for each method can be seen in Table 1. Not included in the table is the single-rate SDC case, which required a step size no greater than $6.25\times 10^{-8}$ s to be stable. We note that the single-rate restriction of $\mathcal{O}(10^{-8})$ is consistent with the analysis of the eigenvalues in the Jacobian of this reaction network (Wartha et al., 2020).

The rate of convergence for each of the methods is shown in Figure 3 and the same data is shown in tabular form in Table 1. We note that the convergence rates for all cases remains very similar if we consider the $L_{2}$ norm instead of $L_{\infty}$ norm. In general, the MRSDC methods exhibit highly non-monotonic convergence rates, despite their theoretical fourth-order accuracy. This problem could possibly be alleviated by further reducing the step size, or increasing the number of SDC iterations, as has been demonstrated in other works using SDC for reacting flow (Nonaka et al., 2012, 2018; Pazner et al., 2016). The rate of convergence of the MRI methods is a balance between the accuracy of the fast and slow partitions. The MRI-4 and MRI-5 methods are asymptotically limited to third-order accuracy due to the third-order accuracy of the slow integrator, however the slope of the convergence plot is better than four in some ranges when the reactions dominate the overall accuracy. The MRI-2 method is asymptotically limited to only second-order as step sizes become small, but the slope of the line is at least three for the ranges of step sizes used in the tests. Overall, the MRI methods exhibit convergence rates consistent with the expected order of accuracy, and appear to reach the asymptotic range of convergence at larger step sizes than the MRSDC approaches.

Figure 4 shows an efficiency diagram comparing the performance of the various integrators, plotting the error in density as a function of wallclock time. The figure confirms the results of Emmett et al. (2014) that the MRSDC methods have a significant performance advantage over the single-rate SDC method on this problem. We also see that in general, the MRI methods produce smaller error for a given wallclock time compared to the SDC approaches. The exception is for the smallest step sizes reported, the MRSDC352 4-iteration performs comparably to the MRI approaches. As discussed in Section 2.2, the performance benefits of MRI approaches arise because the MRI methods need fewer calls to the right-hand-side functions per step than the SDC methods. To further quantify this point, Table 2 presents the $H=6.25\times 10^{-8}$ s data from Figure 4 / Table 1, including the number of function calls, runtime, and error using four different approaches. The error for the MRI approaches is comparable to or smaller than the SDC approaches, in addition to offering runtimes that are notably smaller due to the reduced number of slow and fast evaluations.

The problem domain is a 3D cube $[-0.05,0.05]$ cm on a side with periodic boundaries in each direction. The chemistry mechanism is the same GRI-MECH 3.0 network (Frenklach et al., 1995) as the previous problem (53 species, 325 reactions). The model is integrated over a total time of $t=2.0\times 10^{-6}$ s.

The problem is initialized using the same predetermined realistic distribution of species in 1D for a fuel mix as that used in the previous example. However, here the 1D distribution is swept angularly around the origin with high concentration offset from the middle, resulting in a ring of fuel centered around the center. The temperature is set as $1500$ K everywhere, the pressure is set at $1$ atmosphere everywhere, and the density is calculated from the standard equation of state from the pressure and temperature at each cell. The velocity of the system is $\langle 0.001,0.001,0.001\rangle$ cm/s uniformly across the domain.

The problem is tested using ARK-IMEX, IMEX-MRI3 and IMEX-MRI4. Details of the methods and the partitioning of the model terms are described in Section 3.2. The time stepping is done with constant step sizes, with ten fast steps taken for each slow step in the multirate methods. The multirate methods are thus labeled as “IMEX-MRI3-10” and “IMEX-MRI4-10” in what follows.

For the implicit partitions, the implicit stage solves are done using SUNDIALS’ Newton iteration and matrix-free GMRES solver. No preconditioning is used. To prevent iterations being terminated before reaching tolerance, the maximum number of nonlinear and linear iterations is relaxed from the default to 10 and 100 respectively. However, solver tolerances also need to be set for the nonlinear and linear solvers in order to avoid “over-solving” the algebraic systems of equations as well. The nonlinear solver tolerances were chosen on a per-method and per-timestep size basis. The chosen solver tolerance values are listed in Tables 3 and 4. The reported values were chosen by sweeping each case over a range of tolerances from $10^{-1}$ to $10^{-15}$, with jumps of a factor of $10$ between each, with the reported value being the coarsest one for which both the H2 and H mass fraction error is perturbed by less than a factor of $0.01$. The tolerance value is used as an absolute tolerance in the WRMS norm for the nonlinear solver and the 2-norm for the linear solver over the entire solution vector, i.e. over all species at once. Safety factors of 0.1 and 0.05 are applied for the nonlinear and linear solvers, respectively.

The reference solution is computed using the fifth-order explicit Runge–Kutta method from Cash and Karp (1990) with a step size of $10^{-10}$, which is an order of magnitude smaller than the smallest fast step size used in all the runs.

The grid resolution is 128 cells on a side. Diffusion coefficients are computed in the same manner and library code as the SMC problem, using temperatures, molar fractions, and specific heat at constant pressure. The calculation of the diffusion coefficient is computationally expensive, and dominates the cost of the implicit terms more than the nonlinear solves for the tested configuration of the problem. The coefficients are not updated within the solves, so the computational cost of that portion is not affected by the changes in iteration counts with stiffness. However, in order to make the problem sufficiently stiff to warrant an implicit method, the computed coefficient in each cell is scaled by a factor of 1000. The problem is thus not physically realistic, but is meant to exercise the same code components that are expensive in a reactive-flow problem at a level of stiffness that would be found in highly resolved realistic cases. The scale of the diffusion coefficient combined with the grid resolution results in a problem that cannot be completed successfully using explicitly evolved diffusion at the coarser step sizes of the runs.

The problem is run with 512 total ranks over 16 nodes on the Quartz cluster at Lawrence Livermore National Laboratory.

The results of the runs are shown in Figure 5. The error is shown for the species H2, but as with the SMC problem, the outcome is insensitive to which species is used to measure error.

We see in the plots that there is a significant performance advantage when using multirate over single-rate IMEX methods for this problem. We tried the problem over a range of fast-to-slow step size ratios and found ten-to-one to be the best balance. Due to the expensive diffusion coefficient calculation, reducing the number of coarse steps that must be computed per fast step results in a large time savings.

An interesting phenomenon from Tables 3 and 4 is that the solver tolerances needed by the multirate methods were more stringent than for the single-rate method for the same coarse step size. This seems to be due to the increased accuracy of the MRI solutions at those step sizes, thereby necessitating tighter algebraic solver tolerances.

Table 5 shows the scaling in CPU time and change in nonlinear and linear iterations for the problem as the grid resolution is refined. The results are from using the IMEX-MRI3-10 method over eight time steps, i.e. at about mid-height in Figure 5. Even at a somewhat finer step size than the coarsest case, the problem is moderately stiff so the number of linear iterations grows with the grid resolution. However, the primary cost is the computation of the diffusion coefficients, which are not done within the linear iterations. As a result, the total cost goes up only slowly with stiffness within the range of resolutions tested. If the problem were to continue to grow in size, eventually the solve cost would take over the cost of the diffusion coefficient calculation and a preconditioner would be needed. As mentioned in Section 4.2, AMReX includes the MLABecLaplacian multigrid solver that works over multiple AMR grid levels. Here, we use this preconditioner for this single-level problem, with results comparing the un-preconditioned and preconditioned solvers in Tables 5 and 6, respectively. Due to the overhead in invoking the preconditioner, for coarse problems the un-preconditioned solve is more efficient, but once the grid reaches at least 128 cells per side the preconditioner becomes cost effective. Tables 7 and 8 push this experiment further, by artificially increasing the problem stiffness through scaling the diffusion coefficients by another factor of 10. Here again, preconditioning is effective at grid sizes exceeding 64 cells per side, and becomes critical as the grids are increased to 256 cells per side.

As a proof of concept and demonstration of capability, problem (9) is also implemented in AMReX on multi-level grids. The problem now uses a ball of pre-mixed methane fuel in the center of the domain, but is otherwise the same as from the single-level case. Figure 6 shows a 2D slice of the methane mass fraction at the initial time point, where we have overlaid the coarse and fine grid structures. Unlike what is typically done with AMR, the multirate time stepping here is dictated by the accuracy requirements of the chemistry and diffusion and not by the CFL stability restriction. The advection is computed explicitly, but the velocity is too small in this problem for CFL limitations to restrict step sizes. The diffusion is computed implicitly and thus also does not limit the step size. Therefore, the coarse step sizes for diffusion and advection and fine time step size for reactions are respectively the same for all levels of the multi-level grid.

To demonstrate the advantage of using both multirate time stepping and multi-level grids, versus without or with either option alone, we ran both ARK-IMEX and IMEX-MRI4-10 on single and two-level grids. In the multilevel case, the fine grid only covers the fuel in the center, i.e., for a radius of $0.01$ cm extending from the center of the domain (which covers $[-0.05,0.05]$ per side). In the single level case, the domain contains $256^{3}$ grid cells divided into blocks of 32 cells on a side. In the multilevel case, we use an effective $256^{3}$ resolution with a coarse domain containing $128^{3}$ grid cells, using blocks of 16 cells on a side to maintain the same number of MPI ranks as the single level case, along with another level of refinement. A multilevel setup of our simulation is shown in Figure 6; we note that the fine grids cover roughly 5% of the domain volume. In the figure, the resolution is lower for clarity. We leverage the composite multilevel MLABecLaplacian operator built into AMReX for implicit diffusion. The solver tolerance is $10^{-12}$ for all cases.

Table 9 shows that using either AMR or multirate alone takes considerably less time than single-rate time stepping on the single level grid. However, combining the two techniques results in a considerably more efficient result than with either feature alone.