$$$$

Advection schemes

In AMReX-Hydro, the fundamental algorithm is either a Method-of-Lines (MOL) or Godunov approach. Each method provides functions for two separate usages:

1. Construct values of the normal velocity at the centroid on each cell face, termed “Pre-MAC” or extrapolated velocity. Typically, this velocity is later MAC projected before being used as the advective velocity. (Information on the MAC projection is in the MAC Projection section.)

2. Construct states on faces, termed “edge states” or “Post-MAC.” These are typically later used to make fluxes which are then differenced to create the advective term. (Information on how fluxes and the convective term are constructed from edge states is given in the Computing Fluxes and Constructing the advective term sections.)

Domain boundary conditions affect the computation of these pre- and post-MAC states in the same way for all advection methods, and this is described in the Boundary conditions section. Note that while these routines ensure that the resulting computed states obey the boundary conditions, they also still require input with filled grow cells. All schemes also use the same routines to construct fluxes and then the convective term.

Next, we provide notation, and then detail the available advection schemes in EB-regular, as well as EB-aware form when available.

Note

If a cell and all of its neighbors have volume fraction of 1 (i.e. they are not cut or covered cells), the EB methodology will return exactly the same answer (to machine precision) as the non-EB methodology. Here we define neighbor to mean any cell that would be involved in the calculation of the face-based state, and the extent of the resulting neighborhood varies depending on the order of the slope used.

Notation

Here we use $(i,j,k)$ to denote cell centers (or centroids for EB), and thus $(i-\frac{1}{2},j,k)$ denotes the lower x-face of the $(i,j,k)$-th cell, $(i,j+\frac{1}{2},k)$ denotes the upper y-face of the $(i,j,k)$-th cell, etc.

Super- or subscript $L$ (for left) indicates a state that has been extrapolated from values at lower x indices. $R$ (for right) indicates extrapolation higher x indices. For example, for the x-face located at $(i+1/2,j,k)$, $L$ indicates extrapolation from the $(i,j,k)$-th cell center/centroid, and $R$ extrapolation from the $(i+1,j,k)$-th cell center/centroid.

Similarly, for the y-dimension, $F$ (for forward) indicates a state that has been extrapolated from values at lower y indices and $B$ (for back) indicates extrapolation from higher y indices. And for the third dimension, $D$ (for down) indicates a state that has been extrapolated from values at lower z indices. $U$ (for up) indicates extrapolation from higher z indices.

${\mathit{U}}^{MAC}$ is the MAC-projected velocity at face centers (or centroids for EB).

We define $\epsilon =1.e-8$ in Utils/hydro_constants.H. This is an empirically determined constant that works well for flows where velocities are on the order of 1.

Method of Lines (MOL)

The procedure for computing MAC velocities and edge states with MOL involves extrapolation in space only, and does not involve any time derivatives. All slope computations use second-order limited slopes as described in Slopes.

These alogrithms are applied in the MOL namespace. For API documentation, see Doxygen: MOL Namespace.

Pre-MAC (API ref. MOL::ExtrapVelToFaces)

For computing the pre-MAC edge states to be MAC-projected, we define on every x-face:

$$\begin{array}{r}\begin{array}{rlr}{u}_L& =& u_{i-1,j,k}+\frac{\mathrm{\Delta }x}{2}{u^x}_{i-1,j,k},\\ u_R& =& u_{i,j,k}-\frac{\mathrm{\Delta }x}{2}{u^x}_{i,j,k},\end{array}\end{array}$$

where $u^x$ are the (limited) slopes in the x-direction.

Boundary conditions are applied (as decribed in Boundary conditions). Then, at each face we upwind based on $u_L$ and $u_R$

