High-Order Algorithms for Compressible Reacting Flow with Complex Chemistry

The multicomponent reacting compressible Navier-Stokes equations are given by

where $\rho$ is the density, $u$ is the velocity, $p$ is the pressure, $E$ is the specific energy density (kinetic energy plus internal energy), $T$ is the temperature, $\tau$ is the viscous stress tensor, and $\lambda$ is the thermal conductivity. For each of the chemical species $k$, $Y_{k}$ is the mass fraction, $\mathcal{F}_{k}$ is the diffusive flux, $h_{k}$ is the specific enthalpy, and $\dot{\omega}_{k}$ is the production rate. The system is closed by an equation of state that specifies $p$ as a function of $\rho$, $T$ and $Y_{k}$. For the examples presented here, we assume an ideal gas mixture.

The viscous stress tensor is given by

where $\eta$ is the shear viscosity and $\xi$ is the bulk viscosity.

A mixture model for species diffusion is employed in SMC. We define an initial approximation $\bar{\mathcal{F}}_{k}$ to the species diffusion flux given by

where $X_{k}$ is the mole fraction, and $D_{k}$ is the diffusion coefficient of species $k$. For species transport to be consistent with conservation of mass we require

In general, the fluxes defined by (3) do not satisfy this requirement and some type of correction is required. When there is a good choice for a reference species that has a high molar concentration throughout the entire domain, one can use (3) to define the species diffusion flux (i.e., $\mathcal{F}_{k}\equiv\bar{\mathcal{F}}_{k}$) for all species except the reference, and then use (4) to define the flux of the reference species. This strategy is often used in combustion using $\mathrm{N_{2}}$ as the reference species. An alternative approach [30] is to define a correction velocity $\mbox{\boldmath$V$}_{c}$ given by

and then to set

The SMC code supports both of these operating modes. We will discuss these approaches further in §2.4.

The SMC code discretizes the governing Navier-Stokes equations (1) in space using high-order centered finite-difference methods. For first order derivatives a standard $8^{\rm th}$ order stencil is employed. For second order derivatives with variable coefficients of the form

a novel $8^{\rm th}$ order narrow stencil is employed. This stencil is an extension of those developed by Kamakoti and Pantano [18]. The derivatives in (7) on a uniform grid are approximated by

where

$s=4$, and $M$ is an $8\times 8$ matrix (in general, for a discretization of order $2s$, $M$ is a $2s\times 2s$ matrix). This discretization is conservative because (8) is in flux form. Note, however, that $H_{i+1/2}$ is not a high-order approximation of the flux at $x_{i+1/2}$ even though the flux difference (8) approximates the derivative of the flux to order $2s$. This type of stencil is narrow because the width of the stencil is only $2s+1$ for both $u$ and $a$. In contrast, the width of a wide stencil in which the first-order derivative operator is applied twice is $4s+1$ for $u$ and $2s+1$ for $a$.

Following the strategy in [18], a family of $8^{\rm th}$ order narrow stencils for the second order derivative in (7) has been derived. The resulting $8^{\rm th}$ order stencil matrix is given by

where

and $m_{47}$ and $m_{48}$ are two free parameters. Note that the well-posedness requirement established in [18] is satisfied by the stencil matrix above independent of the two free parameters. In SMC the default values of the free parameters were set to $m_{47}=3557/44100$ and $m_{48}=-2083/117600$ to minimize an upper bound on truncation error that is calculated by summing the absolute value of each term in the truncation error. Additional analysis of truncation errors and numerical experiments on the $8^{\rm th}$ order narrow stencil using various choices of the two parameters are presented in Appendix 7.

The SMC code uses a $6^{\rm th}$ order stencil near physical boundaries where the conditions outside the domain are not well specified (see §2.5). The general $6^{\rm th}$ order stencil matrix is given by

where

and $m_{36}$ is a free parameter. The $6^{\rm th}$ order narrow stencil developed in [18] is a special case of the general $6^{\rm th}$ order stencil when $m_{36}=0.220063$. Like the $8^{\rm th}$ order stencil, the well-posedness and the order of accuracy for the $6^{\rm th}$ order stencil is independent of the free parameter. In SMC, the free parameter is set to $281/3600$ to minimize an upper bound on truncation error that is calculated by summing the absolute value of each term in the truncation error.

