\[\newcommand{\half}{\frac{1}{2}} \newcommand{\nph}{{n + \half}} \newcommand{\nmh}{{n - \frac{1}{2}}} \newcommand{\iphj}{{i+\frac{1}{2},j,k}} \newcommand{\ijph}{{i,j+\frac{1}{2}},k} \newcommand{\imhj}{{i-\frac{1}{2},j,k}} \newcommand{\ijmh}{{i,j-\frac{1}{2}},k} \newcommand{\ijkmh}{{i,j,k-\frac{1}{2}}} \newcommand{\ijkph}{{i,j,k+\frac{1}{2}}} \newcommand{\grad}{\nabla} \newcommand{\del}{\nabla} \newcommand{\AN}{[(U \cdot \nabla)U]^{n+\frac{1}{2}}} \newcommand{\npk}{{n + \frac{p+\half}{R}}} \newcommand{\nak}{{n + \frac{p}{R}}} \newcommand{\nmk}{{n + \frac{p-\half}{R}}} \newcommand{\iph}{i+\half} \newcommand{\imh}{i-\half} \newcommand{\ipmh}{i\pm\half} \newcommand{\jph}{j+\half} \newcommand{\jmh}{j-\half} \newcommand{\jpmh}{j\pm\half} \newcommand{\kph}{k+\half} \newcommand{\kmh}{k-\half} \newcommand{\GMAC}{C \rightarrow E} \newcommand{\DMAC}{E \rightarrow C} \newcommand{\U}{\boldsymbol{U}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\dt}{\Delta t} \newcommand{\shalf}{\sfrac{1}{2}} \newcommand\hathat[1]{\widehat{\widehat{#1}}} \newcommand\hairbow[1]{\overset{\bowtie}{#1}} \newcommand\widebreve[1]{\overset{\smile}{#1}} \newcommand{\vol}{\mathcal{V}} \newcommand{\area}{\mathcal{A}}\]

Helper functions

Notation

  • \(dx\), \(dy\), \(dz\) : cell sizes in the x-, y- and z- directions, respectively.

  • \(\vol_{i,j,k}\) : Volume of cell-centered element \((i,j,k)\).

  • \(\area\) : Area of a cell face. For example, \(\area_{i-\half ,j,k}\) represents the area of the lower x-face of the \((i,j,k)\)-th cell,

For problems with embedded boundaries, we also define

  • \(\kappa_{i,j,k}\) : Volume fraction of cell \((i,j,k)\). Data are in the range of \([0,1]\) with zero representing covered cells and one for regular cells, so that for regular cells \(\kappa_{i,j,k} \vol_{i,j,k} = dx\ dy\ dz\)

  • \(\alpha\) : Area fractions. For example, \(\alpha_{i-\half ,j,k}\) corresponds to the the lower x-face of the \((i,j,k)\)-th cell. Data are in the range of \([0,1]\) with zero representing a covered face and one an un-cut face, so that for regular cells \(\alpha_{i-\half ,j,k} \area_{i-\half ,j,k} = dy\ dz\).

Note that EB is only an option for Cartesian geometries as of this writing.

Slopes

AMReX-Hydro’s implementation of the piecewise linear method provides several options, and leverages slopes routines from AMReX. For cells where this calculation would involve all regular cells (i.e. no cut or covered cells), there are second-order and fourth-order stencils along with options to apply limiters. Note that the piecewise parabolic and BDS methods have their own routines for formulating slopes (and these are housed within AMReX-Hydro).

For (EB)Godunov, the default is monotonicity-limited fourth-order slopes. Default limiting is as described in Colella (1985) [Colella and Glaz, 1985], where limiting is done on each component of the velocity individually.

For (EB)MOL, the default is monotonicity-limited second-order slopes. Default is the second order Monotonized Central (MC) limiter (van Leer, 1977 [Van Leer, 1977]), where limiting is applied direction by direction.

NOTE ON BOUNDARY CONDITIONS:

When periodic or Neumann BCs are imposed, schemes can be applied without any change since the ghost cells outside the domain are filled by either periodicity or by extrapolation.

For Dirichlet BCs, the BC value stored in the first ghost cell outside the domain is considered as located directly on the boundary, despite the fact that it is stored in what is otherwise considered a cell-centered array. We then utilize one-sided differencing schemes that use ONLY values from inside or on the domain boundary.

EB Slopes

The procedure for problems with embedded boundaries is detailed below, and attempts to use standard (non-EB) stencils wherever possible.

  • First, AMReX-Hydro attempts to compute the slope of the desired order (4 for Godunov and 2 for MOL) using standard stencils and limiters.

  • For cases where a fourth-order slope is desired but the stencil would require the use of cut cells (as happens in EBGodunov), the code next attempts to use the standard second-order slope and limiter.

  • If the standard second-order slope calculation would require the use of cut cells, then the slope computation will use a least squares approach, involving a linear fit to the at-most 26 (or 8 in 2D) nearest neighbors, with the function going through the centroid of cell (i,j,k) to the face centroid. This does not assume that the cell centroids, where the data is assumed to live, are the same as cell centers. This least-squares slope is then multiplied by a limiter based on the work of Barth-Jespersen that enforces no new maxima or minima when the state is predicted to the face centroids.

\[\newcommand{\half}{\frac{1}{2}} \newcommand{\nph}{{n + \half}} \newcommand{\nmh}{{n - \frac{1}{2}}} \newcommand{\iphj}{{i+\frac{1}{2},j,k}} \newcommand{\ijph}{{i,j+\frac{1}{2}},k} \newcommand{\imhj}{{i-\frac{1}{2},j,k}} \newcommand{\ijmh}{{i,j-\frac{1}{2}},k} \newcommand{\ijkmh}{{i,j,k-\frac{1}{2}}} \newcommand{\ijkph}{{i,j,k+\frac{1}{2}}} \newcommand{\grad}{\nabla} \newcommand{\del}{\nabla} \newcommand{\AN}{[(U \cdot \nabla)U]^{n+\frac{1}{2}}} \newcommand{\npk}{{n + \frac{p+\half}{R}}} \newcommand{\nak}{{n + \frac{p}{R}}} \newcommand{\nmk}{{n + \frac{p-\half}{R}}} \newcommand{\iph}{i+\half} \newcommand{\imh}{i-\half} \newcommand{\ipmh}{i\pm\half} \newcommand{\jph}{j+\half} \newcommand{\jmh}{j-\half} \newcommand{\jpmh}{j\pm\half} \newcommand{\kph}{k+\half} \newcommand{\kmh}{k-\half} \newcommand{\GMAC}{C \rightarrow E} \newcommand{\DMAC}{E \rightarrow C} \newcommand{\U}{\boldsymbol{U}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\dt}{\Delta t} \newcommand{\shalf}{\sfrac{1}{2}} \newcommand\hathat[1]{\widehat{\widehat{#1}}} \newcommand\hairbow[1]{\overset{\bowtie}{#1}} \newcommand\widebreve[1]{\overset{\smile}{#1}} \newcommand{\vol}{\mathcal{V}} \newcommand{\area}{\mathcal{A}}\]

Computing Fluxes

Doxygen links ComputeFluxes and EB_ComputeFluxes .

AMReX-Hydro has the option to compute intesive or extensive, i.e. area-weighted, fluxes. Extensive fluxes are always used for problems using R-Z geometry, and AMReX has functions for computing the position dependent cell volume and face areas (see Geometry::GetFaceArea and Geometry::GetVolume). We first give the formulas for the EB-regular case, and then those used when embedded boundaries are present.

Intensive

Intensive fluxes are computed from the advective velocity, \(\U^{MAC}\), and edge states, \(s\), by defining

(16)\[F_{i-\frac{1}{2},j,k}^{x,n+\frac{1}{2}} = u^{MAC}_{i-\frac{1}{2},j,k}\; s_{i-\frac{1}{2},j,k}^{n+\frac{1}{2}}\]

on all x-faces,

(17)\[F_{i,j-\frac{1}{2},k}^{y,n+\frac{1}{2}} = v^{MAC}_{i,j-\frac{1}{2},k}\; s_{i,j-\frac{1}{2},k}^{n+\frac{1}{2}}\]

on all y-faces, and

(18)\[F_{i,j,k-\frac{1}{2}}^{z,n+\frac{1}{2}} = w^{MAC}_{i,j,k-\frac{1}{2}}\; s_{i,j,k-\frac{1}{2}}^{n+\frac{1}{2}}\]

on all z-faces.


When embedded boundaries are present, intensive fluxes are computed as

\[F_{i-\frac{1}{2},j,k}^{x,n+\frac{1}{2}} = \alpha_{i-\frac{1}{2},j,k} \; u^{MAC}_{i-\frac{1}{2},j,k} \; s_{i-\frac{1}{2},j,k}^{n+\frac{1}{2}}\]

on all x-faces,

\[F_{i,j-\frac{1}{2},k}^{y,n+\frac{1}{2}} = \alpha_{i,j-\frac{1}{2},k} \; v^{MAC}_{i,j-\frac{1}{2},k} \; s_{i,j-\frac{1}{2},k}^{n+\frac{1}{2}}\]

on all y-faces,

\[F_{i,j,k-\frac{1}{2}}^{z,n+\frac{1}{2}} = \alpha_{i,j,k-\frac{1}{2}} \; w^{MAC}_{i,j,k-\frac{1}{2}}\; s_{i,j,k-\frac{1}{2}}^{n+\frac{1}{2}}\]

on all z-faces. Here \(\alpha_{i-\frac{1}{2},j,k}\) is the area fraction of the lower x-face of cell-centered element \((i,j,k)\), etc.

Extensive

Extensive fluxes are computed as

\[F_{i-\frac{1}{2},j,k}^{x,n+\frac{1}{2}} = \area_{i-\frac{1}{2},j,k} \; u^{MAC}_{i-\frac{1}{2},j,k} \; s_{i-\frac{1}{2},j,k}^{n+\frac{1}{2}}\]

on all x-faces,

\[F_{i,j-\frac{1}{2},k}^{y,n+\frac{1}{2}} = \area_{i,j-\frac{1}{2},k} \; v^{MAC}_{i,j-\frac{1}{2},k} \; s_{i,j-\frac{1}{2},k}^{n+\frac{1}{2}}\]

on all y-faces,

\[F_{i,j,k-\frac{1}{2}}^{z,n+\frac{1}{2}} = \area_{i,j,k-\frac{1}{2}} \; w^{MAC}_{i,j,k-\frac{1}{2}}\; s_{i,j,k-\frac{1}{2}}^{n+\frac{1}{2}}\]

on all z-faces.

where \(\area_{i-\frac{1}{2},j,k}\) is the area of the lower x-face of cell-centered element \((i,j,k)\), etc.


For EB, we simply scale area the by the area fraction in the above equations. For example, we use \(\alpha_{i-\frac{1}{2},j,k} \area_{i-\frac{1}{2},j,k}\) in place of \(\area_{i-\frac{1}{2},j,k}\), etc.

\[\newcommand{\half}{\frac{1}{2}} \newcommand{\nph}{{n + \half}} \newcommand{\nmh}{{n - \frac{1}{2}}} \newcommand{\iphj}{{i+\frac{1}{2},j,k}} \newcommand{\ijph}{{i,j+\frac{1}{2}},k} \newcommand{\imhj}{{i-\frac{1}{2},j,k}} \newcommand{\ijmh}{{i,j-\frac{1}{2}},k} \newcommand{\ijkmh}{{i,j,k-\frac{1}{2}}} \newcommand{\ijkph}{{i,j,k+\frac{1}{2}}} \newcommand{\grad}{\nabla} \newcommand{\del}{\nabla} \newcommand{\AN}{[(U \cdot \nabla)U]^{n+\frac{1}{2}}} \newcommand{\npk}{{n + \frac{p+\half}{R}}} \newcommand{\nak}{{n + \frac{p}{R}}} \newcommand{\nmk}{{n + \frac{p-\half}{R}}} \newcommand{\iph}{i+\half} \newcommand{\imh}{i-\half} \newcommand{\ipmh}{i\pm\half} \newcommand{\jph}{j+\half} \newcommand{\jmh}{j-\half} \newcommand{\jpmh}{j\pm\half} \newcommand{\kph}{k+\half} \newcommand{\kmh}{k-\half} \newcommand{\GMAC}{C \rightarrow E} \newcommand{\DMAC}{E \rightarrow C} \newcommand{\U}{\boldsymbol{U}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\dt}{\Delta t} \newcommand{\shalf}{\sfrac{1}{2}} \newcommand\hathat[1]{\widehat{\widehat{#1}}} \newcommand\hairbow[1]{\overset{\bowtie}{#1}} \newcommand\widebreve[1]{\overset{\smile}{#1}} \newcommand{\vol}{\mathcal{V}} \newcommand{\area}{\mathcal{A}}\]

Constructing the advective term

AMReX-Hydro provides the option to compute the advective term either in conservative (\(\grad \cdot \U s\)) or convective form (\(\U \cdot \grad s\)). The exact equations differ slightly depending on whether the fluxes are intensive or extensive (required for R-Z coordinates).

Intensive fluxes

If the variable, \(s\) is to be updated conservatively, we construct

\[\begin{split}\nabla \cdot \left({\bf u} s\right)^{n+\frac{1}{2}} = & \frac{1}{dx} \left(F_{i+\frac{1}{2},j,k}^{x,n+\frac{1}{2}} - F_{i-\frac{1}{2},j,k}^{x,n+\frac{1}{2}}\right) + \\ & \frac{1}{dy} \left(F_{i,j+\frac{1}{2},k}^{y,n+\frac{1}{2}} - F_{i,j-\frac{1}{2},k}^{y,n+\frac{1}{2}}\right) + \\ & \frac{1}{dz} \left(F_{i,j,k+\frac{1}{2}}^{z,n+\frac{1}{2}} - F_{i,j,k-\frac{1}{2}}^{z,n+\frac{1}{2}}\right)\end{split}\]

while if \(s\) is to be updated in convective form, we construct

\[\left({\bf u}\cdot \nabla s\right)^{n+\frac{1}{2}} = \nabla \cdot \left({\bf u} s\right)^{n+\frac{1}{2}} - s_{i,j,k}^{n+\frac{1}{2}} \; \left(DU\right)^{MAC}\]

where

\[\begin{split}\left(DU\right)^{MAC} = \; & \frac{1}{dx} \left(u^{MAC}_{i+\frac{1}{2},j,k} - u^{MAC}_{i-\frac{1}{2},j,k}\right) + \\ & \frac{1}{dy} \left(v^{MAC}_{i,j-\frac{1}{2},k} - v^{MAC}_{i,j-\frac{1}{2},k}\right) + \\ & \frac{1}{dz} \left(w^{MAC}_{i,j,k-\frac{1}{2}} - w^{MAC}_{i,j,k-\frac{1}{2}}\right)\end{split}\]

and

\[s_{i,j,k}^{{n+\frac{1}{2}}} = \frac{1}{6} \left( s_{i-\frac{1}{2},j,k}^{{n+\frac{1}{2}}} + s_{i+\frac{1}{2},j,k}^{{n+\frac{1}{2}}} + s_{i,j-\frac{1}{2},k}^{{n+\frac{1}{2}}} + s_{i,j-\frac{1}{2},k}^{{n+\frac{1}{2}}} + s_{i,j,k-\frac{1}{2}}^{{n+\frac{1}{2}}} + s_{i,j,k-\frac{1}{2}}^{{n+\frac{1}{2}}} \right)\]

Extensive fluxes

Given extensive fluxes we must divide by the cell volume in constructing the advective term. Here, we give the form for EB and note that the EB-regular case is obtained by setting both \(\alpha = 1\) and \(\kappa=1\). If the variable, \(s\) is to be updated conservatively, on all cells with \(\kappa_{i,j,k} > 0\) we construct

\[\begin{split}\nabla \cdot ({\bf u}s)^{n+\frac{1}{2}} = \frac{1}{ \kappa_{i,j,k} \vol_{i,j,k}}[ & ( F_{i+\frac{1}{2},j,k}^{{x,n+\frac{1}{2}}} -F_{i-\frac{1}{2},j,k}^{{x,n+\frac{1}{2}}}) + \\ & ( F_{i,j+\frac{1}{2},k}^{{y,n+\frac{1}{2}}} -F_{i,j-\frac{1}{2},k}^{{y,n+\frac{1}{2}}}) + \\ & ( F_{i,j,k+\frac{1}{2}}^{{z,n+\frac{1}{2}}} -F_{i,j,k-\frac{1}{2}}^{{z,n+\frac{1}{2}}}) ]\end{split}\]

while if \(s\) is to be updated in convective form, we construct

\[({\bf u}\cdot \nabla s)^{n+\frac{1}{2}} = \nabla \cdot ({\bf u}s)^{n+\frac{1}{2}} - s_{i,j,k}^{{n+\frac{1}{2}}} (DU)^{MAC}\]

where

\[\begin{split}(DU)^{MAC} = \frac{1}{\kappa_{i,j,k} \vol_{i,j,k}} [ & (\alpha_{i+\frac{1}{2},j,k} \area_{i+\frac{1}{2},j,k} u^{MAC}_{i+\frac{1}{2},j,k}- \alpha_{i-\frac{1}{2},j,k} \area_{i-\frac{1}{2},j,k} u^{MAC}_{i-\frac{1}{2},j,k}) + \\ & (\alpha_{i,j+\frac{1}{2},k} \area_{i,j+\frac{1}{2},k} v^{MAC}_{i,j-\frac{1}{2},k}- \alpha_{i,j-\frac{1}{2},k} \area_{i,j-\frac{1}{2},k} v^{MAC}_{i,j-\frac{1}{2},k}) + \\ & (\alpha_{i,j,k+\frac{1}{2}} \area_{i,j,k+\frac{1}{2}} w^{MAC}_{i,j,k-\frac{1}{2}}- \alpha_{i,j,k-\frac{1}{2}} \area_{i,j,k-\frac{1}{2}} w^{MAC}_{i,j,k-\frac{1}{2}}) ]\end{split}\]
\[s_{i,j,k}^{{n+\frac{1}{2}}} = (1/6) ( s_{i-\frac{1}{2},j,k}^{{n+\frac{1}{2}}} + s_{i+\frac{1}{2},j,k}^{{n+\frac{1}{2}}} + s_{i,j-\frac{1}{2},k}^{{n+\frac{1}{2}}} + s_{i,j-\frac{1}{2},k}^{{n+\frac{1}{2}}} + s_{i,j,k-\frac{1}{2}}^{{n+\frac{1}{2}}} + s_{i,j,k-\frac{1}{2}}^{{n+\frac{1}{2}}} )\]
\[\newcommand{\half}{\frac{1}{2}} \newcommand{\nph}{{n + \half}} \newcommand{\nmh}{{n - \frac{1}{2}}} \newcommand{\iphj}{{i+\frac{1}{2},j,k}} \newcommand{\ijph}{{i,j+\frac{1}{2}},k} \newcommand{\imhj}{{i-\frac{1}{2},j,k}} \newcommand{\ijmh}{{i,j-\frac{1}{2}},k} \newcommand{\ijkmh}{{i,j,k-\frac{1}{2}}} \newcommand{\ijkph}{{i,j,k+\frac{1}{2}}} \newcommand{\grad}{\nabla} \newcommand{\del}{\nabla} \newcommand{\AN}{[(U \cdot \nabla)U]^{n+\frac{1}{2}}} \newcommand{\npk}{{n + \frac{p+\half}{R}}} \newcommand{\nak}{{n + \frac{p}{R}}} \newcommand{\nmk}{{n + \frac{p-\half}{R}}} \newcommand{\iph}{i+\half} \newcommand{\imh}{i-\half} \newcommand{\ipmh}{i\pm\half} \newcommand{\jph}{j+\half} \newcommand{\jmh}{j-\half} \newcommand{\jpmh}{j\pm\half} \newcommand{\kph}{k+\half} \newcommand{\kmh}{k-\half} \newcommand{\GMAC}{C \rightarrow E} \newcommand{\DMAC}{E \rightarrow C} \newcommand{\U}{\boldsymbol{U}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\dt}{\Delta t} \newcommand{\shalf}{\sfrac{1}{2}} \newcommand\hathat[1]{\widehat{\widehat{#1}}} \newcommand\hairbow[1]{\overset{\bowtie}{#1}} \newcommand\widebreve[1]{\overset{\smile}{#1}} \newcommand{\vol}{\mathcal{V}} \newcommand{\area}{\mathcal{A}}\]

Enforcing inflow-outflow solvability

This routine enforces solvability for inflow-outflow boundaries, which have both inflow and outflow cells. A flux correction factor, \(\alpha_\text{fcf}\) is introduced which scales the outflow velocities to match the inflow:

\[\sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area} + \alpha_\text{fcf} \sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area} = 0\]

The new flux-conserving velocities to be used for the MAC/nodal projections, \({\bf u}_\text{fc}\), are calculated from the correction factor as:

\[\alpha_\text{fcf} = \frac{-\sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area}}{\sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area}},\]
\[{\bf u}_\text{fc} = \alpha_\text{fcf} \cdot {\bf u}, \ \forall {\bf x} \in \partial\Omega_\text{out}.\]

It must be noted that this routine currently only accounts for boundaries with the math BC BCType::direction_dependent, which is to be used for an inflow-outflow boundary. It does not compute or correct the boundary velocities over other math BC types such as those representing pure inflow or pure outflow.