The PeleLMeX Model
In this section, we present the actual model that is evolved numerically by PeleLMeX, 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 PeleLMeX, 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 PeleLMeX in later sections.
Overview of PeleLMeX
PeleLMeX is the non-subcycling version of PeleLM a parallel, adaptive mesh refinement (AMR) code that solves the reacting Navier-Stokes equations in the low Mach number regime.
In a nutshell, PeleLMeX features include:
- Physics & numerics :
Finite volume, block-structured AMR approach
2D-Cartesian, 2D-Axisymmetric and 3D support
Combustion (transport, kinetics, thermodynamics) models based on Cantera and EGLib through PelePhysics
Complex geometries using Embedded Boundaries (EB)
Second-order projection methodology for enforcing the low Mach number constraint
Time advance based on a spectral-deferred corrections approach that conserves species mass and energy and evolves on the equation of state
Several higher-order Godunov integration schemes for advection
Temporally implicit viscosity, species mass diffusion, thermal conductivity, chemical kinetics
Closed chamber algorithm enable time-varying background pressure changes
Lagrangian spray description and Hybrid Method of Moments Soot modeling using PelePhysics (formerly these were part of PeleMP).
Mathematical background
PeleLMeX 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. PeleLMeX also utilizes the source code and algorithmic infrastructure of AMReX-Hydro. PeleLMeX borrows heavily from PeleLM. The core algorithms in PeleLM 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
PeleLMeX 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)\). Note that \(p_0\) is spatially constant, but can be time-dependent as described in the next section.
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 source terms, \(S_{\text{ext}}\), are user defined external forcing terms. For example, we have used \(\boldsymbol{S}_{\text{ext},\rho\boldsymbol{u}} = \rho\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,
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 takes 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.
For the standard ideal gas EOS,
Therefore, the divergence constraint on velocity becomes:
However, it can be shown that
and
Thus, the terms containing \(S_{\text{ext},\rho}\) cancel and the divergence constraint for the standard ideal gas EOS simplifies to:
In addition to the flow equations, PeleLMeX can also solve for a set of quantities that are neither advected nor diffused, satisfying:
Confined domain ambient pressure
In unconfined domains, the ambient pressure will remain constant in time, but for confined domains, this is not the case. Above, we assumed that \(p_0\) was constant. If \(p_0\) is a function of time, the pressure derivative term must be restored in the velocity divergence constraint as:
where \(\theta \equiv 1/(\Gamma_1 p_0)\), with \(\Gamma_1 = \partial ln(p)/\partial ln(\rho)|_s\) being the first adiabatic exponent. \(\Gamma_1\) depends on the composition and is not a constant. Both \(\theta\) and \(\S\) can be decomposed into mean and fluctuating components and the above equation can be rewritten as:
where \(\overline \theta\) and \(\overline S\) are the mean values of \(\theta\) and \(S\) over the domain, and \(\delta \theta\) and \(\delta S\) are the perturbations off their respective means that both integrate to zero over the domain, by definition. This equation can be simplified by integrating over the domain volume:
Since the perturbations integrate to zero over the domain volume, the mean values are constants, and \(p_0\) is only a function of time, the above simplifies to:
Solving for \(dp_0/dt\) yields an evolution equation of \(p_0\):
where we have used the divergence theorem to convert the volume integral into a surface integral over the domain boundaries: \(\int_V \nabla \cdot \boldsymbol{u} dV = \int_A \boldsymbol{u} dA\). The above pressure evolution is accompanied by a modified velocity constraint:
The above equations hold for any fully enclosed or partially enclosed domain where there can be mass flowing into or out of the domain, but the next flowrate is non-zero and therefore the pressure should be adjusted in time. In a perfectly enclosed domain, where there is no mass in or out of the system, \(\int_A \boldsymbol{u} dA = 0\) and the pressure evolution is simplified to:
and simplified velocity constraint,
PeleLMeX Algorithm
An overview of PeleLMeX time-advance function is provided in Fig. 1 and details are provided in the following subsections.
Fig. 1 : Flowchart of the PeleLMeX advance function.
The three steps of the low Mach number projection scheme described below are referenced to better emphasize how the thermodynamic solve is closely weaved into the fractional step approach. Striped boxes indicate where the Godunov procedure is employed while the four different linear solves are highlighted.
Low Mach number projection scheme
PeleLMeX 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. The projection scheme is based on a fractional step approach where, for purely incompressible flow, the velocity is first advanced in time using the momentum equation (Step 1) and subsequently projected to enforce the divergence constraint (Step 3). When considering variable density flows, the scalar thermodynamic advance is performed between these two steps (Step 2), but within the SDC context, Step 1 and Step 2 are effectively interlaced. The three major steps of the algorithm (Almgren et al. 1998, Day and Bell, 2000, Nonaka et al. 2012):
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}\). This provisional field, \(U^{{\rm ADV},*}\), fails to satisfy the divergence constraint. We apply a discrete projection (termed MAC-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, \(S^{MAC}\) 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 (face-centered) 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 using a spectral deferred correction (SDC) approach, the details of which can be found in PeleLM documentation. An even more detailed version of the algorithm is available in Nonaka et al., 2018.
We begin by computing the diffusion terms \(D^n\) at \(t^n\) that will be needed throughout the SDC iterations. 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,(k=0)}\), 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}\) starting at \(k=0\):
if \(k>0\), compute the lagged (previous \(k\) iteration) transport properties, diffusion terms \(D^{n+1,(k)}\) and divergence constraint \(\widehat S^{n+1,(k)}\)
construct the MAC-projection RHS \(S^{MAC}\), combining \(t^n\) and \(t^{n+1,(k)}\) estimates of \(\widehat S\), and the pressure correction \(\chi\) (Nonaka et al, 2018):
\[S^{MAC} = \frac{1}{2}(\widehat S^n + \widehat S^{n+1,(k)}) + \sum_{i=0}^k \frac{1}{p_{therm}^{n+1,(i)}}\frac{p_{therm}^{n+1,(i)}-p_0}{\Delta t}\]In this update, it is optional whether to update the \(\widehat S^{n+1}\) term on every SDC iteration, or to simply compute it for \(k = 0\) and then hold it constant, with the \(\chi\) correction iterations accounting for changes during the SDC iterations. The latter strategy has been observed to improve convergence in some cases.
Perform Step 1 to obtain the time-centered, staggered \(U^{ADV}\)
Use a 2nd Godunov integrator to predict species time-centered edge states, \((\rho Y_m)^{n+1/2,(k+1)}\) and their advection terms \(A_m^{n+1/2,(k+1)}\) using \(U^{ADV}\). Source terms for this prediction include explicit diffusion forcing, \(D^{n}\), and an iteration-lagged reaction term, \(I_R^{(k)}\). Since the diffusion and chemistry will not affect the new-time density, we can already compute \(\rho^{n+1,(k+1)}\). This will be needed in the trapezoidal-in-time diffusion solves. We also compute \(A_h^{n+1/2,(k+1)}\): we could also use a Godunov scheme, but 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 non-mononoticity in the temperature field. Using the fact that temperature should be smoother and monotonic through the flame, we instead predict temperature with the Godunov scheme and use face-centered \(T\), \(\rho = \sum (\rho Y_m)\) and \(Y_m = (\rho Y_m)/\rho\) to evaluate \(h\) instead of extrapolating. We can then evaluate the enthalpy advective flux divergence, \(A_h^{n+1/2,(k+1)}\), for \(\rho h\).
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. Note that the provisional species diffusion fluxes reads \(\widetilde{\boldsymbol{\cal F}}_{m,{\rm AD}}^{(k)} = -\rho^n D_{m,mix}^n \nabla \widetilde X_{m,{\rm AD}}^{(k)}\). 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}}^{(k)}\) 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):
\[\frac{\rho^{n+1,(k+1)}\widetilde Y_{m,{\rm AD}}^{n+1,(k+1)} - (\rho Y_m)^n}{\Delta t} = A_m^{{n+1/2,(k+1)}} + \widetilde D_{m,AD}^{n+1,(k+1)} + \frac{1}{2}(D_m^n - D_m^{n+1,(k)}) + I_{R,m}^{(k)}\]The resulting \(\rho^{n+1,(k+1)}\widetilde Y_{m,{\rm AD}}^{n+1,(k+1)}\) are used to compute the implicit (conservatively-corrected) species diffusion fluxes and implicit diffusion term \(D_{m,AD}^{n+1,(k+1)}\), which is employed to get a final AD updated \(\rho^{n+1,(k+1)}\widetilde Y_{m,{\rm AD}}^{n+1,(k+1)}\). Next, we compute the time-advanced enthalpy, \(h_{\rm AD}^{n+1,(k+1)}\). Much like for the diffusion of the \(\rho Y_m\), the \(\nabla T\) driving force leads to a nonlinear, coupled Crank-Nicolson update for \(\rho h\). We define an alternative linearized strategy 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):
\[\begin{split} \frac{\rho^{n+1,(k+1)} h_{{\rm AD}}^{n+1,(k+1)} - (\rho h)^n}{\Delta t} = A_h^{n+1/2,(k+1)} + D_{T,AD}^{n+1,(k+1)} + H_{AD}^{n+1,(k+1)} \\ + \frac{1}{2} \Big( D_T^n - D_T^{n+1,(k)} + H^n - H^{n+1,(k)} \Big)\end{split}\]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}\):
\[\begin{split}\rho^{n+1,(k+1)} C_p^{(k+1),\ell} \delta T^{(k+1),\ell+1} -\Delta t \, \nabla \cdot \lambda^{(k)} \nabla (\delta T^{(k+1),\ell +1}) = \rho^n h^n - \rho^{n+1,(k+1)} \\ h_{AD}^{(k+1),\ell} + \Delta t \Big( A_h^{n+1/2,(k+1)} + D_{T,AD}^{(k+1),\ell} + H_{AD}^{(k+1),\ell} \Big) + \frac{\Delta t}{2} \Big( D_T^n - D_T^{n+1,(k)} + H^n - H^{n+1,(k)} \Big)\end{split}\]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}\). 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)}\)).
Based on the updates above, we define an effective contribution of advection and diffusion to the update of \(\rho Y_m\) and \(\rho h\):
\[\begin{split}&&Q_{m}^{n+1,(k+1)} = A_m^{n+1/2,(k+1)} + D_{m,AD}^{(n+1,k+1)} + \frac{1}{2}(D_m^n - D_m^{n+1,(k)}) \\ &&Q_{h}^{n+1,(k+1)} = A_h^{n+1/2,(k+1)} + D_{T,AD}^{n+1,(k+1)} + H_{AD}^{n+1,(k+1)} + \frac{1}{2}(D_T^n - D_T^{n+1,(k)} + H^n - H^{n+1,(k)} )\end{split}\]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:
\[\begin{split}\frac{\partial(\rho Y_m)}{\partial t} &=& Q_{m}^{n+1,(k+1)} + \rho\dot\omega_m(Y_m,T),\label{eq:MISDC VODE 3}\\ \frac{\partial(\rho h)}{\partial t} &=& Q_{h}^{n+1,(k+1)}.\label{eq:MISDC VODE 4}\end{split}\]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,
\[I_{R,m}^{(k+1)} = \frac{(\rho Y_m)^{n+1, (k+1)} - (\rho Y_m)^n}{\Delta t} - Q_{m}^{n+1,(k+1)}.\]If \(k=k_{max}-1\), the time-advancement of the thermodynamic variables is complete, set \((\rho Y_m,\rho h)^{n+1} = (\rho Y_m,\rho h)^{n+1,(k+1)}\).
Before moving to Step 3, the new time viscosity and instantaneous divergence constraint \(\widehat S^{n+1}\) are evaluated.
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 - \frac{2}{3} \mathcal{I} \, \nabla \cdot U^{n+1,*}]\) 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.
Advection schemes
PeleLMeX relies on the AMReX-Hydro implementation of the 2nd-order Godunov method, with several variants available. The basis of the Godunov approach is to extrapolate the cell-centered quantity of interest (\(U\), \(\rho Y\), \(\rho h\)) to cell faces using a second-order Taylor series expansion in space and time. As detailed in AMReX-Hydro documentation, the choice of the slope order and limiting scheme define the exact variant of the Godunov method. Of particular interest for combustion applications, where sharp gradients of intermediate chemical species are found within flame fronts, the Godunov_BDS approach provides a bound-preserving advection scheme which greatly limits the appearance of over-/under-shoots, often leading to critical failure of the stiff chemical kinetic integration.
Note that in the presence of EB, only the Godunov_PLM variant is available.
AMR extension
In contrast with PeleLM, PeleLMeX do not rely a on subcycling approach to advance the AMR hierarchy. This difference is illustrated in the figure below comparing the multi-level time-stepping approach in both codes:
PeleLM will recursively advance finer levels, halving the time step size (when using a refinement ratio of 2) at each level. For instance, considering a 3 levels simulation, PeleLM advances the coarse Level0 over a \(\Delta t_0\) step, then Level1 over a \(\Delta t_1\) step and Level2 over two \(\Delta t_2\) steps, performing an interpolation of the Level1 data after the first Level2 step. At this point, a synchronization step is performed to ensure that the fluxes are conserved at coarse-fine interface and a second Level1 step is performed, followed by the same two Level2 steps. At this point, two synchronizations are needed between the two pairs of levels.
In order to get to the same physical time, PeleLMeX will perform 4 time steps of size similar to PeleLM’s \(\Delta t_2\), advancing all the levels at once. The coarse-fine fluxes consistency is this time ensured by averaging down the face-centered fluxes from fine to coarse levels. Additionally, the state itself is averaged down at the end of each SDC iteration.
In practice, PeleLM will perform a total of 7 single-level advance steps, while PeleLMeX will perform 4 multi-level ones to reach the same physical time, advancing the coarser levels at a smaller CFL number whereas PeleLM maintain a fixed CFL at all the level. It might seem that PeleLMeX is thus performing extra work, but because it ignore fine-covered regions, PeleLMeX do not need to perform the expensive (and often very under-resolved) chemistry integration in fine-covered areas. An exact evaluation of the benefits and drawbacks of each approach is under way.
Geometry with Embedded Boundaries
PeleLMeX relies on AMReX’s implementation of the Embedded Boundaries (EB) approach to represent geometrical objects. In this approach, the underlying computational mesh is uniform and block-structured, but the boundary of the irregular-shaped computational domain conceptually cuts through this mesh. Each cell in the mesh becomes labeled as regular, cut or covered, and the finite-volume based discretization methods traditionally used in AMReX applications need to be modified to incorporate these cell shapes. AMReX provides the necessary EB data structures, including volume and area fractions, surface normals and centroids, as well as local connectivity information. The fluxes described in the projection scheme section are then modified to account for the aperture opening between adjacent cells and the additional EB-fluxes are included when constructing the cell flux divergences.
A common problem arising with EB is the presence of the small cut-cells which can either introduce undesirable constraint on the explicit time step size or lead to numerical instabilities if not accounterd for. PeleLMeX relies on a combination of classical flux redistribution (FRD) (Pember et al, 1998) and state redistribution (SRD) (Giuliani et al., 2022) to circumvent the issue. In particular, explicit advective fluxes \(A^{n+1/2,(k+1)}\) are treated using SRD while explicit diffusion fluxes \(D^{n}\) and SDC iteration-lagged \(D^{n+1,(k)}\) are treated with FRD. Note that implicit diffusion fluxes are not redistributed as AMReX’s linear operators are EB-aware.
The use of AMReX’s multigrid linear solver introduces constraint on the complexity of the geometry PeleLMeX is able to handle. The efficiency of the multigrid approach relies on generating coarse version of the linear problem. If the geometry includes thin elements (such as tube or plate) or narrow channels, coarsening of the geometry is rapidly limited by the occurrence of multi-cut cells (not supported by AMReX) and the linear solvers are no longer able to robustly tackle projections and implicit diffusion solves. AMReX include an interface to HYPRE which can help circumvent the issue by sending the coarse-level geometry directly to HYPRE algebraic multigrid solvers. More details on how to use HYPRE is provided in control Section.
Large Eddy Simulation
To provide closure for the unresolved turbulent stress/flux terms in Large Eddy Simulation (LES), PeleLMeX supports the constant-coefficient Smagorinsky, WALE and Sigma models for turbulent transport of momentum, species, and energy. These models are based on a gradient transport assumption, resulting in terms analogous to the molecular transport of these quantities, but with modified turbulent transport coefficients. The basis of all these algebraic closures is to model the subgrid scale (sgs) viscosity with:
where \(C_m\) is a model constant, \(\Delta\) is the subgrid length scale (typically the grid size) and \(\mathcal{D}(\boldsymbol{u})\) is a differential operator acting on the filtered velocity. \(\mu_t\) is then calculated using the filtered density \(\mu_t = \overline{\rho} \cdot \nu_{t}\). In the Smagorinsky model, \(C_m = C_s = 0.18\) and \(\mathcal{D}\) writes:
where filtered quantities are indicated with an overbar, Favre-filtered quantities are indicated with a tilde. In the WALE model, \(C_m = C_w = 0.60\) and \(\mathcal{D}\) writes:
When using the Sigma model, \(C_m = C_{\sigma} = 1.35\) and \(\mathcal{D}\) writes:
where the \(\sigma\) are the singular values of the filtered velocity gradient tensor.
The total diffusive transport of momentum from both viscous and turbulent stresses is then computed as
The thermal conducivity and species diffusivities are similarly modified with turbulent contributions, \(\lambda_t = \mu_t \widetilde{c_p} / Pr_t\) and \((\rho D)_t = \mu_t/Sc_t\). The solution algorithm is unchanged other than the addition of these turbulent coefficients to the corresponding molecular transport properties. Nominal values for the model coefficient \(Sc_t = Pr_t = 0.7\).
Limitations: Because the turbulent transport coefficients are nonlinear functions of the velocity field, the treatment of the diffusion terms is not fully implicit when LES models are active. While the implicit solves as described above are kept in place to ensure numerical stability, the turbulent transport coefficients are evaluated only at the old timestep, with the old turbulent values also used to approximate the values at the new timestep. Additionally, the present implementation cannot be used with EFIELD.