The $8^{\rm th}$ order narrow stencil presented in § 2.2 requires significantly more floating-point operations (FLOPs) than the standard wide stencil, but has the benefit of less communication. The computational cost of a single evaluation the narrow stencil, for $N$ cells in one dimension, is $112N+111$ FLOPs (this is the cost of operations in (9) and (8), excluding the division by $\Delta x^{2}$). For the $8^{\rm th}$ order wide stencil with two steps of communication, the computational cost on $N$ cells is $23N+8$ FLOPs. Hence, the narrow stencil costs $89N+103$ more FLOPs than the wide stencil on $N$ cells in one dimension. As for communication, the wide stencil requires $64$ more bytes than the narrow stencil for $N$ cells with 8 ghost cells in one dimension. Here we assume that an 8-byte double-precision floating-point format is used. Ignoring latency, the difference in time for the two algorithms is then

where $F$ is the floating-point operations per second (FLOPS) and $B$ is the network bandwidth. To estimate this difference on current machines, we used the National Energy Research Scientific Computing Center’s newest supercomputer Edison as a reference. Each compute node on Edison has a peak performance of 460.8 GFLOPs (Giga-FLOPs per second), and the MPI bandwidth is about 8 GB/s ¹¹1The configuration of Edison is available at http://www.nersc.gov/users/computational-systems/edison/configuration/.. For $N=64$ in one dimension, the run time cost for the extra FLOPs of the narrow stencil is about $1.3\times 10^{-8}\,\mathrm{s}$, which is about 1.6 times of the extra communication cost of the wide stencil, $8.0\times 10^{-9}\,\mathrm{s}$. It should be emphasized that the estimates here are crude for several reasons. We neglected the efficiency of both floating-point and network operations for simplicity. In counting FLOPs, we did not consider a fused multiply-add instruction or memory bandwidth. In estimating communication costs, we did not consider network latency, the cost of data movement if MPI messages are aggregated to reduce latency, or communication hiding through overlapping communication and computation. Moreover, the efficiency of network operations is likely to drop when the number of nodes used increases, whereas the efficiency of FLOPs alone is independent of the number of total nodes. Nevertheless, as multicore systems continue to evolve, the ratio of FLOPS to network bandwidth is most likely to increase. Thus, we expect that the narrow stencil will cost less than the wide stencil on future systems, although any net benefit will depend on the systems configuration. As a data point, we note that both stencils are implemented in SMC, and for a test with $512^{3}$ cells using 512 MPI ranks and 12 threads on 6144 cores of Edison, a run using the narrow stencil took about the same time as using the wide stencil. Further measurements are presented in §5.2.

In §2.1, we introduced the two approaches employed in SMC to enforce the consistency of the mixture model of species diffusion with respect to mass conservation. The first approach uses a reference species and is straightforward to implement. We note that the last term in (1d) can be treated as $-\nabla\cdot\sum_{k}\mathcal{F}_{k}(h_{k}-h_{\mathrm{ref}})$, where $h_{\mathrm{ref}}$ is the species enthalpy for the reference species, because of (4). The second approach uses a correction velocity (5) and involves the computation of $\nabla\cdot(Y_{k}\mbox{\boldmath$V$}_{c})$ and $\nabla\cdot(h_{k}Y_{k}\mbox{\boldmath$V$}_{c})$. Although one might attempt to compute the correction velocity according to $V_{c,i+1/2}=\sum_{\ell}H_{\ell,i+1/2}$, where $H_{\ell,i+1/2}$ for species $\ell$ is computed using (9), this would reduce the accuracy for correction velocity to second-order since $H_{i+1/2}$ is only a second-order approximation to the flux at $i+1/2$. To maintain $8^{\rm th}$ order accuracy without incurring extra communication and too much extra computation, we employ the chain rule for the correction terms. For example, we compute $\nabla\cdot(h_{k}Y_{k}\mbox{\boldmath$V$}_{c})$ as $h_{k}Y_{k}\nabla\cdot\mbox{\boldmath$V$}_{c}+\nabla(h_{k}Y_{k})\cdot\mbox{\boldmath$V$}_{c}$, and use $\nabla\cdot\mbox{\boldmath$V$}_{c}=\sum_{\ell}\nabla\cdot\bar{\mathcal{F}}_{\ell}$, where $\nabla\cdot\bar{\mathcal{F}}_{\ell}$ is computed using the $8^{\rm th}$ order narrow stencil for second order derivatives. The velocity correction form arises naturally in multicomponent transport theory; namely, it represents the first term in an expansion of the full multicomponent diffusion matrix. (See Giovangigli [13].) For that reason, velocity correction is the preferred option in SMC.