$$\begin{array}{r}{u}_{i-\frac{1}{2},j,k}=\{\begin{array}{ll}0,& \mathrm{if}{u}_L<0\mathrm{and}{u}_R>0\mathrm{else}\\ u_L,& \mathrm{if}{u}_L+u_R\ge \epsilon \mathrm{else}\\ u_R,& \mathrm{if}{u}_L+u_R\le -\epsilon \mathrm{else}\\ 0\end{array}\end{array}$$

We similarly compute $v_{i,j-\frac{1}{2},k}$ on y-faces and $w_{i,j,k-\frac{1}{2}}$ on z-faces.

Post-MAC (API ref. MOL::ComputeEdgeState)

Once we have the MAC-projected velocities, we extrapolate all quantities to faces as above:

$$\begin{array}{r}\begin{array}{rlr}{s}_L& =& s^{i-1,j,k}+\frac{\mathrm{\Delta }x}{2}{s^x}_{i-1,j,k},\\ s_R& =& s^{i,j,k}-\frac{\mathrm{\Delta }x}{2}{s^x}_{i,j,k},\end{array}\end{array}$$

where $s^x$ are the (limited) slopes in the x-direction.

Boundary conditions are applied (as decribed in Boundary conditions). Then, at each face, we upwind based on $u_{i-\frac{1}{2},j,k}^{MAC}$

$$\begin{array}{r}{s}_{i-\frac{1}{2},j,k}=\{\begin{array}{ll}{s}_L,& \mathrm{if}{u}_{i-\frac{1}{2},j,k}^{MAC}\ge \epsilon \mathrm{else}\\ s_R,& \mathrm{if}{u}_{i-\frac{1}{2},j,k}^{MAC}\le -\epsilon \mathrm{else}\\ \frac{1}{2}(s_L+s_R),\end{array}\end{array}$$

Method of Lines with Embedded Boundaries (EBMOL)

AMReX-Hydro has also implemented an embedded boundary (EB) aware version of the MOL algorithm discussed above. All slope computations use second-order limited slopes as described in EB Slopes.

Pre-MAC (API ref. EBMOL::ExtrapVelToFaces)

For computing the pre-MAC edge states to be MAC-projected, we define on every x-face with non-zero area fraction:

$$\begin{array}{r}\begin{array}{rlr}{u}_L& =& u_{i-1,j,k}+{\delta }_x{u^x}_{i-1,j,k}+{\delta }_y{u^y}_{i-1,j,k}+\delta z{u^z}_{i-1,j,k},\\ u_R& =& u_{i,j,k}-{\delta }_x{u^x}_{i,j,k}-{\delta }_y{u^y}_{i,j,k}-\delta z{u^z}_{i,j,k},\end{array}\end{array}$$

where we calculate $u^x$, $u^y$ and $u^z$ as described in EB Slopes, and ${\delta }_x$, ${\delta }_y$, and ${\delta }_z$ are the components of the distance vector from the cell centroid to the face centroid of the face at $(i-\frac{1}{2},j,k).$

Boundary conditions are applied (as decribed in Boundary conditions). Then, at each face we upwind based on $u_L$ and $u_R$

$$\begin{array}{r}{u}_{i-\frac{1}{2},j,k}=\{\begin{array}{ll}0,& \mathrm{if}{u}_L<0\mathrm{and}{u}_R>0\mathrm{else}\\ u_L,& \mathrm{if}{u}_L+u_R\ge \epsilon \mathrm{else}\\ u_R,& \mathrm{if}{u}_L+u_R\le -\epsilon \mathrm{else}\\ 0\end{array}\end{array}$$

We similarly compute $v_{i,j-\frac{1}{2},k}$ on y-faces and $w_{i,j,k-\frac{1}{2}}$ on z-faces.

Post-MAC (API ref. EBMOL::ComputeEdgeState)

Once we have the MAC-projected velocities, we predict all quantities to faces with non-zero area fractions as above:

$$\begin{array}{r}\begin{array}{rlr}{s}_L& =& s_{i-1,j,k}+{\delta }_x{s^x}_{i-1,j,k}+{\delta }_y{s^y}_{i-1,j,k}+\delta z{s^z}_{i-1,j,k},\\ s_R& =& s_{i,j,k}-{\delta }_x{s^x}_{i,j,k}-{\delta }_y{s^y}_{i,j,k}-\delta z{s^z}_{i,j,k},\end{array}\end{array}$$

where we calculate $s^x$, $s^y$ and $s^z$ as described in EB Slopes, and ${\delta }_x$, ${\delta }_y$, and ${\delta }_z$ are the components of the distance vector from the cell centroid to the face centroid of the face at $(i-\frac{1}{2},j,k).$

Boundary conditions are applied (as decribed in Boundary conditions). Then, at each face we then upwind based on $u_{i-\frac{1}{2},j,k}^{MAC}$

$$\begin{array}{r}{s}_{i-\frac{1}{2},j,k}=\{\begin{array}{ll}{s}_L,& \mathrm{if}{u}_{i-\frac{1}{2},j,k}^{MAC}\ge \epsilon \mathrm{else}\\ s_R,& \mathrm{if}{u}_{i-\frac{1}{2},j,k}^{MAC}\le -\epsilon \mathrm{else}\\ \frac{1}{2}(s_L+s_R),\end{array}\end{array}$$

$$$$

Godunov Methods

AMReX-Hydro’s implementation uses dimenensionally unsplit algorithms with full corner coupling in 3D, with the option to use either piecewise linear (PLM) [Saltzman, 1994, Colella, 1990] or piecewise parabolic (PPM) [Colella and Woodward, 1984, Miller and Colella, 2002] reconstructions of the state.

These alogrithms are applied in the Godunov namespace. For API documentation, see Doxygen: Godunov Namespace.

Pre-MAC (API ref. Godunov::ExtrapVelToFaces)

We extrapolate the normal velocities to cell faces using a second-order Taylor series expansion in space and time. For each face, we extrapolate the normal velocity component from the centers of the cells on either side to the face, creating left (L) and right (R) states. For face $(i+1/2,j,k)$ this gives

(1)$$\begin{array}{rl}{\tilde{u}}_{i+\frac{1}{2},j,k}^{L,n+\frac{1}{2}}\approx & u_{i,j,k}^n+\frac{dx}{2}{u}_x+\frac{dt}{2}{u}_t\\ =& u_{i,j,k}^n+(\frac{dx}{2}-u_{i,j,k}^n\frac{dt}{2})(u_x^{n,lim}{)}_{i,j,k}\\ & +\frac{dt}{2}(-(\widehat{v{u}_y}{)}_{i,j,k}-(\widehat{w{u}_z}{)}_{i,j,k}+F_{x,i,j,k}^n)\end{array}$$

extrapolated from $(i,j,k)$, and

(2)$$\begin{array}{rl}{\tilde{u}}_{i+\frac{1}{2},j,k}^{R,n+\frac{1}{2}}\approx & u_{i+1,j,k}^n-\frac{dx}{2}{u}_x+\frac{dt}{2}{u}_t\\ =& u_{i+1,j,k}^n-(\frac{dx}{2}+u_{i+1,j,k}^n\frac{dt}{2})(u_x^{n,lim}{)}_{i+1,j,k}\\ & +\frac{dt}{2}(-(\widehat{v{u}_y}{)}_{i+1,j,k}-(\widehat{w{u}_z}{)}_{i+1,j,k}+F_{x,i+1,j,k}^n)\end{array}$$

extrapolated from $(i+1,j,k).$ Here, $\mathit{F}$ is the sum of forcing terms, which typically might include viscous, gravitational, and the pressure gradient terms.

If the parameter use_ppm is false, the first derivatives normal to the face (in this case $u_x^{n,lim}$) are evaluated using a monotonicity-limited fourth-order slope approximation as described in Slopes. Otherwise the Piecewise Parabolic Method (PPM) [Colella and Woodward, 1984, Miller and Colella, 2002] is used to compute

$$\begin{array}{rl}{\hat{u}}_{i+\frac{1}{2},j,k}^L=& u_{i,j,k}^n+(\frac{dx}{2}-u_{i,j,k}^n\frac{dt}{2})(u_x^{n,lim}{)}_{i,j,k}\\ {\hat{u}}_{i+\frac{1}{2},j,k}^R=& u_{i+1,j,k}^n-(\frac{dx}{2}+u_{i+1,j,k}^n\frac{dt}{2})(u_x^{n,lim}{)}_{i+1,j,k}\end{array}$$

The transverse derivative terms ($\widehat{v{u}_y}$ and $\widehat{w{u}_z}$ in this case) are evaluated using a three step process: first extrapolating all velocity components to the transverse faces from the cell centers on either side (either via PLM or PPM), applying boundary conditions and then choosing between these states using the upwinding procedure defined below. In particular, in the $y$ direction we define

$$\begin{array}{r}\begin{array}{rl}{\hat{\mathit{U}}}_{i,j+\frac{1}{2},k}^F=& {\mathit{U}}_{i,j,k}^n+(\frac{dy}{2}-\frac{dt}{2}{v}_{i,j,k}^n)({\mathit{U}}_y^{n,lim}{)}_{i,j,k}\\ {\hat{\mathit{U}}}_{i,j+\frac{1}{2},k}^B=& {\mathit{U}}_{i,j+1,k}^n-(\frac{dy}{2}+\frac{dt}{2}{v}_{i,j+1,k}^n)({\mathit{U}}_y^{n,lim}{)}_{i,j+1,k}\end{array}\end{array}$$

Values are similarly traced from $(i,j,k)$ and $(i,j,k+1)$ to the $(i,j,k+\frac{1}{2})$ faces to define ${\hat{\mathit{U}}}_{i,j,k+\frac{1}{2}}^D$ and ${\hat{\mathit{U}}}_{i,j,k+\frac{1}{2}}^U$, respectively.

If the parameter use_forces_in_trans is true, the forcing terms ($\mathit{F}$ in Eqs. (1) and (2)) are added to $\hat{\mathit{U}}$ now. Otherwise, they are included later in the the formation of the edge state.

Next, boundary conditions are enforced on domain faces as decribed in Boundary conditions #2. Note that this means face-based values lying within the physical boundary but not exactly on the boundary face (e.g. values located on the y-faces of ghost cells abutting the x-boundary but not on the y-boundary) do not have boundary conditions enforced at this point.

Now, we define a normal advective velocity on the face (suppressing the $(i,j+\frac{1}{2},k)$ spatial indices on front and back states here and in the next equation):

$$\begin{array}{r}{\hat{v}}_{i,j+\frac{1}{2},k}^{adv}=\{\begin{array}{ll}{\hat{v}}^F& \text{if }{\hat{v}}^F>0,{\hat{v}}^F+{\hat{v}}^B>\epsilon \\ 0& \text{if }{\hat{v}}^F\le 0,{\hat{v}}^B\ge 0\text{ or }|{\hat{v}}^F+{\hat{v}}^B|<\epsilon \\ {\hat{v}}^B& \text{if }{\hat{v}}^B<0,{\hat{v}}^F+{\hat{v}}^B<\epsilon \text{ .}\end{array}\end{array}$$

We now upwind $\hat{\mathit{U}}$ based on ${\hat{v}}_{i,j+\frac{1}{2},k}^{adv}$:

$$\begin{array}{r}{\hat{\mathit{U}}}_{i,j+\frac{1}{2},k}=\{\begin{array}{ll}{\hat{\mathit{U}}}^F& \text{if }{\hat{v}}_{i,j+\frac{1}{2},k}^{adv}>\epsilon \\ {\hat{\mathit{U}}}^B& \text{if }{\hat{v}}_{i,j+\frac{1}{2},k}^{adv}<-\epsilon \\ \frac{1}{2}({\hat{\mathit{U}}}^F+{\hat{\mathit{U}}}^B)& \text{otherwise}\end{array}\end{array}$$

In 3D, we complete the intermediate transverse-face centered states by accounting for transverse corner coupling following [Saltzman, 1994, Colella, 1990]. For example, for $u$ predicted to y-faces we modify with a factor of the z-derivative at cell centers

$$\begin{array}{rl}{\breve{u}}_{i,j+\frac{1}{2},k}^F=& {\hat{u}}_{i,j+\frac{1}{2},k}^F-\frac{dt}{3}{\left({\hat{w}}^{adv}{u}_z\right)}_{i,j,k}\\ =& {\hat{u}}_{i,j+\frac{1}{2},k}^F-\frac{dt}{3}\left(\frac{{\hat{u}}_{i,j,k+\frac{1}{2}}{\hat{w}}_{i,j,k+\frac{1}{2}}^{adv}-{\hat{u}}_{i,j,k-\frac{1}{2}}{\hat{w}}_{i,j,k-\frac{1}{2}}^{adv}}{dz}\right)\\ & +\frac{dt}{3}\left(\frac{{\hat{w}}_{i,j,k+\frac{1}{2}}^{adv}-{\hat{w}}_{i,j,k-\frac{1}{2}}^{adv}}{dz}\right)u_{i,j,k}\\ {\breve{u}}_{i,j+\frac{1}{2},k}^B=& {\hat{u}}_{i,j+\frac{1}{2},k}^B-\frac{dt}{3}{\left({\hat{w}}^{adv}{u}_z\right)}_{i,j+1,k}\\ =& {\hat{u}}_{i,j+\frac{1}{2},k}^B-\frac{dt}{3}\left(\frac{{\hat{u}}_{i,j+1,k+\frac{1}{2}}{\hat{w}}_{i,j+1,k+\frac{1}{2}}^{adv}-{\hat{u}}_{i,j+1,k-\frac{1}{2}}{\hat{w}}_{i,j+1,k-\frac{1}{2}}^{adv}}{dz}\right)\\ & +\frac{dt}{3}\left(\frac{{\hat{w}}_{i,j+1,k+\frac{1}{2}}^{adv}-{\hat{w}}_{i,j+1,k-\frac{1}{2}}^{adv}}{dz}\right)u_{i,j+1,k}\end{array}$$

and then apply boundary conditions on domian faces before upwinding according to ${\hat{v}}_{i,j+\frac{1}{2},k}^{adv}$ as was done above for $\hat{\mathit{U}}$. ${\stackrel{⌣}{\mathit{U}}}_{i,j-\frac{1}{2},k},{\stackrel{⌣}{\mathit{U}}}_{i,j,k+\frac{1}{2}}$ and ${\stackrel{⌣}{\mathit{U}}}_{i,j,k-\frac{1}{2}}$ are constructed in a similar manner. For 2D, we take $\stackrel{⌣}{\mathit{U}}=\hat{\mathit{U}}$.

We use these upwind values to form the transverse derivatives in Eqs. (1) and (2) :

(3)$$(\widehat{v{u}_y}{)}_{i,j,k}=\frac{1}{2dy}({\hat{v}}_{i,j+\frac{1}{2},k}^{adv}+{\hat{v}}_{i,j-\frac{1}{2},k}^{adv})({\breve{u}}_{i,j+\frac{1}{2},k}-{\breve{u}}_{i,j-\frac{1}{2},k})$$

$$(\widehat{w{u}_z}{)}_{i,j,k}=\frac{1}{2dz}({\hat{w}}_{i,j,k+\frac{1}{2}}^{adv}+{\hat{w}}_{i,j,k-\frac{1}{2}}^{adv})({\breve{u}}_{i,j,k+\frac{1}{2}}-{\breve{u}}_{i,j,k-\frac{1}{2}})$$

We now have all the terms needed to form $\tilde{u}$. If use_forces_in_trans is false, the forcing terms were not included in the computation of the transverse deriviates and are instead included at this point. We apply boundary conditions on domain faces, including preventing backflow (as decribed in Boundary conditions #2 & 3).

The normal velocity at each face is then determined by an upwinding procedure based on the states predicted from the cell centers on either side. The procedure is similar to that described above, i.e. (suppressing the ($i+\frac{1}{2},j,k$) indices)

$$\begin{array}{r}{\tilde{u}}_{i+\frac{1}{2},j,k}^{n+\frac{1}{2}}=\{\begin{array}{ll}{\tilde{u}}^{L,n+\frac{1}{2}}& \text{if }{\tilde{u}}^{L,n+\frac{1}{2}}>0\text{ and }{\tilde{u}}^{L,n+\frac{1}{2}}+{\tilde{u}}^{R,n+\frac{1}{2}}>\epsilon \\ 0& \text{if }{\tilde{u}}^{L,n+\frac{1}{2}}\le 0,{\tilde{u}}^{R,n+\frac{1}{2}}\ge 0\text{ or }|{\tilde{u}}^{L,n+\frac{1}{2}}+{\tilde{u}}^{R,n+\frac{1}{2}}|<\epsilon \\ {\tilde{u}}^{R,n+\frac{1}{2}}& \text{if }{\tilde{u}}^{R,n+\frac{1}{2}}<0\text{ and }{\tilde{u}}^{L,n+\frac{1}{2}}+{\tilde{u}}^{R,n+\frac{1}{2}}<\epsilon \end{array}\end{array}$$

We follow a similar procedure to construct ${\tilde{v}}_{i,j+\frac{1}{2},k}^{n+\frac{1}{2}}$ and ${\tilde{w}}_{i,j,k+\frac{1}{2}}^{n+\frac{1}{2}}$. We refer to this unique value of normal velocity on each face as ${\mathit{U}}^{MAC,\ast }$.

Post-MAC (API ref. Godnuov::ComputeEdgeState)

Here, the face-centered advective velocity field, which we will call ${\mathit{U}}^{MAC}$, is already known. Typically, this is the MAC projection of ${\mathit{U}}^{MAC,\ast }$, which would have been computed via the Pre-MAC procedure detailed above. We extrapolate all quantities to faces as above:

(4)$$\begin{array}{rl}{\tilde{s}}_{i+\frac{1}{2},j,k}^{L,n+\frac{1}{2}}& \approx s_{i,j,k}^n+\frac{dx}{2}{s}_x+\frac{dt}{2}{s}_t\\ & =s_{i,j,k}^n+(\frac{dx}{2}-s_{i,j,k}^n\frac{dt}{2})(s_x^{n,lim}{)}_{i,j,k}\\ & +\frac{dt}{2}(-(\widehat{v{s}_y}{)}_{i,j,k}-(\widehat{w{s}_z}{)}_{i,j,k}+f_{x,i,j,k}^n)\end{array}$$

extrapolated from $(i,j,k)$, and

(5)$$\begin{array}{rl}{\tilde{s}}_{i+\frac{1}{2},j,k}^{R,n+\frac{1}{2}}& \approx s_{i+1,j,k}^n-\frac{dx}{2}{s}_x+\frac{dt}{2}{s}_t\\ & =s_{i+1,j,k}^n-(\frac{dx}{2}+s_{i+1,j,k}^n\frac{dt}{2})(s_x^{n,lim}{)}_{i+1,j,k}\\ & +\frac{dt}{2}(-(\widehat{v{s}_y}{)}_{i+1,j,k}-(\widehat{w{s}_z}{)}_{i+1,j,k}+f_{x,i+1,j,k}^n)\end{array}$$

extrapolated from $(i+1,j,k).$ Details of how these terms are computed are analogous to the Pre-MAC case.

At each face we then upwind based on $u_{i-\frac{1}{2},j,k}^{MAC}$

$$\begin{array}{r}{\tilde{s}}_{i-\frac{1}{2},j,k}^{n+\frac{1}{2}}=\{\begin{array}{ll}{\tilde{s}}_L,& \mathrm{if}{u}_{i-\frac{1}{2},j,k}^{MAC}\ge \epsilon \mathrm{else}\\ {\tilde{s}}_R,& \mathrm{if}{u}_{i-\frac{1}{2},j,k}^{MAC}\le -\epsilon \mathrm{else}\\ \frac{1}{2}({\tilde{s}}_L+{\tilde{s}}_R),\end{array}\end{array}$$

Godunov with Embedded Boundaries (EBGodunov)

AMReX-Hydro contains an embedded boundary (EB) aware version of the Godunov algorithm discussed above, although with fewer options available. This EB implementation employs a piecewise linear method (PLM) [], and always includes any forcing terms after the computation of the transverse terms. EBGodunov attempts to use fourth-order limited slopes wherever possible, as described in EB Slopes.

Pre-MAC (API ref. EBGodunov::ExtrapVelToFaces)

We extrapolate the normal velocities to cell faces using a second-order Taylor series expansion in space and time. For each face with a non-zero area fraction, we extrapolate the normal velocity component from the centroids of the cells on either side to the face centroid, creating left (L) and right (R) states. For face $(i+1/2,j,k)$ this gives

(6)$${\tilde{u}}_{i+\frac{1}{2},j,k}^{L,\frac{1}{2}}={\hat{u}}_{i+\frac{1}{2},j,k}^L+\frac{dt}{2}(-(\widehat{v{u}_y}{)}_{i,j,k}-(\widehat{w{u}_z}{)}_{i,j,k}+F_{x,i,j,k}^n)$$

extrapolated from $(i,j,k)$, where

(7)$${\hat{u}}_{i+\frac{1}{2},j,k}^L=u_{i,j,k}^n+({\delta }_x-\frac{dt}{2}{u}_{i,j,k}^n){u^x}_{i,j,k}+{\delta }_y{u^y}_{i,j,k}+{\delta }_z{u^z}_{i,j,k}$$

and

(8)$${\tilde{u}}_{i+\frac{1}{2},j,k}^{R,\frac{1}{2}}={\hat{u}}_{i+\frac{1}{2},j,k}^R+\frac{dt}{2}(-(\widehat{v{u}_y}{)}_{i+1,j,k}-(\widehat{w{u}_z}{)}_{i+1,j,k}+F_{x,i+1,j,k}^n)$$

extrapolated from $(i+1,j,k),$ where

(9)$${\hat{u}}_{i+\frac{1}{2},j,k}^R=u_{i+1,j,k}^n+({\delta }_x-\frac{dt}{2}{u}_{i,j,k}^n){u^x}_{i+1,j,k}+{\delta }_y{u^y}_{i+1,j,k}+{\delta }_z{u^z}_{i+1,j,k}$$

Here, $F$ is the sum of forcing terms; ${\delta }_x,$ ${\delta }_y$ and ${\delta }_z$ are the components of the distance vector from the cell centroid to the face centroid of the $x$-face at $(i-\frac{1}{2},j,k)$; and the slopes $(u^x,u^y,u^z)$ are calculated as described in the EB Slopes section.

The transverse derivative terms ( $\widehat{v{u}_y}$ and $\widehat{w{u}_z}$ in this case) are evaluated by first extrapolating all velocity components to the face centroids of the transverse faces from the cell centers on either side, applying boundary conditions on domain faces, and then choosing between these states using the upwinding procedure defined below. In particular, in the $y$ direction we define ${\hat{\mathit{U}}}_{i,j+\frac{1}{2},k}^F$ and ${\hat{\mathit{U}}}_{i,j+\frac{1}{2},k}^B$ analogously to how we defined ${\hat{u}}_{i+\frac{1}{2},j,k}^R$ and ${\hat{u}}_{i+\frac{1}{2},j,k}^L$, but here on the y-faces and including all three velocity components. Values are similarly traced from $(i,j,k)$ and $(i,j,k+1)$ to the $(i,j,k+\frac{1}{2})$ faces to define ${\hat{\mathit{U}}}_{i,j,k+\frac{1}{2}}^D$ and ${\hat{\mathit{U}}}_{i,j,k+\frac{1}{2}}^U$, respectively.

In this upwinding procedure we first define a normal advective velocity on the face (suppressing the $(i,j+\frac{1}{2},k)$ spatial indices on front and back states here and in the next equation):

$$\begin{array}{r}{\hat{v}}_{i,j+\frac{1}{2},k}^{adv}=\{\begin{array}{ll}{\hat{v}}^F& \text{if }{\hat{v}}^F>0,{\hat{v}}^F+{\hat{v}}^B>\epsilon \\ 0& \text{if }{\hat{v}}^F\le 0,{\hat{v}}^B\ge 0\text{ or }|{\hat{v}}^F+{\hat{v}}^B|<\epsilon \\ {\hat{v}}^B& \text{if }{\hat{v}}^B<0,{\hat{v}}^F+{\hat{v}}^B<\epsilon \text{ .}\end{array}\end{array}$$

We now upwind $\hat{\mathit{U}}$ based on ${\hat{v}}_{i,j+\frac{1}{2},k}^{adv}$:

$$\begin{array}{r}{\hat{\mathit{U}}}_{i,j+\frac{1}{2},k}=\{\begin{array}{ll}{\hat{\mathit{U}}}^F& \text{if }{\hat{v}}_{i,j+\frac{1}{2},k}^{adv}>0\\ \frac{1}{2}({\hat{\mathit{U}}}^F+{\hat{\mathit{U}}}^B)& \text{if }{\hat{v}}_{i,j+\frac{1}{2},k}^{adv}=0\\ {\hat{\mathit{U}}}^B& \text{if }{\hat{v}}_{i,j+\frac{1}{2},k}^{adv}<0\end{array}\end{array}$$

In 3D, we modify the intermediate transverse-face centered states to accounting for transverse corner coupling. For example, for $u$ predicted to y-faces we add a factor of the z-derivative

$$\begin{array}{rl}{\breve{u}}_{i,j+\frac{1}{2},k}^F=& {\hat{u}}_{i,j+\frac{1}{2},k}^F-\frac{dt}{3}{\left({\hat{w}}^{adv}{u}_z\right)}_{i,j,k}\\ =& {\hat{u}}_{i,j+\frac{1}{2},k}^F-\frac{dt}{3}[\left(\frac{{\alpha }_{i,j,k+\frac{1}{2}}{\hat{u}}_{i,j,k+\frac{1}{2}}{\hat{w}}_{i,j,k+\frac{1}{2}}^{adv}-{\alpha }_{i,j,k-\frac{1}{2}}{\hat{u}}_{i,j,k-\frac{1}{2}}{\hat{w}}_{i,j,k-\frac{1}{2}}^{adv}}{dz{\kappa }_{i,j,k}}\right)\\ & -\left(\frac{{\alpha }_{i,j,k+\frac{1}{2}}{\hat{w}}_{i,j,k+\frac{1}{2}}^{adv}-{\alpha }_{i,j,k-\frac{1}{2}}{\hat{w}}_{i,j,k-\frac{1}{2}}^{adv}}{dz{\kappa }_{i,j,k}}\right)u_{i,j,k}]\\ {\breve{u}}_{i,j+\frac{1}{2},k}^B=& {\hat{u}}_{i,j+\frac{1}{2},k}^B-\frac{dt}{3}{\left({\hat{w}}^{adv}{u}_z\right)}_{i,j+1,k}\\ =& {\hat{u}}_{i,j+\frac{1}{2},k}^B-\frac{dt}{3}[\left(\frac{{\alpha }_{i,j+1,k+\frac{1}{2}}{\hat{u}}_{i,j+1,k+\frac{1}{2}}{\hat{w}}_{i,j+1,k+\frac{1}{2}}^{adv}-{\alpha }_{i,j+1,k-\frac{1}{2}}{\hat{u}}_{i,j+1,k-\frac{1}{2}}{\hat{w}}_{i,j+1,k-\frac{1}{2}}^{adv}}{dz{\kappa }_{i,j+1,k}}\right)\\ & -\left(\frac{{\alpha }_{i,j+1,k+\frac{1}{2}}{\hat{w}}_{i,j+1,k+\frac{1}{2}}^{adv}-{\alpha }_{i,j+1,k-\frac{1}{2}}{\hat{w}}_{i,j+1,k-\frac{1}{2}}^{adv}}{dz{\kappa }_{i,j+1,k}}\right)u_{i,j+1,k}]\end{array}$$

where $alpha$ are cell face area fractions and $kappa$ is the cell volume fraction. Then we apply boundary conditions on domain faces before upwinding according to ${\hat{v}}_{i,j+\frac{1}{2},k}^{adv}$ as was done above for $\hat{\mathit{U}}$. ${\stackrel{⌣}{\mathit{U}}}_{i,j-\frac{1}{2},k},{\stackrel{⌣}{\mathit{U}}}_{i,j,k+\frac{1}{2}}$ and ${\stackrel{⌣}{\mathit{U}}}_{i,j,k-\frac{1}{2}}$ are constructed in a similar manner. For 2D, we take $\stackrel{⌣}{\mathit{U}}=\hat{\mathit{U}}$.

If any of the four faces that contribute to the transverse derivatives for a particular cell have zero area, all of the transverse and forcing terms are identically set to 0. For example, when constructing ${\tilde{u}}_{i+\frac{1}{2},j,k}^{L,n+\frac{1}{2}}$, if any of the areas $a_{i,j+\frac{1}{2},k},a_{i,j-\frac{1}{2},k},a_{i,j,k-\frac{1}{2}}$ or $a_{i,j,k+\frac{1}{2}}$ are zero, then we simply define

(10)$${\tilde{u}}_{i+\frac{1}{2},j,k}^{L,n+\frac{1}{2}}={\hat{u}}_{i+\frac{1}{2},j,k}^L$$

and use this in the upwinding step given by Eq. (11).

For cut faces (i.e. faces intersecting the EB), $\stackrel{⌣}{\mathit{U}}$ and ${\hat{\mathit{U}}}_{ad}$ are linearly extrapolated from face centroids to face centers. Then the transverse derivatives are constructed from these face centered values using the same formulas as for the non-EB case (Eq. (3)).

The normal velocity at each face is then determined by an upwinding procedure based on the states predicted from the cells on either side. The procedure is similar to that described above, i.e. (suppressing the ($i+\frac{1}{2},j,k$) indices)

(11)$$\begin{array}{r}{\tilde{u}}_{i+\frac{1}{2},j,k}^{n+\frac{1}{2}}=\{\begin{array}{ll}{\tilde{u}}^{L,n+\frac{1}{2}}& \text{if }{\tilde{u}}^{L,n+\frac{1}{2}}>0\text{ and }{\tilde{u}}^{L,n+\frac{1}{2}}+{\tilde{u}}^{R,n+\frac{1}{2}}>\epsilon \\ 0& \text{if }{\tilde{u}}^{L,n+\frac{1}{2}}\le 0,{\tilde{u}}^{R,n+\frac{1}{2}}\ge 0\text{ or }|{\tilde{u}}^{L,n+\frac{1}{2}}+{\tilde{u}}^{R,n+\frac{1}{2}}|<\epsilon \\ {\tilde{u}}^{R,n+\frac{1}{2}}& \text{if }{\tilde{u}}^{R,n+\frac{1}{2}}<0\text{ and }{\tilde{u}}^{L,n+\frac{1}{2}}+{\tilde{u}}^{R,n+\frac{1}{2}}<\epsilon \end{array}\end{array}$$

We follow a similar procedure to construct ${\tilde{v}}_{i,j+\frac{1}{2},k}^{n+\frac{1}{2}}$ and ${\tilde{w}}_{i,j,k+\frac{1}{2}}^{n+\frac{1}{2}}$. We refer to these unique values of normal velocity on each face as ${\mathit{U}}^{MAC,\ast }$.

Post-MAC (API ref. EBGondunov::ComputeEdgestate)

Here, the face-centered advective velocity field, which we will call ${\mathit{U}}^{MAC}$, is already known. Typically, this is the MAC projection of ${\mathit{U}}^{MAC,\ast }$, which would have been computed via the Pre-MAC procedure detailed above. Now we predict all quantities to faces. Let the scalar $s$ represent any advected quantities as well as all three velocity components. We now extrapolate $s$ from cell centroids to face centroids as described in Sec. Pre-MAC (API ref. EBGodunov::ExtrapVelToFaces). For example, on face $(i+1/2,j,k)$ we define

(12)$${\tilde{s}}_{i+\frac{1}{2},j,k}^{L,n+\frac{1}{2}}={\hat{s}}_{i+\frac{1}{2},j,k}^L+\frac{dt}{2}(-(\widehat{v{s}_y}{)}_{i,j,k}-(\widehat{w{s}_z}{)}_{i,j,k}+f_{x,i,j,k}^n)$$

extrapolated from $(i,j,k)$, where

(13)$${\hat{s}}_{i+\frac{1}{2},j,k}^L=s_{i,j,k}^n+({\delta }_x-\frac{dt}{2}{u}_{i,j,k}^n){s^x}_{i,j,k}+{\delta }_y{s^y}_{i,j,k}+{\delta }_z{s^z}_{i,j,k}$$

and

$${\tilde{s}}_{i+\frac{1}{2},j,k}^{R,n+\frac{1}{2}}={\hat{s}}_{i+\frac{1}{2},j,k}^R+\frac{dt}{2}(-(\widehat{v{s}_y}{)}_{i+1,j,k}-(\widehat{w{s}_z}{)}_{i+1,j,k}+f_{x,i+1,j,k}^n)$$

extrapolated from $(i+1,j,k),$ where

(14)$${\hat{u}}_{i+\frac{1}{2},j,k}^R=u_{i+1,j,k}^n+({\delta }_x-\frac{dt}{2}{u}_{i,j,k}^n){s^x}_{i+1,j,k}+{\delta }_y{s^y}_{i+1,j,k}+{\delta }_z{s^z}_{i+1,j,k}$$

Here again ${\delta }_x,$ ${\delta }_y$ and ${\delta }_z$ are the components of the distance vector from the cell centroid to the face centroid of the $x$-face at $(i-\frac{1}{2},j,k)$, and the slopes $(u^x,u^y,u^z)$ are calculated as decribded in the EB Slopes section.

The transverse terms are computed exactly as described earlier for the Pre-MAC case, except for the upwinding process. Where we previously used the predicted states themselves to upwind, we now use the component of the advective velocity normal to the face in question.

We upwind ${\tilde{s}}_{i+\frac{1}{2},j,k}^{L,n+\frac{1}{2}}$ and ${\tilde{s}}_{i+\frac{1}{2},j,k}^{L,n+\frac{1}{2}}$ using the normal component of ${\mathit{U}}^{MAC}$ to define ${\tilde{s}}_{i+\frac{1}{2},j,k}^{n+\frac{1}{2}}.$ Again, suppressing the subscripts, we define

(15)$$\begin{array}{r}{\tilde{s}}^{n+\frac{1}{2}}=\{\begin{array}{ll}{\tilde{s}}^{L,n+\frac{1}{2}}& \text{if }{u}^{MAC}>\epsilon \\ {\tilde{s}}^{R,n+\frac{1}{2}}& \text{if }{u}^{MAC}<\epsilon \\ \frac{1}{2}({\tilde{s}}^{L,n+\frac{1}{2}}+{\tilde{s}}^{R,n+\frac{1}{2}})& \text{otherwise}\end{array}\end{array}$$

$$$$

BDS Algorithm

The Bell-Dawson-Shubin (BDS) algorithm is a higher order Godunov method for scalar conservation laws in multiple dimensions. Satisfying the maximum principal for constant coefficient linear advection, the BDS routine provides accurate resolution for smooth problems while avoiding undershoot and overshoot for non-smooth profiles. Additional details and comparisons to other schemes can be found in [Nonaka et al., 2011] and references therein.

This implementation of BDS closely follows the Godunov approach and leverages some of the same code. The difference appears in the computation of edge states given face-centered velocities, i.e. the Post-MAC computation. Currently, periodic, Dirichlet, and outflow (extrapolation) boundary conditions are supported. Embedded boundaries are not supported within BDS at this time. If additional functionality is desired, or if questions remain after reading this guide, further help is available by submitting an issue through Github

Pre-MAC

The BDS routine follows the Godunov PLM method to extrapolate velocities to cell faces, see Godunov Methods: Pre-MAC.

Post-MAC

In the notation below, $s$ is a scalar field of the form $s=s(x,y,z,t)$ and $\mathit{U}=(u,v,w)$ represents a known face-centered velocity field, typically the projected velocity field from the Pre-MAC step (${\mathit{U}}^{MAC}$). $s_{ijk}^n$ represents the average value of $s$ over the cell with index $(ijk)$ at time $t^n$. At each face the normal velocity (e.g., $u_{i+1/2,j,k}$) is assumed constant over the time step.

Obtaining the edge states is a two step process:

• Step 1: Construct a limited piecewise trilinear (bilinear in 2D) representation of the solution in each grid cell of the form,

$$\begin{array}{r}\begin{array}{rcl}{s}_{ijk}(x,y,z)& =& s_{ijk}+s_{x,ijk}\cdot (x-x_i)+s_{y,ijk}\cdot (y-y_j)+s_{z,ijk}\cdot (z-z_k)\\ & & +s_{xy,ijk}\cdot (x-x_i)(y-y_j)+s_{xz,ijk}\cdot (x-x_i)(z-z_k)\\ & & +s_{yz,ijk}\cdot (y-y_j)(z-z_k)+s_{xyz,ijk}\cdot (x-x_i)(y-y_j)(z-z_k).\end{array}\end{array}$$

• Step 2: Construct edge states $s_{i+1/2,j,k}$, etc. by integrating the limited piecewise trilinear (bilinear in 2D) profiles over the space-time region determined by the characteristic domain of dependence of the face. We enforce no inflow at an outflow face as described in the post-MAC Boundary Conditions Section.