The PeleLM Model
In this section, we present the actual model that is evolved numerically by PeleLM, and the numerical algorithms to do it. There are many control parameters to customize the solution strategy and process, and in order to actually set up and run specific problems with PeleLM, the user must specific the chemical model, and provide routines that implement initial and boundary data and refinement criteria for the adaptive mesh refinement. We discuss setup and control of PeleLM in later sections.
Overview of PeleLM
PeleLM evolves chemically reacting low Mach number flows with block-structured adaptive mesh refinement (AMR). The code depends upon the AMReX library to provide the underlying data structures, and tools to manage and operate on them across massively parallel computing architectures. PeleLM also borrows heavily from the source code and algorithmic infrastructure of the IAMR. IAMR implements an AMR integration for the variable-density incompressible Navier-Stokes equations. PeleLM extends IAMR to include complex coupled models for generalized thermodynamic relationships, multi-species transport and chemical reactions. The core algorithms in PeleLM (and IAMR) are described in the following papers:
A conservative, thermodynamically consistent numerical approach for low Mach number combustion. I. Single-level integration, A. Nonaka, J. B. Bell, and M. S. Day, Combust. Theor. Model., 22 (1) 156-184 (2018)
A Deferred Correction Coupling Strategy for Low Mach Number Flow with Complex Chemistry, A. Nonaka, J. B. Bell, M. S. Day, C. Gilet, A. S. Almgren, and M. L. Minion, Combust. Theory and Model, 16 (6) 1053-1088 (2012)
Numerical Simulation of Laminar Reacting Flows with Complex Chemistry, M. S. Day and J. B. Bell, Combust. Theory Model 4 (4) 535-556 (2000)
An Adaptive Projection Method for Unsteady, Low-Mach Number Combustion, R. B. Pember, L. H. Howell, J. B. Bell, P. Colella, W. Y. Crutchfield, W. A. Fiveland, and J. P. Jessee, Comb. Sci. Tech., 140 123-168 (1998)
A Conservative Adaptive Projection Method for the Variable Density Incompressible Navier-Stokes Equations, A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, and M. L. Welcome, J. Comp. Phys., 142 1-46 (1998)
The low Mach number flow equations
PeleLM solves the reacting Navier-Stokes flow equations in the low Mach number regime, where the characteristic fluid velocity is small compared to the sound speed, and the effect of acoustic wave propagation is unimportant to the overall dynamics of the system. Accordingly, acoustic wave propagation can be mathematically removed from the equations of motion, allowing for a numerical time step based on an advective CFL condition, and this leads to an increase in the allowable time step of order \(1/M\) over an explicit, fully compressible method (\(M\) is the Mach number). In this mathematical framework, the total pressure is decomposed into the sum of a spatially constant (ambient) thermodynamic pressure \(P_0\) and a perturbational pressure, \(\pi({\vec x})\) that drives the flow. Under suitable conditions, \(\pi/P_0 = \mathcal{O} (M^2)\).
The set of conservation equations specialized to the low Mach number regime is a system of PDEs with advection, diffusion and reaction (ADR) processes that are constrained to evolve on the manifold of a spatially constant \(P_0\):
where \(\rho\) is the density, \(\boldsymbol{u}\) is the velocity, \(h\) is the mass-weighted enthalpy, \(T\) is temperature and \(Y_m\) is the mass fraction of species \(m\). \(\dot{\omega}_m\) is the molar production rate for species \(m\), the modeling of which will be described later in this section. \(\tau\) is the stress tensor, \(\boldsymbol{\mathcal{Q}}\) is the heat flux and \(\boldsymbol{\mathcal{F}}_m\) are the species diffusion fluxes. These transport fluxes require the evaluation of transport coefficients (e.g., the viscosity \(\mu\), the conductivity \(\lambda\) and the diffusivity matrix \(D\)) which are computed using the library EGLIB, as will be described in more depth in the diffusion section. The momentum source, \(\boldsymbol{F}\), is an external forcing term. For example, we have used \(\boldsymbol{F}\) to implement a long-wavelength time-dependent force to establish and maintain quasi-stationary turbulence.
These evolution equations are supplemented by an equation of state for the thermodynamic pressure. For example, the ideal gas law,
can be used, although PeleLM will soon support other more general expressions, such as Soave-Redlich-Kwong. In the above, \(W_m\) and \(W\) are the species \(m\), and mean molecular weights, respectively. To close the system we also require a relationship between enthalpy, species and temperature. We adopt the definition used in the CHEMKIN standard,
where \(h_m\) is the species \(m\) enthalpy. Note that expressions for \(h_m(T)\) see <section on thermo properties> incorporate the heat of formation for each species.
Neither species diffusion nor reactions redistribute the total mass, hence we have \(\sum_m \boldsymbol{\mathcal{F}}_m = 0\) and \(\sum_m \dot{\omega}_m = 0\). Thus, summing the species equations and using the definition \(\sum_m Y_m = 1\) we obtain the continuity equation:
This, together with the conservation equations form a differential-algebraic equation (DAE) system that describes an evolution subject to a constraint. A standard approach to attacking such a system computationally is to differentiate the constraint until it can be recast as an initial value problem. Following this procedure, we set the thermodynamic pressure constant in the frame of the fluid,
and observe that if the initial conditions satisfy the constraint, an evolution satisfying the above will continue to satisfy the constraint over all time. Expanding this expression via the chain rule and continuity:
The constraint here take the form of a condition on the divergence of the flow. Note that the actual expressions to use here will depend upon the chosen models for evaluating the transport fluxes.
Transport Fluxes
Expressions for the transport fluxes appearing above can be approximated in the Enskog-Chapman expansion as:
where \(\mu\) is the shear viscosity, \(\kappa\) is the bulk viscosity, and \(\lambda'\) is the partial thermal conductivity. In the full matrix diffusion model, the vector of \(m\) species diffusion velocities, \(\boldsymbol{V_m}\), is given by:
where \({D}_{m,j}\) is the diffusion matrix, and \(\boldsymbol{\theta}\) are thermal diffusion coefficients associated with the Soret (mass concentration flux due to an energy gradient) and Dufour (the energy flux due to a mass concentration gradient) effects. The \(m\) species transport driving force due to composition gradients, \(\boldsymbol{d_m}\), is given by:
Alternatively (as in the transport library, EGLIB) the thermal diffusion ratios \(\boldsymbol{\chi}\) may be preferred and the diffusion velocities and energy flux recast as:
where \({D} \boldsymbol{\chi} = \boldsymbol{\theta}\).
As can be seen, the expression for these fluxes relies upon several transport coefficients that need to be evaluated. However, in the present framework several effects are neglected, thus simplifying the fluxes evaluation.
The PeleLM Equation Set
The full diffusion model couples together the advance of all thermodynamics fields, including a dense matrix transport operator that is cumbersome to deal with computationally, while also being generally viewed as an overkill for most practical combustion applications – particularly those involving turbulent fluid dynamics. For PeleLM, we make the following simplifying assumptions:
The bulk viscosity, \(\kappa\), is usually negligible, compared to the shear viscosity,
The low Mach limit implies that there are no spatial gradients in the thermodynamic pressure,
The mixture-averaged diffusion model is assumed,
Dufour and Soret effects are negligible
With these assumptions, the conservation equations take the following form:
with
where \(\boldsymbol{d_m} = \nabla X_m\) and \(\widetilde{D_{m,k}} = Y_m D_{m,k}\). Using these expressions, we can write an equation for \(T\) that is needed in order to evaluate the right-hand side of the divergence constraint:
where \(C_p = \partial h/\partial T\) is the specific heat of the mixture at constant pressure. For an ideal gas, the constraint then becomes:
The mixture-averaged transport coefficients discussed above (\(\mu\), \(\lambda\) and \(D_{m,mix}\)) can be evaluated from transport properties of the pure species. We follow the treatment used in the EGLib library, based on the theory/approximations developed by Ern and Givangigli (however, PeleLM uses a recoded version of these routines that are thread safe and vectorize well on suitable processors).
The following choices are currently implemented in PeleLM
The viscosity, \(\mu\), is estimated based on one step of the conjugate gradient method, using temperature dependent ratios of collisions integrals (EGZE3).
The conductivity, \(\lambda\), is based on an empirical mixture formula (EGZL1):
with
The flux diffusion matrix is approximated using the diagonal of the flux diffusion vector \(\rho \widetilde{ \Upsilon}\), where:
and the \(\mathcal{D}_{m,j}\) are the binary diffusion coefficients of the pair (m,j). This leads to a mixture-averaged approximation that is similar to that of Hirschfelder-Curtiss (EGZVR1):
Note that with these definitions, there is no guarantee that \(\sum \boldsymbol{\mathcal{F}}_{m} = 0\), as required for mass conservation. An arbitrary correction flux, consistent with the mixture-averaged diffusion approximation, is added in PeleLM to enforce conservation.
The pure species and mixture transport properties are evaluated with (thread-safe, vectorized) EGLib functions, which require as input polynomial fits of the logarithm of each quantity versus the logarithm of the temperature.
\(q_m\) represents \(\eta_m\), \(\lambda_m\) or \(D_{m,j}\). These fits are generated as part of a preprocessing step managed by the tool FUEGO based on the formula (and input data) discussed above. The role of FUEGO is to preprocess the model parameters for transport as well as chemical kinetics and thermodynamics.
Chemical kinetics and the reaction source term
Chemistry in combustion systems involves the \(N_s\) species interacting through a set of \(M_r\) elementary reaction steps, expressed as
where \([X_m]\) is the molar concentration of species \(m\), and \(\nu_{m,j}'\), \(\nu_{m,j}''\) are the stoichiometric coefficients on the reactant and product sides of reaction \(j\), associated with \(m\). For such a system, the rate of reaction \(j\) (\(R_j\)) can be expressed in terms of the the forward (\(k_{f,j}\)) and backward (\(k_{r,j}\)) rate coefficients,
The net molar production rate, :math:` dot{omega}_m` of species \(m\) is obtained by collating the rate of creation and destruction over reactions:
where \(\nu_{m,j} =\nu_{m,j}'' - \nu_{m,j}'\). Expressions for the reaction rates coefficients \(k_{(f,r),j}\) depend on the type of reaction considered. We use the CHEMKIN modified Arrhenius reaction format:
where \(A\) is the pre-exponential (frequency) factor, \(\beta\) is the temperature exponent and \(E_a\) is the activation energy. The CHEMKIN format additionally allows for a number of specializations of this format to represent pressure dependencies and third-body enhancements – see the CHEMKIN Manual or Cantera website for additional information.
Most fundamental Arrhenius reactions are bidirectional, and typically only the forward rates are specified. In this case, the balance of forward and reverse rates are dictacted by equilibrium thermodynamics, via the equilibrium constant, \(K_{c,j}\). In a low Mach system, \(K_{c,j}\) is a function only of temperature and the thermodynamic properties of the reactants and products of reaction \(j\),
\(\Delta H_j\) and \(\Delta S_j\) are the change in enthalpy and entropy of the reaction \(j\), and \(P_0\) is the ambient thermodynamic pressure.
Species production rates are evaluated via functions that are generated as part of a preprocessing step managed by the tool FUEGO.
Thermodynamic properties
Currently, expressions for the thermodynamic properties in PeleLM follow those of CHEMKIN, which assume a mixture of ideal gases. Species enthalpies and entropies are thus functions of only temperature (for perfect gases, they are independent of pressure) and are given in terms of polynomial fits to the species molar heat capacities (\(C_{p,\cdot}\)),
where, in the standard CHEMKIN framework (the 7-coefficients NASA format), \(N =5\) and
Accordingly, the standard-state molar enthalpy of species \(m\) is given by:
Note that the standard specifies that the heat of formation for the molecule is included in this expression. Similarly, the standard-state molar entropy is written as:
For each species \(m\), in the model the user must specify the 7 \(k\) coefficients \(a_{k,m}\). All other required thermodynamic properties are then determined (see, e.g., the CHEMKIN manual for additional details). Thermodynamic properties of the species, and those of the mixture, are evaluated via functions that are generated as part of a preprocessing step managed by the tool FUEGO.
FUEGO chemistry preprocessing
A typical model for PeleLM contains all the information associated with the CHEMKIN parameterization of the Arrhenius reaction set, as well as fitting coefficients for the thermodynamic relationships, and the specification of the species including data required to compute pure-species transport properties. In the combustion community, this information is communicated for each complete model –or mechanism, through multiple text files that conform to the CHEMKIN standards. The CHEMKIN driver code (or equivalent) can then be used to ingest the large number of parameters contained in these files and provide a set of functions for evaluating all the properties and rates required. Earlier versions of PeleLM linked to the CHEMKIN codes directly (and thereby assumed that all problems consisted of a mixture of ideal gases). However, evaluations were not very efficient because the functions stepped through generic expressions that included a large number of conditional statements and unused generality. Direct evaluation of these complex expressions allows for a much more efficient code that optimizes well with modern compilers. This is important because an appreciable fraction of PeleLM runtime is spent in these functions. Performance issues notwithstanding, customized evaluators will be necessary to extend PeleLM to a larger class of (real) gas models outside the CHEMKIN standard, such as SRK, that are already part of the PeleC code capabilities (PeleC shares use of PelePhysics for combustion model specification).
For these reasons, PeleLM no longer uses CHEMKIN functions directly, but instead relies on a preprocessing tool, FUEGO, to generate highly efficient C code implementations of the necessary thermodynamic, transport and kinetics evaluations. The source code generated from FUEGO is linked into the PeleLM executable, customizing each executable for a specific model at compile time. The implementation source code files can also be linked conveniently to post-processing analysis tools. The FUEGO processing tool, and the functions necessary to interface the generated functions to PeleLM are distributed in the auxiliary code package, PelePhysics. Included in the PelePhysics distribution is a broad set of models for the combustion of hydrogen, carbon-monoxide, methane, heptane, \(n\)-dodecane, dimethyl ether, and others, as well as instructions for users to extend this set using FUEGO, based on their own CHEMKIN-compliant inputs. PelePhysics also provides support for simpler gama-law equations-of-state, and simple/constant transport properties.
The PeleLM temporal integration
The temporal discretization in PeleLM combines a modified spectral deferred correction (SDC) coupling of chemistry and transport with a density-weighted approximate projection method for low Mach number flow. The projection method enforces a constrained evolution of the velocity field, and is implemented iteratively in such a way as to ensure that the update simultaneously satisfies the equation of state and discrete conservation of mass and total enthalpy. A time-explicit approach is used for advection; faster diffusion and chemistry processes are treated time-implicitly, and iteratively coupled together within the deferred corrections strategy. The integration algorithm, discussed in the following sections, is second-order accurate in space and time, and is implemented in the context of a subcycled approach for a nested hierarchy of mesh levels, where each level consists of logically rectangular patches of rectangular cells. All cells at a level have the same size in all coordinates.
Due to the complexity of the PeleLM algorithm, it is best presented in a number of passes. Focusing first on the single-level advance, we begin with a general discussion of the SDC-based time step iteration, which is designed to couple together the various physics processes. We then describe the projection steps used to enforce the constraint in the context of this iterative update. Next, we dive a little deeper into precisely how the advance of the thermodynamic components of the state is sequenced. There are a few crucial nuances to the formulation/sequencing of the energy advection, energy diffusion, conservative corrections to the species diffusion fluxes, and of the projection that can then be discussed in the context of overall single-level time step. Finally, with all these aspects defined, we give an overview of the modifications necessary to support the AMR subcycling strategy.
SDC preliminaries
The basic idea of SDC is to write the solution of an ODE
as an integral,
where we suppress explicit dependence of \(F\) and \(\phi\) on \(t\) for notational simplicity. Given an approximation \(\phi^{(k)}(t)\) to \(\phi(t)\), one can then define a residual,
Defining the error as \(\delta^{(k)}(t) = \phi(t) - \phi^{(k)}(t)\), one can then show that
In SDC algorithms, the integral in the above equation is evaluated with a higher-order quadrature rule. By using a low-order discretization of the integral one can construct an iterative scheme that improves the overall order of accuracy of the approximation by one per iteration, up to the order of accuracy of the underlying quadrature rule used to evaluate the integral. Specifically, if we let \(\phi^{(k)}\) represent the current approximation and define \(\phi^{(k+1)} = \phi^{(k)} + \delta^{(k)}\) to be the iterative update, then arrive at the update equation,
where a low-order discretization (e.g., forward or backward Euler) is used for the first integral and a higher-order quadrature is used to evaluate the second integral. For our reacting flow model, the underlying projection methodology for the time-advancement of velocity is second-order, so we require the use of second-order (or higher) numerical quadrature for the second integral.
MISDC Correction Equations
We based the time advance here on a variant of SDC, referred to as MISDC, in which \(F\) is decomposed into distinct processes, each treated separately with methods appropriate to its own time scale. Here, we write
to refer to advection, diffusion, and reaction processes. For this construction we assume that we are given an approximate solution \(\phi^{(k)}\) that we want to improve. A series of correction equations is developed to update \(\phi^{(k)}\) that uses relatively simple second-order discretizations of \(A(\phi)\) and \(D(\phi)\) but a high-accuracy treatment of \(R(\phi)\). In our approach, \(A(\phi^{(k)})\) is piecewise-constant over each time step, and is evaluated using a second-order Godunov procedure. The Godunov procedure computes a time-centered advection term at \(t^{n+1/2}\), and incorporates an explicit diffusion source term and an iteratively lagged reaction source term, i.e.,
where \(I_R^{(k-1)}\) is the effective contribution due to reactions from the previous iteration, i.e.,
where \(\Delta t^n = t^{n+1} - t^n\). Here \(I_R^{(k-1)}\) is computed from a high-accuracy integration of the reaction kinetics equations, augmented with piecewise constant-in-time representation of advection and diffusion. Details of this procedure are given below.
We also represent \(D(\phi^{(k)})\) as piecewise constant over the time step, using a mid-point rule:
In the spirit of MISDC, we solve correction equations for the individual processes sequentially. We begin by discretizing the update equation, but only including the advection and diffusion terms in the correction integral,
Thus, \(\phi_{\rm AD}^{(k+1)}(t)\) represents an updated approximation of the solution after correcting the advection and diffusion terms only. For the first integral, we use an explicit update for the advection term and a backward Euler discretization for the diffusion term. For the second integral, we represent \(F\) in terms of \(A\), \(D\), and \(R\) and use the definition of \(A^{(k)}\), \(D^{(k)}\), and \(I_R^{(k-1)}\) to obtain a discrete update for \(\phi_{\rm AD}^{n+1,(k+1)}\):
This equation simplifies to the following backward Euler type linear system, with the right-hand-side consisting of known quantities:
After computing \(\phi_{\rm AD}^{n+1,(k+1)}\), we complete the update by solving a correction equation for the reaction term. Standard MISDC approaches would formulate the reaction correction equation as
and use a backward Euler type discretization for the integral of the reaction terms. Here, to address stiffness issues with detailed chemical kinetics, we will instead formulate the correction equation for the reaction as an ODE, which is treated separately with an ODE integrator package. In particular, by differentiating the SDC update we obtain
which we then advance with the ODE integrator over \(\Delta t\) to obtain \(\phi^{n+1,(k+1)}\). After the integration, we can evaluate \(I_R^{(k+1)}\), which is required for the next iteration
Summarizing, the variant of SDC used in the single-level time-step of PeleLM integrates the \(A\), \(D\) and \(R\) components of the discretization scheme in an iterative fashion, and each process incorporates a source term that is constructed using a lagged approximation of the other processes. In the case of the implicit diffusion, an additional source term arises from the SDC formulation. If the SDC iterations were allowed to fully converge, all the process advanced implicitly would be implicitly coupled to all others. Moreover, each process is discretized using methods that are tailored specifically to the needs of that operator. In the next section, we give more details for each of the components, including how and where the velocity projections play a role.
Data centering, \(A\)-\(D\)-\(R\), and the projections
PeleLM implements a finite-volume, Cartesian grid discretization approach with constant grid spacing, where \(U\), \(\rho\), \(\rho Y_m\), \(\rho h\), and \(T\) represent cell averages, and the pressure field, \(\pi\), is defined on the nodes of the grid, and is temporally constant on the intervals over the time step. There are three major steps in the algorithm:
Step 1: (Compute advection velocities) Use a second-order Godunov procedure to predict a time-centered velocity, \(U^{{\rm ADV},*}\), on cell faces using the cell-centered data (plus sources due to any auxiliary forcing) at \(t^n\), and the lagged pressure gradient from the previous time interval, which we denote as \(\nabla \pi^{n-1/2}\). The provisional field, \(U^{{\rm ADV},*}\), fails to satisfy the divergence constraint. We apply a discrete projection by solving the elliptic equation with a time-centered source term:
for \(\phi\) at cell-centers, where \(D^{{\rm FC}\rightarrow{\rm CC}}\) represents a cell-centered divergence of face-centered data, and \(G^{{\rm CC}\rightarrow{\rm FC}}\) represents a face-centered gradient of cell-centered data, and \(\rho^n\) is computed on cell faces using arithmetic averaging from neighboring cell centers. Also, \(\widehat S\) refers to the RHS of the constraint equation, with adjustments to be discussed in the next section – these adjustments are computed to ensure that the final update satisfied the equation of state. The solution, \(\phi\), is then used to define
After the MAC-projection, \(U^{\rm ADV}\) is a second-order accurate, staggered grid vector field at \(t^{n+1/2}\) that discretely satisfies the constraint. This field is the advection velocity used for computing the time-explicit advective fluxes for \(U\), \(\rho h\), and \(\rho Y_m\).
Step 2: (Advance thermodynamic variables) Integrate \((\rho Y_m,\rho h)\) over the full time step. The details of this are presented in the next subsection.
Step 3: (Advance the velocity) Compute an intermediate cell-centered velocity field, \(U^{n+1,*}\) using the lagged pressure gradient, by solving
where \(\tau^{n+1,*} = \mu^{n+1}[\nabla U^{n+1,*} +(\nabla U^{n+1,*})^T - 2\mathcal{I}\widehat S^{n+1}/3]\) and \(\rho^{n+1/2} = (\rho^n + \rho^{n+1})/2\), and \(F\) is the velocity forcing. This is a semi-implicit discretization for \(U\), requiring a linear solve that couples together all velocity components. The time-centered velocity in the advective derivative, \(U^{n+1/2}\), is computed in the same way as \(U^{{\rm ADV},*}\), but also includes the viscous stress tensor evaluated at \(t^n\) as a source term in the Godunov integrator. At this point, the intermediate velocity field \(U^{n+1,*}\) does not satisfy the constraint. Hence, we apply an approximate projection to update the pressure and to project \(U^{n+1,*}\) onto the constraint surface. In particular, we compute \(\widehat S^{n+1}\) from the new-time thermodynamic variables and an estimate of \(\dot\omega_m^{n+1}\), which is evaluated directly from the new-time thermodynamic variables. We project the new-time velocity by solving the elliptic equation,
for nodal values of \(\phi\). Here, \(L^{{\rm N}\rightarrow{\rm N}}\) represents a nodal Laplacian of nodal data, computed using the standard bilinear finite-element approximation to \(\nabla\cdot(1/\rho^{n+1/2})\nabla\). Also, \(D^{{\rm CC}\rightarrow{\rm N}}\) is a discrete second-order operator that approximates the divergence at nodes from cell-centered data and \(G^{{\rm N}\rightarrow{\rm CC}}\) approximates a cell-centered gradient from nodal data. Nodal values for \(\widehat S^{n+1}\) required for this equation are obtained by interpolating the cell-centered values. Finally, we determine the new-time cell-centered velocity field using
and the new time-centered pressure using \(\pi^{n+1/2} = \phi\).
Thus, there are three different types of linear solves required to advance the velocity field. The first is the MAC solve in order to obtain face-centered velocities used to compute advective fluxes. The second is the multi-component cell-centered solver is used to obtain the provisional new-time velocities. Finally, a nodal solver is used to project the provisional new-time velocities so that they satisfy the constraint.
Thermodynamic Advance
Here we describe the details of Step 2 above, in which we iteratively advance \((\rho Y_m,\rho h)\) over the full time step. We begin by computing the diffusion operators at \(t^n\) that will be needed throughout the iteration. Specifically, we evaluate the transport coefficients \((\lambda,C_p,\mathcal D_m,h_m)^n\) from \((Y_m,T)^n\), and the provisional diffusion fluxes, \(\widetilde{\boldsymbol{\cal F}}_m^n\). These fluxes are conservatively corrected (i.e., adjusted to sum to zero by adding a mass-weighted “correction velocity”) to obtain \({\boldsymbol{\cal F}}_m^n\) such that \(\sum {\boldsymbol{\cal F}}_m^n = 0\). Finally, we copy the transport coefficients, diffusion fluxes and the thermodynamic state from \(t^n\) as starting values for \(t^{n+1}\), and initialize the reaction terms, \(I_R\) from the values used in the previous step. The following sequence is then repeated for each iteration, \(k<k_{max}\)
Step 2-I: Use a second-order Godunov integrator to predict species time-centered edge states, \((\rho Y_m)^{n+1/2,(k+1)}\) and their advection terms at \(t^{n+1/2}\), \((A_m^{n+1/2,(k+1)})\). Source terms for this prediction include explicit diffusion forcing, \(D^{n}\), and an iteration-lagged reaction term, \(I_R^{(k)}\). Since the remaining steps of the algorithm for this iteration (including diffusion and chemistry advances) will not affect the new-time density for this iteration, we can already compute \(\rho^{n+1,(k+1)}\). This will be needed in the trapezoidal-in-time diffusion solves.
In addition to predicting \(\rho\) and \(\rho Y_m\) to the faces to compute advective fluxes, we also need \(\rho h\) there. We could use a Godunov scheme as well, however, because \(h\) contains the heat of formation scaled to an arbitrary reference state, it is not generally monotonic through flames. Also, because the equation of state is generally nonlinear, this will often lead to numerically-generated n on-mononoticity in the temperature field. An analytically equivalent approach, based on the fact that temperature should be smoother and monotonic through the flame, is to instead predict temperature with the Godunov scheme to the cell faces directly. Then, using \(T\), \(\rho = \sum (\rho Y_m)\) and \(Y_m = (\rho Y_m)/\rho\) on the cell faces directly, we can evaluate \(h\) instead of extrapolating. We can then evaluate the enthalpy advective flux divergence, \(A_h^{n+1/2,(k+1)}\), for \(\rho h\).
Step 2-II: Update the transport coefficients (if necessary) with the most current cell-centered thermodynamic state, then interpolate those values to the cell faces. We now compute provisional, time-advanced species mass fractions, \(\widetilde Y_{m,{\rm AD}}^{n+1,(k+1)}\), by solving a backward Euler type correction equation for the Crank-Nicolson update, using Eq. (1). Note that the provisional species diffusion fluxes reads \(\widetilde{\boldsymbol{\cal F}}_{m,{\rm AD}}^{(0)} = -\rho^n D_{m,mix}^n \nabla \widetilde X_{m,{\rm AD}}^{(0)}\). This expression couples together all of the species mass fractions (\(Y_m\)) in the update of each, even for the mixture-averaged model. Computationally, it is much more tractable to write this as a diagonal matrix update with a lagged correction by noting that \(X_m = (W/W_m)Y_m\). Using the chain rule, \(\widetilde{\boldsymbol{\cal F}}_{m,{\rm AD}}^{(0)}\) then has components proportional to \(\nabla Y_m\) and \(\nabla W\). The latter is lagged in the iterations, and is typically very small. In the limit of sufficient iterations, diffusion is driven by the true form of the the driving force, \(d_m\), but in this form, each iteration involves decoupled diagonal solves (following the SDC formalism used above):
where
By iteratively lagging the \(\nabla W\) term (and \(D_{m,mix}\)), this equation is a scalar, time-implicit, parabolic and linear for the updated \(\widetilde Y_{m,{\rm AD}}^{n+1,(k+1)}\) (and requires a linear solve). The form of this solve, from a software perspective, is identical to that of the MAC projection discussed above.
Once all the species equations are updated, we compute \({\boldsymbol{\cal F}}_{m,{\rm AD}}^{n+1,(k+1)}\), which are conservatively corrected versions of \(\widetilde{\boldsymbol{\cal F}}_{m,{\rm AD}}^{n+1,(k+1)}\), and then the species mass fractions are updated too, using
where
Next, we compute the time-advanced enthalpy, \(h_{\rm AD}^{n+1,(k+1)}\). Much like for the diffusion of the species densities, \(\rho Y_m\), where a \(\nabla X_m\) driving force leads to a nonlinear, coupled Crank-Nicolson update; the enthalpy diffuses with a \(\nabla T\) driving force. We define an alternative linearized strategy. We begin by following the same SDC-correction formalism used for the species, and write the nonlinear update for \(\rho h\) (noting that there is no reaction source term here):
where
However, since we cannot compute \(h_{{\rm AD}}^{n+1,(k+1)}\) directly, we solve this iteratively based on the approximation \(h_{{\rm AD}}^{(k+1),\ell+1} \approx h_{{\rm AD}}^{(k+1),\ell} + C_{p}^{(k+1),\ell} \delta T^{(k+1),\ell+1}\), with \(\delta T^{(k+1),\ell+1} = T_{{\rm AD}}^{(k+1),\ell+1} - T_{{\rm AD}}^{(k+1),\ell}\), and iteration index, \(\ell\) = 1:\(\,\ell_{MAX}\). The enthalpy update equation is thus recast into a linear equation for \(\delta T^{(k+1);\ell+1}\)
where \(H_{AD}^{(k+1),\ell} = - \nabla \cdot \sum h_m(T_{AD}^{(k+1),\ell}) \, {\boldsymbol{\cal F}}_{m,AD}^{n+1,(k+1)}\) and \(D_{T,AD}^{(k+1),\ell} = \nabla \cdot \lambda^{(k)} \, \nabla T_{AD}^{(k+1),\ell}\). Note that again the solve for this Crank-Nicolson update has a form that is identical to that of the MAC projection discussed above. After each iteration, update \(T_{{\rm AD}}^{(k+1),\ell+1} = T_{{\rm AD}}^{(k+1),\ell} + \delta T^{(k+1),\ell+1}\) and re-evaluate \((C_p ,h_m)^{(k+1),\ell+1}\) using \((T_{{\rm AD}}^{(k+1),\ell+1}, Y_{m,{\rm AD}}^{n+1,(k+1)}\)).
Step 2-III: Based on the updates above, we define an effective contribution of advection and diffusion to the update of \(\rho Y_m\) and \(\rho h\):
that we treat as piecewise-constant source terms to advance \((\rho Y_m,\rho h)^n\) to \((\rho Y_m,\rho h)^{n+1,(k+1)}\). The ODE system for the reaction part over \(\Delta t^n\) then takes the following form:
After the integration is complete, we make one final call to the equation of state to compute \(T^{n+1,(k+1)}\) from \((Y_m,h)^{n+1,(k+1)}\). We also can compute the effect of reactions in the evolution of \(\rho Y_m\) using,
If \(k<k_{\rm max}-1\), set \(k=k+1\) and return to Step 2-I. Otherwise, the time-advancement of the thermodynamic variables is complete, and set \((\rho Y_m,\rho h)^{n+1} = (\rho Y_m,\rho h)^{n+1,(k+1)}\). If \(k+1 = k_{max}\), Step 2 of our algorithm is complete.
Modifications for AMR
The framework to manage adaptive mesh refinement (AMR) used in PeleLM borrows heavily from the AMReX library, and the IAMR code; the reader is referred to documentation of both of these components in order to understand the distributed, logically rectangular data structures used, and the recursive time-stepping strategy for advancing a hierarchy of nested grid levels.
Summarizing, there is a bulk-synchronous advance of each level over its respective time step, \(dt\), followed recursively by a number of (sub-)steps of the next-finer AMR level. Each fine level advanced is over an interval \((1/R) dt\), if the fine cells are a factor of \(R\) smaller, and in this scenario, the coarser level provides Dirichlet boundary condition data for the fine-level advances. Note that the levels are properly nested so that the finer level is fully contained within the coarser level, except perhaps at physical boundaries, where their edges can be coincident - thus, the fine level has sufficient boundary data for a well-posed advance.
After two adjacent levels in the hierarchy reach the same physical time, a synchronization operation is performed to ensure that the coarse data is consistent with the volume integral of the fine data that covers it, and the fluxes across of the coarse-fine interface are those of the fine solution. The latter of these two operations can be quite complex, as it must correct coarse-grid errors committed by each of the operators used to perform the original advance. It may also be non-local, in that cells far away from the coarse-fine interface may need to incorporate flux increments due to the mismatched coarse and fine solutions. Formally, the synchronzation is a bilevel correction that should be computed as a sequence of two-level solves. However, this would lead to the same amound of work that was required to create the original (pre-sync) data. We assume that the corrections computed for the synchronization are smooth enough to be well represented by an increment on the coaser of the two-levels, and interpolated to the finer grid. Note that the transport coefficients are not updated to account for the state changes during the synchronization.
Generically, the synchronization procedure in PeleLM follows that described for the IAMR code, but with modifications to explicitly enforce that the sum of the species diffusion correction fluxes is zero, that the nonlinear enthalpy update is solved (similar to described above for the single-level advance), and the corection for the advection velocity is adjusted iteratively so that the final synchronized state satisfies the EOS.
There are several components in PeleLM that contribute to the flux mismatch at the coarse-fine interface. The first component arise from the face-centered velocity, \(U^{ADV,\ell}\), used to advect the scalars at each AMR level \(\ell\), since the field satifies a divergence constraint on the coarse and fine levels separately. We compute a velocity mismatch
(where \(d\) is the number of spatial dimensions) along the coarse-fine boundary. We then solve the elliptic projection equation
where \(\delta_{\chi}^{\ell}\) is incremented iteratively to enfoce the final state to satisfy the EOS and compute the correction velocity
which is the increment of velocity required to carry advection fluxes needed to correct the errors made by advancing the coarse state with the wrong velocities.
The second part of the mismatch arises because the advective and diffusive fluxes on the coarse grid were computed without explicitly accounting for the fine grid, while on the fine grid the fluxes were computed using coarse-grid Dirichlet boundary data. We define the flux discrepancies on the coarser level \(\ell\) of the pair of levels considered:
where \(\boldsymbol{\cal F}\) is the total (advective + diffusive) flux through a face on the coarse-fine interface prior the synchronization operations. Since all operations are performed on the coarse level we will drop the \(\ell\) in the following.
Since mass is conserved, corrections to density, \(\delta \rho^{sync}\) on the coarse grid associated with mismatched advection fluxes may be computed explicitly
We can compute the post-sync new-time value of density, \(\rho^{n+1} = \rho^{n+1,p} + \delta \rho^{sync}\), where \(^p\) denotes pre-sync quantities. The synchronization correction of a state variable \(\delta ( \rho \phi )^{sync}\) (where \(\phi \in (Y_m, h)\)) is obtained by subtracting the pre-sync state value \(( \rho \phi )^{n+1,p}\) from the corrected one \(( \rho \phi )^{n+1}\), both of which expressed from an SDC iteration update (see Eq. (2) and (3)) but with the divergence of the correction velocity fluxes (\(U^{ADV,corr}\)) and fluxes mismatch (\(\delta \boldsymbol{\cal F}\)) included in the advection and diffusion corrected operators.
Given the corrected density we can decompose the sync corrections \(\delta ( \rho \phi )^{sync} = \phi^{n+1,p} \delta \rho^{sync} + \rho^{n+1} \delta \phi^{sync}\) and obtain the linear system for \(\delta \phi^{sync}\) since the fluxes mismatch contain implicit diffusion fluxes from the Crank-Nicolson discretization. For species \(m\) the implicit system reads:
where \(\widetilde{\cal F_{m}}\) is the species correction flux due to the sync correction, \(\widetilde{\delta Y_m^{sync}}\). However, as in the single-level algorithm, the species fluxes must be corrected to sum to zero. These adjusted fluxes are then used to recompute a \(\delta Y_m^{sync}\), which is then used via the expression above to compute \(\delta (\rho Y_m)^{sync}\), the increment to the species mass densities.
In order to get the equation for the enthalpy sync correction, we operate as for species mass fractions. We will present the details of the method. The SDC advection-diffusion udpate for pre-sync enthalpy is (3) (now including the superscript \(p\)) and its corrected counterpart reads:
where
and with \(\delta T^{sync} = T^{n+1,(k+1)} - T^{n+1,(k+1),p}\). Subtracting the pre-sync eq. (3) (with the superscript \(p\)) from the above equation (5) and gathering the flux mismatches and correction velocity fluxes in \(S_h^{sync}\) we obtain:
where
and the updated (\(n+1\)) fluxes (without the \(*\)) only now contain the implicit contribution, i.e:
To go further we note that:
where \(\delta h_m^{sync} = h_m(T^{n+1,(k+1)}) - h_m(T^{n+1,(k+1),p})\) and \(\delta {\boldsymbol{\cal F}}_m^{sync}\) is the species flux increment due to the species sync correction appearing on the LHS of eq. (4). Eq. (6) is the equation for the sync correction. At this point, we can drop the SDC iteration index \(k+1\) for simplicity (all \(k\) related quantities are contained in \(S_h^{sync}\)). Note that the evaluation of the transport properties is relatively expensive, such that we don’t want to update the conductivity in \(D_T^{n+1}\) since a lagged (pre-sync) version is sufficient for second-order accuracy. However we do want to use an updated version of \(h_m\).
Just as in the level advance, we cannot compute \(h^{n+1}\) directly, so we solve this iteratively based on the approximation \(h^{n+1,\eta+1} \approx h^{n+1,\eta} + C_{p}^{n+1,\eta} \Delta T^{\eta+1}\), with \(\Delta T^{\eta+1} = T^{n+1,\eta+1} - T^{n+1,\eta}\), and iteration index, \(\eta\) = 1:\(\,\eta_{MAX}\). The sync equation is thus recast into a linear equation for \(\Delta T^{\eta+1}\), and we lag the \(H\) terms in iteration \(\eta\),
After each :math:eta iteration, update \(T^{n+1,\eta+1} = T^{n+1,\eta} + \Delta T^{\eta+1}\), \(\delta T^{sync,\eta+1} = T^{n+1,\eta+1} - T^{(n+1,p)}\), and re-evaluate \((C_p,h_m)^{n+1,\eta+1}\) using \((T^{n+1,\eta+1}, Y_m^{n+1}\)). Iterations are continued until the norm of \(\Delta T^{\eta+1}\) drops below a tolerance threshold. Then set \(T^{n+1} = T^{(n+1,p)} + \delta T^{sync,\eta_{MAX}}\), and compute \(h^{n+1} = h(T^{n+1},Y_m^{n+1})\).