The SMC code uses characteristic boundary conditions to treat physical boundaries. Previous studies in the literature have shown that characteristic boundary conditions can successfully lower acoustic reflections at open boundaries without suffering from numerical instabilities. Thompson [36, 37] applied the one-dimensional approximation of the characteristic boundary conditions to the hyperbolic Euler equations, and the one-dimensional formulation was extended by Poinsot and Lele [29] to the viscous Navier-Stokes equations. For multicomponent reacting flows, Sutherland and Kennedy [34] further improved the treatment of boundary conditions by the inclusion of reactive source terms in the formulation. In SMC, we have adopted the formulation of characteristic boundary conditions proposed in [41, 39], in which multi-dimensional effects are also included.

The $8^{\rm{th}}$ order stencils used in SMC require four cells on each side of the point where the derivative is being calculated. Near the physical boundaries where the conditions outside the domain are unknown, lower order stencils are used for derivatives with respect to the normal direction. At the fourth and third cells from the boundary we use $6^{\rm{th}}$ and $4^{\rm{th}}$ order centered stencils respectively, for both advection and diffusion. At cells next to the boundary we reduce to centered second-order for diffusion and a biased third-order discretization for advection. Finally, for the boundary cells, a completely biased second-order stencil is used for diffusion and a completely biased third-order stencil is used for advection.

SDC methods for ODEs were first introduced in [10] and have been subsequently refined and extended in [25, 26, 16, 14, 5, 22]. SDC methods construct high-order solutions within one time step by iteratively approximating a series of correction equations at collocation nodes using low-order sub-stepping methods.

Consider the generic ODE initial value problem

where $t\in[0,T]$; $u_{0},\,u(t)\in{\mathbb{C}}^{N}$; and $f:{\mathbb{R}}\times{\mathbb{C}}^{N}\rightarrow{\mathbb{C}}^{N}$. SDC methods are derived by considering the equivalent Picard integral form of (12) given by

A single time step $[t_{n},t_{n+1}]$ is divided into a set of intermediate sub-steps by defining $M+1$ collocation points ${\tau}_{m}\in[t_{n},t_{n+1}]$ such that $t_{n}={\tau}_{0}<{\tau}_{1}<\cdots<{\tau}_{M}=t_{n+1}$. Then, the integrals of $f\bigl{(}u(t),t\bigr{)}$ over each of the intervals $[t_{n},{\tau}_{m}]$ are approximated by

where $U_{j}\approx u({\tau}_{j})$, $\Delta t=t_{n+1}-t_{n}$, and $q_{mj}$ are quadrature weights. The quadrature weights that give the highest order of accuracy given the collocation points ${\tau}_{m}$ are obtained by computing exact integrals of the Lagrange interpolating polynomial over the collocation points ${\tau}_{m}$. Note that each of the integral approximations $I_{n}^{m}$ (that represent the integral of $f$ from $t_{n}$ to ${\tau}_{m}$) in (14) depends on the function values $f(U_{m},{\tau}_{m})$ at all of the collocation nodes ${\tau}_{m}$.

To simplify notation, we define the integration matrix ${\bm{Q}}$ to be the $M\times(M+1)$ matrix consisting of entries $q_{mj}$; and the vectors ${\bm{U}}\equiv[U_{1},\cdots,U_{M}]$ and ${\bm{F}}\equiv[f(U_{0},t_{0}),\cdots,f(U_{M},t_{M})]$. With these definitions, the Picard equation (13) within the time step $[t_{n},t_{n+1}]$ is approximated by

where ${\bm{U}}_{0}=U_{0}\bm{1}$. Note again that the integration matrix ${\bm{Q}}$ is dense so that each entry of ${\bm{U}}$ depends on all other entries of ${\bm{U}}$ (through ${\bm{F}}$) and $U_{0}$. Thus, (15) is an implicit equation for the unknowns in ${\bm{U}}$ at all of the quadrature points. Finally, we note that the solution of (15) corresponds to the collocation or implicit Runge-Kutta solution of (12) over the nodes ${\tau}_{m}$. Hence SDC can be considered as an iterative method for solving the spectral collocation formulation.

The SDC scheme used here begins by spreading the initial condition $u_{n}$ to each of the collocation nodes so that the provisional solution ${\bm{U}}^{0}$ is given by ${\bm{U}}^{0}=[U_{0},\cdots,U_{0}]$. Subsequent iterations (denoted by $k$ superscripts) proceed by applying the node-to-node integration matrix ${\bm{S}}$ to ${\bm{F}}^{k}$ and correcting the result using a forward-Euler time-stepper between the collocation nodes. The node-to-node integration matrix ${\bm{S}}$ is used to approximate the integrals $I_{m}^{m+1}=\int_{{\tau}_{m}}^{{\tau}_{m+1}}f(u(s),s)\,ds$ (as opposed to the integration matrix ${\bm{Q}}$ which approximates $I_{n}^{m}$) and is constructed in a manner similar to the integration matrix ${\bm{Q}}$. The update equation corresponding to the forward-Euler sub-stepping method for computing ${\bm{U}}^{k+1}$ is given by

where $S^{k,m}$ is the $m^{\rm th}$ row of ${\bm{S}}{\bm{F}}^{k}$. The process of computing (16) at all of the collocation nodes ${\tau}_{m}$ is referred to as an SDC sweep or an SDC iteration. The accuracy of the solution generated after $k$ SDC iterations done with such a first-order method is formally $O(\Delta t^{k})$ as long as the spectral integration rule (which is determined by the choice of collocation nodes ${\tau}_{m}$) is at least order $k$.

Here, a method of lines discretization based on the spatial discretization presented in §2 is used to reduce the partial differential equations (PDEs) in question to a large system of ODEs.

The computational cost of the SDC scheme used here is determined by the number of nodes $M+1$ chosen and the number of SDC iterations $K$ taken per time step. The number of resulting function evaluations $N_{\rm evals}$ is $N_{\rm evals}=KM$. For example, with 3 Gauss-Lobatto ($M=2$) nodes, the resulting integration rule is $4^{\rm th}$ order. To achieve this formal order of accuracy we perform $K=4$ SDC iterations per time step and hence we require $N_{\rm evals}=8$ function evaluations (note that $F(U_{0},0)$ is recycled from the previous time step, and hence the first time-step requires 9 function evaluations).

By increasing the number of SDC nodes and adjusting the number of SDC iterations taken, we can construct schemes of arbitrary order. For example, to construct a $6^{\rm th}$ order scheme one could use a 5 point Gauss-Lobatto integration rule (which is $8^{\rm th}$ order accurate), but only take 6 SDC iterations. For high order explicit schemes used here we have observed that the SDC iterations sometimes converge to the collocation solution in fewer iterations than is formally required. This can be detected by computing the SDC residual

where $\bm{q}$ is the last row of ${\bm{Q}}$, and comparing successive residuals. If $|R^{k-1}|/|R^{k}|$ is close to one, the SDC iterations have converged to the collocation solution and subsequent iterations can be skipped.

Multirate SDC (MRSDC) methods use a hierarchy of SDC schemes to integrate systems that contain processes with disparate time-scales more efficiently than fully coupled methods [5, 22]. A traditional MRSDC method integrates processes with long time-scales with fewer SDC nodes than those with short time-scales. They are similar in spirit to sub-cycling schemes in which processes with short characteristic time-scales are integrated with smaller time steps than those with long characteristic time-scales. MRSDC schemes, however, update the short time-scale components of the solutions at the long time-scale time steps during each MRSDC iteration, as opposed to only once for typical sub-cycling schemes, which results in tighter coupling between processes throughout each time-step [5, 22]. The MRSDC schemes presented here generalize the structure and coupling of the quadrature node hierarchy relative to the schemes first presented in [5, 22].

Before presenting MRSDC in detail, we note that the primary advantage of MRSDC schemes lies in how the processes are coupled. For systems in which the short time-scale process is moderately stiff, using more SDC nodes than the long time-scale process helps alleviate the stiffness of the short time-scale processes. However, extremely stiff processes are best integrated between SDC nodes with, for example, a variable order/step-size backward differentiation formula (BDF) scheme. In this case, the short time-scale processes can in fact be treated with fewer SDC nodes than long time-scale processes.

Consider the generic ODE initial value problem with two disparate time-scales

where we have split the generic right-hand-side $f$ in (12) into two components: $f_{1}$ and $f_{2}$. A single time step $[t_{n},t_{n+1}]$ is first divided into a set of intermediate sub-steps by defining $M_{1}+1$ collocation points such that ${\tau}_{m}\in[t_{n},t_{n+1}]$. The vector of these $M_{1}+1$ collocation points ${\tau}_{m}$ is denoted $\bm{t}_{1}$, and these collocation points will be used to integrate the $f_{1}$ component. The time step is further divided by defining $M_{2}+1$ collocation points ${\tau}_{p}\in[t_{n},t_{n+1}]$ such that $\bm{t}_{1}\subset\bm{t}_{2}$, where $\bm{t}_{2}$ is the vector of these new ${\tau}_{p}$ collocation points. These collocation points will be used to integrate the $f_{2}$ component. That is, a hierarchy of collocation points is created with $\bm{t}_{1}\subset\cdots\subset\bm{t}_{N}$. The splitting of the right-hand-side in (18) should be informed by the numerics at hand. Traditionally, slow time-scale processes are split into $f_{1}$ and placed on the $\bm{t}_{1}$ nodes, and fast time-scale processes are split into $f_{2}$ and places on the $\bm{t}_{2}$ nodes. However, if the fast time-scale processes in the MRSDC update equation (20) are extremely stiff then it is more computationally effective to advance these processes using stiff integration algorithms over a longer time scale. In this case, we can swap the role of of the processes within MRSDC and treat the explicit update on the faster time scale (i.e., with the $\bm{t}_{2}$ nodes). As such, we hereafter refer to the collocation points in $\bm{t}_{1}$ as the “coarse nodes” and those in $\bm{t}_{2}$ as the “fine nodes”.

The SDC nodes $\bm{t}_{1}$ for the coarse component are typically chosen to correspond to a formal quadrature rule like the Gauss-Lobatto rule. These nodes determine the overall order of accuracy of the resulting MRSDC scheme. Subsequently, there is some flexibility in how the nodes for fine components are chosen. For example, the points in $\bm{t}_{2}$ can be chosen to correspond to a formal quadrature rule

1. (a)

   over the entire time step $[t_{n},t_{n+1}]$;

2. (b)

   over each of the intervals defined by the points in $\bm{t}_{1}$;

3. (c)

   over several sub-intervals of the intervals defined by the points in $\bm{t}_{1}$;

Note that in the first case, a quadrature rule that nests (or embeds) itself naturally is required (such as the Clenshaw-Curtis rule), whereas in the second and third cases the nodes in $\bm{t}_{2}$ usually do not correspond to a formal quadrature rule over the entire time step $[t_{n},t_{n+1}]$. In the second and third cases we note that the basis used to construct the ${\bm{Q}}$ and ${\bm{S}}$ integration matrices are typically piecewise polynomial – this allows us to shrink the effective time-step size of the fine process without constructing an excessively high-order polynomial interpolant. These three possibilities are depicted in Fig. 1.

Once the SDC nodes for the desired MRSDC scheme have been chosen, the node-to-node and regular integration matrices (${\bm{S}}$ and ${\bm{Q}}$, respectively) can be constructed. While single-rate SDC schemes, as presented in §3.1, require one set of ${\bm{S}}$ and ${\bm{Q}}$ matrices, the two-component MRSDC scheme used here requires two sets of ${\bm{S}}$ and ${\bm{Q}}$ matrices. These are denoted with subscripts $i,j$ so that the action of ${\bm{S}}_{i,j}$ approximates the node-to-node integrals between nodes in $\bm{t}_{i}$ of the ${\bm{F}}_{j}$ component that “lives” on the $\bm{t}_{j}$ nodes, and similarly for the ${\bm{Q}}_{i,j}$ matrices. Note that both ${\bm{S}}_{i,j}$ and ${\bm{Q}}_{i,j}$ are of size $M_{i}\times(M_{j}+1)$.

With these definitions, the Picard equation (13) within the time step $[t_{n},t_{n+1}]$ is approximated by

where ${\bm{U}}$ is the vector of approximate solution values $U_{q}$ over the fine nodes in $\bm{t}_{2}$. The MRSDC scheme used in SMC has two operating modes:

1. 1.

   The coarse ($f_{1}$) and fine ($f_{2}$) components of (18) correspond to the advection-diffusion and reaction parts of (1), respectively, and the forward-Euler sub-stepping scheme is used for both components.

2. 2.

   The coarse ($f_{1}$) and fine ($f_{2}$) components of (18) correspond to the reaction and advection-diffusion parts of (1), respectively. The advection-diffusion component is treated with a forward-Euler sub-stepper while the reaction component is treated with the variable-order BDF scheme in VODE [6].

In the first mode, the scheme is purely explicit whereas in the second mode the algorithm uses an implicit treatment of reactions. Which version is most appropriate for a given simulation depends on the relative stiffness of chemistry and fluid dynamics. When the chemical time scales are not significantly disparate from the fluid mechanical time scale, the fully explicit version will be the most efficient. When the chemistry is extremely stiff, the second mode is preferred.

The MRSDC update equation used in SMC for mode 1 over the nodes in $\bm{t}_{2}$ is given by

where $S_{2,j}^{k,q}$ is the $q^{\rm th}$ row of ${\bm{S}}_{2,j}{\bm{F}}_{j}$, and the $p^{\rm th}$ node of the fine component corresponds to the closest coarse node to the left of the fine node $q$. That is, the coarse component $f_{1}$ is held constant, and equal to its value at the left coarse node $p=\lfloor q/(M_{2}/M_{1})\rfloor$, during the fine updating steps. The update equation for mode 2 is

where $f_{1}(U^{k+1}_{p+1},t_{p})$ is computed with VODE.

The computational cost of an MRSDC scheme is determined by the number of nodes $M_{\ell}+1$ chosen for each multirate component $\ell$ and the number of SDC iterations $K$ taken per time step. The number of resulting function evaluations $N_{\rm evals}^{\ell}$ per component $\ell$ is $N_{\rm evals}^{\ell}=KM_{\ell}$. For example, with 3 Gauss-Lobatto ($M=2$) nodes for the coarse component and 5 Gauss-Lobatto nodes for the fine component between each of the coarse nodes (so that there are 9 nodes in $\bm{t}_{2}$ and hence $M_{2}=8$), the number of coarse function evaluations over $K=4$ iterations is $N_{\rm evals}^{1}=8$, and the number of fine function evaluations is $N_{\rm evals}^{2}=32$.

The reacting compressible Navier-Stokes equations (1) have two characteristic time-scales: one over which advection and diffusion processes evolve, and another over which reactions evolve. The reaction time-scale is usually shorter than the advection/diffusion time-scale. As such, the coarse and fine MRSDC nodes typically correspond to the advection/diffusion and reaction parts of (1), respectively. However, as is the case for the Dimethyl Ether Jet problem shown in §5.3, the reaction terms in (1) can become so stiff that it becomes impractical to integrate them explicitly. In these cases, SMC employs the variable-order BDF schemes in VODE [6] to integrate the reaction terms implicitly. This effectively reverses the time-step constraint of the system: the reaction terms can be integrated using VODE (which internally takes many steps) with a larger time-step than the advection-diffusion terms. As such, the roles of the coarse and fine MRSDC nodes can be reversed so that the coarse and fine MRSDC nodes correspond to the reaction and advection-diffusion terms of (1), respectively, thereby amortizing the cost of the implicit solves by using fewer SDC nodes to integrate the stiff chemistry.