Godunov Interface States

These are working notes for the Godunov step in MAESTROeX and VARDEN.

MAESTROeX Notation

• For 2D, $\mathbf{U}=(u,w)$ and $\tilde{\mathbf{U}}=(\tilde{u},\tilde{w})$. Note that $u=\tilde{u}$. We will use the shorthand $\mathbf{i}=(x,r)$.

• For 3D plane parallel, $\mathbf{U}=(u,v,w)$ and $\tilde{\mathbf{U}}=(\tilde{u},\tilde{v},\tilde{w})$. Note that $u=\tilde{u}$ and $v=\tilde{v}$. We will use the shorthand $\mathbf{i}=(x,y,r)$.

• For 3D spherical, $\mathbf{U}=(u,v,w)$ and $\tilde{\mathbf{U}}=(\tilde{u},\tilde{v},\tilde{w})$. We will use the shorthand $\mathbf{i}=(x,y,z)$.

Computing $\mathbf{U}$ From $\tilde{\mathbf{U}}$

For plane-parallel problems, in order to compute $w$ from $\tilde{w}$, we use simple averaging

(163)$$w_{\mathbf{i}}^n={\tilde{w}}_{\mathbf{i}}^n+\frac{w_{0,r-\frac{1}{2}}+w_{0,r+\frac{1}{2}}}{2}.$$

For spherial problems, in order to compute $\mathbf{U}$ from $\tilde{\mathbf{U}}$, we first map $w_0$ onto $w_0^{\mathrm{M}\mathrm{A}\mathrm{C}}$ using put_w0_on_edges, where $w_0^{\mathrm{M}\mathrm{A}\mathrm{C}}$ only contains normal velocities at each face. Then we construct $\mathbf{U}$ by using

(164)$$u_{\mathbf{i}}={\tilde{u}}_{\mathbf{i}}+\frac{w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}}{2},$$

(165)$$v_{\mathbf{i}}={\tilde{v}}_{\mathbf{i}}+\frac{w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}}{2},$$

(166)$$w_{\mathbf{i}}={\tilde{w}}_{\mathbf{i}}+\frac{w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}}}{2}.$$

To compute full edge-state velocities, simply add $w_0$ (for plane-parallel) or w0mac to the perturbational velocity directly since only edge-based quantities are involved.

Computing $\partial w_0/\partial r$

For plane-parallel problems, the spatial derivatives of $w_0$ are given by the two-point centered difference:

(167)$${\left(\frac{\partial w_0}{\partial r}\right)}_{\mathbf{i}}=\frac{w_{0,r+\frac{1}{2}}-w_{0,r-\frac{1}{2}}}{h}.$$

For spherical problems, we compute the radial bin centered gradient using

(168)$${\left(\frac{\partial w_0}{\partial r}\right)}_r=\frac{w_{0,r+\frac{1}{2}}-w_{0,r-\frac{1}{2}}}{\mathrm{\Delta }r}.$$

Then we put $\partial w_0/\partial r$ onto a Cartesian grid using put_1d_array_on_cart_3d_sphr.

Computing ${\tilde{\mathbf{U}}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$ in MAESTROeX

In advance_premac, we call mkutrans, to compute ${\tilde{\mathbf{U}}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$. We will only compute the normal component of velocity at each face. These transverse velocities do not contain $w_0$, so immediately following the call to mkutrans, we call addw0 to compute ${\mathbf{U}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$ from ${\tilde{\mathbf{U}}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$.

The evolution equation for the perturbational velocity is:

(169)$$\frac{\partial \tilde{\mathbf{U}}}{\partial t}=-\mathbf{U}\cdot \mathrm{\nabla }\tilde{\mathbf{U}}\underset{\text{forcing terms}}{\underbrace{-(\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r)\frac{\partial w_0}{\partial r}{\mathbf{e}}_r-\frac{1}{\rho }\mathrm{\nabla }\pi +\frac{1}{{\rho }_0}\frac{\partial {\pi }_0}{\partial r}{\mathbf{e}}_r-\frac{(\rho -{\rho }_0)}{\rho }g{\mathbf{e}}_r}}.$$

We extrapolate each velocity component to edge-centered, time-centered locations. For example,

(170)$$\begin{array}{r}\begin{array}{rlr}{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& =& {\tilde{u}}_{\mathbf{i}}^n+\frac{h}{2}\frac{\partial {\tilde{u}}_{\mathbf{i}}^n}{\partial x}+\frac{\mathrm{\Delta }t}{2}\frac{\partial {\tilde{u}}_{\mathbf{i}}^n}{\partial t}\\ & =& {\tilde{u}}_{\mathbf{i}}^n+\frac{h}{2}\frac{\partial {\tilde{u}}_{\mathbf{i}}^n}{\partial x}+\frac{\mathrm{\Delta }t}{2}(-{\tilde{u}}_{\mathbf{i}}^n\frac{\partial {\tilde{u}}_{\mathbf{i}}^n}{\partial x}-{\tilde{w}}_{\mathbf{i}}^n\frac{\partial {\tilde{u}}_{\mathbf{i}}^n}{\partial r}+\text{forcing terms})\end{array}\end{array}$$

We are going to use a 1D Taylor series extrapolation in space and time. By 1D, we mean that we omit any spatial derivatives that are not in the direction of the extrapolation. We also omit the underbraced forcing terms. We also use characteristic tracing.

(171)$${\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}={\tilde{u}}_{\mathbf{i}}^n+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}min(0,{\tilde{u}}_{\mathbf{i}}^n)]\partial {\tilde{u}}_{\mathbf{i}}^n$$

2D Cartesian Case

We predict $\tilde{u}$ to x-faces using a 1D extrapolation:

(172)$$\begin{array}{r}\begin{array}{rlr}{\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& =& {\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x}^n+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}max(0,u_{\mathbf{i}-{\mathbf{e}}_x}^n)]{\mathrm{\Delta }}_x{\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x}^n,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& =& {\tilde{u}}_{\mathbf{i}}^n-[\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}min(0,u_{\mathbf{i}}^n)]{\mathrm{\Delta }}_x{\tilde{u}}_{\mathbf{i}}^n.\end{array}\end{array}$$

We pick the final trans states using a Riemann solver:

(173)$$\begin{array}{r}{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}=\{\begin{array}{ll}0,& ({\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}\le 0\mathrm{A}\mathrm{N}\mathrm{D}{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}\ge 0)\mathrm{O}\mathrm{R}|{\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}+{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}|<ϵ,\\ {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}+{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}>0,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}+{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}<0,\end{array}\end{array}$$

We predict $\tilde{w}$ to r-faces using a 1D extrapolation:

(174)$$\begin{array}{r}\begin{array}{rlr}{\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& =& {\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_r}^n+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}max(0,w_{\mathbf{i}-{\mathbf{e}}_r}^n)]{\mathrm{\Delta }}_r{\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_r}^n,\\ {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& =& {\tilde{w}}_{\mathbf{i}}^n-[\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}min(0,w_{\mathbf{i}}^n)]{\mathrm{\Delta }}_r{\tilde{w}}_{\mathbf{i}}^n.\end{array}\end{array}$$

We pick the final $\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}$ states using a Riemann solver, noting that we upwind based on the full velocity.

(175)$$\begin{array}{r}{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}=\{\begin{array}{ll}0,& (w_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}\le 0\mathrm{A}\mathrm{N}\mathrm{D}{w}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}\ge 0)\mathrm{O}\mathrm{R}|w_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}+w_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}|<ϵ,\\ {\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}+w_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}>0,\\ {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}+w_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}<0,\end{array}\end{array}$$

3D Cartesian Case

We use the exact same procedure in 2D and 3D to compute ${\tilde{u}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$ and ${\tilde{w}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$. The procedure for computing ${\tilde{v}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$ is analogous to computing ${\tilde{u}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$. We predict $\tilde{v}$ to y-faces using the 1D extrapolation:

(176)$$\begin{array}{r}\begin{array}{rlr}{\tilde{v}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}& =& {\tilde{v}}_{\mathbf{i}-{\mathbf{e}}_y}^n+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}max(0,v_{\mathbf{i}-{\mathbf{e}}_y}^n)]{\mathrm{\Delta }}_y{\tilde{v}}_{\mathbf{i}-{\mathbf{e}}_y}^n,\\ {\tilde{v}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}& =& {\tilde{v}}_{\mathbf{i}}^n-[\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}min(0,v_{\mathbf{i}}^n)]{\mathrm{\Delta }}_y{\tilde{v}}_{\mathbf{i}}^n,\end{array}\end{array}$$

(177)$$\begin{array}{r}{\tilde{v}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}=\{\begin{array}{ll}0,& (v_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}\le 0\mathrm{A}\mathrm{N}\mathrm{D}{v}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}\ge 0)\mathrm{O}\mathrm{R}|v_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}+v_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}|<ϵ,\\ {\tilde{v}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y},& v_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}+v_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}>0,\\ {\tilde{v}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y},& v_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}+v_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}<0.\end{array}\end{array}$$

3D Spherical Case

We predict the normal components of velocity to the normal faces using a 1D extrapolation. The equations for all three directions are identical to those given in the 2D and 3D plane-parallel sections. As in the plane-parallel case, make sure that the advection velocities, as well as the upwind velocity, is done with the full velocity, not the perturbational velocity.

Computing ${\tilde{\mathbf{U}}}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }$ in MAESTROeX

In advance_premac, we call velpred to compute ${\tilde{\mathbf{U}}}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }$. We will only compute the normal component of velocity at each face.

For reference, here is the perturbational velocity equation from before:

(178)$$\frac{\partial \tilde{\mathbf{U}}}{\partial t}=-\mathbf{U}\cdot \mathrm{\nabla }\tilde{\mathbf{U}}\underset{\text{forcing terms}}{\underbrace{-(\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r)\frac{\partial w_0}{\partial r}{\mathbf{e}}_r\underset{\text{terms included in }{\mathbf{f}}_{\tilde{\mathbf{U}}}}{\underbrace{-\frac{1}{\rho }\mathrm{\nabla }\pi +\frac{1}{{\rho }_0}\frac{\partial {\pi }_0}{\partial r}{\mathbf{e}}_r-\frac{(\rho -{\rho }_0)}{\rho }g{\mathbf{e}}_r}}}}.$$

Note that the $\partial w_0/\partial r$ term is treated like a forcing term, but it is not actually part of ${\mathbf{f}}_{\tilde{\mathbf{U}}}$. We make use of the 1D extrapolations used to compute ${\tilde{\mathbf{U}}}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$ (${\tilde{u}}_{L/R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}$, ${\tilde{v}}_{L/R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}$, and ${\tilde{w}}_{L/R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}$), as well as the “$\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}$” states (${\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$, ${\tilde{v}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$, and ${\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$)

2D Cartesian Case

1. Predict $\tilde{u}$ to r-faces using a 1D extrapolation.

2. Predict $\tilde{u}$ to x-faces using a full-dimensional extrapolation.

3. Predict $\tilde{w}$ to x-faces using a 1D extrapolation.

4. Predict $\tilde{w}$ to r-faces using a full-dimensional extrapolation.

Predict $\tilde{u}$ to r-faces using a 1D extrapolation:

(179)$$\begin{array}{r}\begin{array}{rlr}{\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& =& {\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_r}^n+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}max(0,w_{\mathbf{i}-{\mathbf{e}}_r}^n)]{\mathrm{\Delta }}_r{\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_r}^n,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& =& {\tilde{u}}_{\mathbf{i}}-[\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}min(0,w_{\mathbf{i}}^n)]{\mathrm{\Delta }}_r{\tilde{u}}_{\mathbf{i}}^n.\end{array}\end{array}$$

Upwind based on $w^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$:

(180)$$\begin{array}{r}{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}=\{\begin{array}{ll}\frac{1}{2}({\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}+{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}),& \left|w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}\right|<ϵ\\ {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}>0,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}<0.\end{array}\end{array}$$

Predict $\tilde{u}$ to x-faces using a full-dimensional extrapolation,

(181)$$\begin{array}{r}\begin{array}{rlr}{\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }& =& {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}-\frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+w_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}-{\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r})+\frac{\mathrm{\Delta }t}{2}{f}_{\tilde{u},\mathbf{i}-{\mathbf{e}}_x},\\ & & \\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }& =& {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}-\frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})+\frac{\mathrm{\Delta }t}{2}{f}_{\tilde{u},\mathbf{i}}.\end{array}\end{array}$$

Solve a Riemann problem:

(182)$$\begin{array}{r}{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }=\{\begin{array}{ll}0,& (u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\le 0\mathrm{A}\mathrm{N}\mathrm{D}{u}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\ge 0)\mathrm{O}\mathrm{R}|u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }|<ϵ,\\ {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }>0,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }<0.\end{array}\end{array}$$

Predict $\tilde{w}$ to x-faces using a 1D extrapolation:

(183)$$\begin{array}{r}\begin{array}{rlr}{\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& =& {\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_x}^n+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}max(0,u_{\mathbf{i}-{\mathbf{e}}_x}^n)]{\mathrm{\Delta }}_x{\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_x}^n,\\ {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& =& {\tilde{w}}_{\mathbf{i}}-[\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}min(0,u_{\mathbf{i}}^n)]{\mathrm{\Delta }}_x{\tilde{w}}_{\mathbf{i}}^n.\end{array}\end{array}$$

Upwind based on $u^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$:

(184)$$\begin{array}{r}{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}=\{\begin{array}{ll}\frac{1}{2}({\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}+{\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}\right|<ϵ\\ {\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}>0,\\ {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}<0.\end{array}\end{array}$$

Predict $\tilde{w}$ to r-faces using a full-dimensional extrapolation:

(185)$$\begin{array}{r}\begin{array}{rlr}{\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }={\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& -& \frac{\mathrm{\Delta }t}{4h}(u_{\mathbf{i}-{\mathbf{e}}_r+\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+u_{\mathbf{i}-{\mathbf{e}}_r-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_r+\frac{1}{2}{\mathbf{e}}_x}-{\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_r-\frac{1}{2}{\mathbf{e}}_x})\\ & -& \frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+{\tilde{w}}_{\mathbf{i}-\frac{3}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})(w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{3}{2}{\mathbf{e}}_r})+\frac{\mathrm{\Delta }t}{2}{f}_{\tilde{w},\mathbf{i}-{\mathbf{e}}_r},\\ & & \\ {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }={\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& -& \frac{\mathrm{\Delta }t}{4h}(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{w}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}-{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x})\\ & -& \frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})(w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})+\frac{\mathrm{\Delta }t}{2}{f}_{\tilde{w},\mathbf{i}}.\end{array}\end{array}$$

Solve a Riemann problem:

(186)$$\begin{array}{r}{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }=\{\begin{array}{ll}0,& (w_L^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\le 0\mathrm{A}\mathrm{N}\mathrm{D}{w}_R^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\ge 0)\mathrm{O}\mathrm{R}|w_L^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+w_R^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }|<ϵ,\\ {\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& w_L^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+w_R^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }>0,\\ {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& w_L^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+w_R^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }<0.\end{array}\end{array}$$

3D Cartesian Case

This algorithm is more complicated than the 2D case since we include the effects of corner coupling.

1. Predict $\tilde{u}$ to y-faces using a 1D extrapolation.

2. Predict $\tilde{u}$ to r-faces using a 1D extrapolation.

3. Predict $\tilde{v}$ to x-faces using a 1D extrapolation.

4. Predict $\tilde{v}$ to r-faces using a 1D extrapolation.

5. Predict $\tilde{w}$ to x-faces using a 1D extrapolation.

6. Predict $\tilde{w}$ to y-faces using a 1D extrapolation.

7. Update prediction of $\tilde{u}$ to y-faces by accounting for $r$-derivatives.

8. Update prediction of $\tilde{u}$ to r-faces by accounting for $y$-derivatives.

9. Update prediction of $\tilde{v}$ to x-faces by accounting for $r$-derivatives.

10. Update prediction of $\tilde{v}$ to r-faces by accounting for $x$-derivatives.

11. Update prediction of $\tilde{w}$ to x-faces by accounting for $y$-derivatives.

12. Update prediction of $\tilde{w}$ to y-faces by accounting for $x$-derivatives.

13. Predict $\tilde{u}$ to x-faces using a full-dimensional extrapolation.

14. Predict $\tilde{v}$ to y-faces using a full-dimensional extrapolation.

15. Predict $\tilde{w}$ to r-faces using a full-dimensional extrapolation.

• Predict $\tilde{u}$ to y-faces using a 1D extrapolation.

  (187)$$\begin{array}{r}\begin{array}{rlr}{\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}& =& {\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_y}^n+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}max(0,v_{\mathbf{i}-{\mathbf{e}}_y}^n)]{\mathrm{\Delta }}_y{\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_y}^n,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}& =& {\tilde{u}}_{\mathbf{i}}-[\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}min(0,v_{\mathbf{i}}^n)]{\mathrm{\Delta }}_y{\tilde{u}}_{\mathbf{i}}^n.\end{array}\end{array}$$

  Upwind based on $v^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$:

  (188)$$\begin{array}{r}{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}=\{\begin{array}{ll}\frac{1}{2}({\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}+{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}),& \left|v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}\right|<ϵ\\ {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y},& v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}>0,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y},& v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}<0.\end{array}\end{array}$$

• Predict $\tilde{u}$ to r-faces using a 1D extrapolation.

• Predict $\tilde{v}$ to x-faces using a 1D extrapolation.

• Predict $\tilde{v}$ to r-faces using a 1D extrapolation.

• Predict $\tilde{w}$ to x-faces using a 1D extrapolation.

• Predict $\tilde{w}$ to y-faces using a 1D extrapolation.

• Update prediction of $\tilde{u}$ to y-faces by accounting for $r$-derivatives. The notation ${\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r}$ means state ${\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}$ that has been updated to account for transverse derives in the r-direction.

  (189)$$\begin{array}{r}\begin{array}{rlr}{\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r}& =& {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}-\frac{\mathrm{\Delta }t}{6h}(w_{\mathbf{i}-{\mathbf{e}}_y+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+w_{\mathbf{i}-{\mathbf{e}}_y-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_y+\frac{1}{2}{\mathbf{e}}_r}-{\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_y-\frac{1}{2}{\mathbf{e}}_r}),\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r}& =& {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}-\frac{\mathrm{\Delta }t}{6h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}).\end{array}\end{array}$$

  Upwind based on $v^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}$:

  (190)$$\begin{array}{r}{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r}=\{\begin{array}{ll}\frac{1}{2}({\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r}+{\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r}),& \left|v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}\right|<ϵ\\ {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r},& v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}>0,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r},& v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}<0.\end{array}\end{array}$$

• Update prediction of $\tilde{u}$ to r-faces by accounting for $y$-derivatives.

• Update prediction of $\tilde{v}$ to x-faces by accounting for $r$-derivatives.

• Update prediction of $\tilde{v}$ to r-faces by accounting for $x$-derivatives.

• Update prediction of $\tilde{w}$ to x-faces by accounting for $y$-derivatives.

• Update prediction of $\tilde{w}$ to y-faces by accounting for $x$-derivatives.

• Predict $\tilde{u}$ to x-faces using a full-dimensional extrapolation.

  (191)$$\begin{array}{r}\begin{array}{rlr}{\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }={\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{y|r}-{\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{y|r})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+w_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{r|y}-{\tilde{u}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{r|y})+\frac{\mathrm{\Delta }t}{2}{f}_{u,\mathbf{i}-{\mathbf{e}}_x},\\ & & \\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }={\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{y|r}-{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})({\tilde{u}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{r|y}-{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{r|y})+\frac{\mathrm{\Delta }t}{2}{f}_{u,\mathbf{i}}.\end{array}\end{array}$$

  Solve a Riemann problem:

  (192)$$\begin{array}{r}{\tilde{u}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }=\{\begin{array}{ll}0,& (u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\le 0\mathrm{A}\mathrm{N}\mathrm{D}{u}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\ge 0)\mathrm{O}\mathrm{R}|u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }|<ϵ,\\ {\tilde{u}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }>0,\\ {\tilde{u}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }<0.\end{array}\end{array}$$

• Predict $\tilde{v}$ to y-faces using a full-dimensional extrapolation.

• Predict $\tilde{w}$ to r-faces using a full-dimensional extrapolation. In this step, make sure you account for the $\partial w_0/\partial r$ term before solving the Riemann problem:

  (193)$$\begin{array}{r}\begin{array}{rlr}{\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }& =& {\tilde{w}}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})(w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})\\ {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }& =& {\tilde{w}}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}}+{\tilde{w}}_{\mathbf{i}-\frac{3}{2}{\mathbf{e}}_r}^{\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{N}\mathrm{S}})(w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{3}{2}{\mathbf{e}}_r})\end{array}\end{array}$$

3D Spherical Case

The spherical case is the same as the plane-parallel 3D Cartesian case, except the $\partial w_0/\partial r$ term enters in the full dimensional extrapolation for each direction. As in the plane-parallel case, make sure to upwind using the full velocity.

Computing ${\rho }^{{}^{\prime }\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},X_k^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},(\rho h{)}^{{}^{\prime }\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$, and ${\tilde{\mathbf{U}}}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ in MAESTROeX

We call make_edge_scal to compute ${\rho }^{{}^{\prime }\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},X_k^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},(\rho h{)}^{{}^{\prime }\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$, and ${\tilde{\mathbf{U}}}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ at each edge. The procedure is the same for each quantity, so we shall simply denote the scalar as $s$. We always need to compute ${\rho }^{\prime }$ and $X_k$ to faces, and the choice of energy prediction is as follows:

• For enthalpy_pred_type = 1, we predict $(\rho h{)}^{\prime }$ to faces.

• For enthalpy_pred_type = 2, we predict $h$ to faces.

• For enthalpy_pred_type = 3 and 4, we predict $T$ to faces.

• For enthalpy_pred_type = 5, we predict $h^{\prime }$ to faces.

• For enthalpy_pred_type = 6, we predict $T^{\prime }$ to faces.

We are using enthalpy_pred_type = 1 for now. The equations of motion are:

(194)$$\begin{array}{r}\begin{array}{rlr}\frac{\partial {\rho }^{\prime }}{\partial t}& =& -\mathbf{U}\cdot \mathrm{\nabla }{\rho }^{\prime }\underset{f_{{\rho }^{\prime }}}{\underbrace{-{\rho }^{\prime }\mathrm{\nabla }\cdot \mathbf{U}-\mathrm{\nabla }\cdot \left({\rho }_0\tilde{\mathbf{U}}\right)}},\\ \frac{\partial X_k}{\partial t}& =& -\mathbf{U}\cdot \mathrm{\nabla }{X}_k\text{(no forcing)},\\ \frac{\partial (\rho h{)}^{\prime }}{\partial t}& =& -\mathbf{U}\cdot \mathrm{\nabla }(\rho h{)}^{\prime }\underset{f_{(\rho h{)}^{\prime }}}{\underbrace{-(\rho h{)}^{\prime }\mathrm{\nabla }\cdot \mathbf{U}-\mathrm{\nabla }\cdot [(\rho h{)}_0\tilde{\mathbf{U}}]+(\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r)\frac{\partial p_0}{\partial r}+\mathrm{\nabla }\cdot k_{\mathrm{t}\mathrm{h}}\mathrm{\nabla }T}},\\ & & \\ \frac{\partial \tilde{\mathbf{U}}}{\partial t}& =& -\mathbf{U}\cdot \mathrm{\nabla }\tilde{\mathbf{U}}\underset{\text{forcing terms}}{\underbrace{-(\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r)\frac{\partial w_0}{\partial r}{\mathbf{e}}_r\underset{\text{terms included in }{\mathbf{f}}_{\tilde{\mathbf{U}}}}{\underbrace{-\frac{1}{\rho }\mathrm{\nabla }\pi +\frac{1}{{\rho }_0}\frac{\partial {\pi }_0}{\partial r}{\mathbf{e}}_r-\frac{(\rho -{\rho }_0)}{\rho }g{\mathbf{e}}_r}}}}.\end{array}\end{array}$$

2D Cartesian Case

1. Predict $s$ to r-faces using a 1D extrapolation.

2. Predict $s$ to x-faces using a full-dimensional extrapolation.

3. Predict $s$ to x-faces using a 1D extrapolation.

4. Predict $s$ to r-faces using a full-dimensional extrapolation.

• Predict $s$ to r-faces using a 1D extrapolation:

  (195)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& =& s_{\mathbf{i}-{\mathbf{e}}_r}^n+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}{w}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_r{s}_{\mathbf{i}-{\mathbf{e}}_r}^n,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}& =& s_{\mathbf{i}}-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}{w}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_r{s}_{\mathbf{i}}^n.\end{array}\end{array}$$

  Upwind based on $w^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

  (196)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}),& \left|w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

  Predict $s$ to x-faces using a full-dimensional extrapolation. First, the normal derivative and forcing terms:

  (197)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{\mathbf{i}-{\mathbf{e}}_x}^n+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_x{s}_{\mathbf{i}-{\mathbf{e}}_x}^n+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}-{\mathbf{e}}_x}^n\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{\mathbf{i}}^n-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_x{s}_{\mathbf{i}}^n+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}}^n\end{array}\end{array}$$

  Account for the transverse terms:

  if is_conservative then

  (198)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{2h}[{\left(w^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}-{\left(w^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}]-\frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}-{\mathbf{e}}_x}^n(u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}-u_{\mathbf{i}-\frac{3}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}})\\ & & \\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{2h}[{\left(w^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-{\left(w^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}]-\frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}}^n(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}-u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}})\end{array}\end{array}$$

  else

  (199)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r})\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})\end{array}\end{array}$$

  end if

• Account for the $\partial w_0/\partial r$ term:

  if is_vel and comp = 2 then

  (200)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+{\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})\end{array}\end{array}$$

  end if

• Upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$.

  (201)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

  Predict $s$ to x-faces using a 1D extrapolation:

  (202)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& =& s_{\mathbf{i}-{\mathbf{e}}_x}^n+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_x{s}_{\mathbf{i}-{\mathbf{e}}_x}^n,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& =& s_{\mathbf{i}}-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_x{s}_{\mathbf{i}}^n.\end{array}\end{array}$$

  Upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

  (203)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

  Predict $s$ to r-faces using a full-dimensional extrapolation. First, the normal derivative and forcing terms:

  (204)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{\mathbf{i}-{\mathbf{e}}_r}^n+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}{w}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_r{s}_{\mathbf{i}-{\mathbf{e}}_r}^n+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}-{\mathbf{e}}_r}^n\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{\mathbf{i}}^n-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}{w}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_r{s}_{\mathbf{i}}^n+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}}^n\end{array}\end{array}$$

  Account for the transverse terms: if is_conservative then

  (205)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{2h}[{\left(u^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}-{\mathbf{e}}_r+\frac{1}{2}{\mathbf{e}}_x}-{\left(u^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}-{\mathbf{e}}_r-\frac{1}{2}{\mathbf{e}}_x}]-\frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}-{\mathbf{e}}_r}^n(w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}-w_{\mathbf{i}-\frac{3}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})\\ & & \\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{2h}[{\left(u^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}-{\left(u^{\mathrm{M}\mathrm{A}\mathrm{C}}s\right)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}]-\frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}}^n(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}-w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})\end{array}\end{array}$$

  else

  (206)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}(u_{\mathbf{i}-{\mathbf{e}}_r+\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}+u_{\mathbf{i}-{\mathbf{e}}_r-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_r+\frac{1}{2}{\mathbf{e}}_x}-s_{\mathbf{i}-{\mathbf{e}}_r-\frac{1}{2}{\mathbf{e}}_x})\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}+u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x})\end{array}\end{array}$$

  end if

• Account for the $\partial w_0/\partial r$ term: if is_vel and comp = 2 then

  (207)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+{\tilde{w}}_{\mathbf{i}-\frac{3}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{3}{2}{\mathbf{e}}_r})\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})\end{array}\end{array}$$

  end if

• Upwind based on $w^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

  (208)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}),& \left|w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r},& w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

3D Cartesian Case

This algorithm is more complicated than the 2D case since we include the effects of corner coupling.

1. Predict $s$ to x-faces using a 1D extrapolation.

2. Predict $s$ to y-faces using a 1D extrapolation.

3. Predict $s$ to r-faces using a 1D extrapolation.

4. Update prediction of $s$ to x-faces by accounting for y-derivatives.

5. Update prediction of $s$ to x-faces by accounting for r-derivatives.

6. Update prediction of $s$ to y-faces by accounting for x-derivatives.

7. Update prediction of $s$ to y-faces by accounting for r-derivatives.

8. Update prediction of $s$ to r-faces by accounting for x-derivatives.

9. Update prediction of $s$ to r-faces by accounting for y-derivatives.

10. Predict $s$ to x-faces using a full-dimensional extrapolation.

11. Predict $s$ to y-faces using a full-dimensional extrapolation.

12. Predict $s$ to r-faces using a full-dimensional extrapolation.

• Predict $s$ to x-faces using a 1D extrapolation.

  (209)
  (210)

  Upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

  (211)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

• Predict $s$ to y-faces using a 1D extrapolation.

• Predict $s$ to r-faces using a 1D extrapolation.

• Update prediction of $s$ to x-faces by accounting for y-derivatives. The notation $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}$ means “state $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}$ that has been updated to account for the transverse derivatives in the $y$-direction”.

  if is_conservative then

  (212)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}-\frac{\mathrm{\Delta }t}{3h}[(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}],\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}-\frac{\mathrm{\Delta }t}{3h}[(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}].\end{array}\end{array}$$

  else

  (213)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}-\frac{\mathrm{\Delta }t}{6h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}),\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}-\frac{\mathrm{\Delta }t}{6h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}).\end{array}\end{array}$$

  end if

• Upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

  (214)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

• Update prediction of $s$ to x-faces by accounting for r-derivatives.

• Update prediction of $s$ to y-faces by accounting for x-derivatives.

• Update prediction of $s$ to y-faces by accounting for r-derivatives.

• Update prediction of $s$ to r-faces by accounting for x-derivatives.

• Update prediction of $s$ to r-faces by accounting for y-derivatives.

• Predict $s$ to x-faces using a full-dimensional extrapolation.

  if is_conservative then

  (215)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& -& \frac{\mathrm{\Delta }t}{2h}[(s^{y|r}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-(s^{y|r}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}]\\ & -& \frac{\mathrm{\Delta }t}{2h}[(s^{r|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}-(s^{r|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}]\\ & -& \frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}-{\mathbf{e}}_x}(u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}-u_{\mathbf{i}-\frac{3}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}})+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}-{\mathbf{e}}_x},\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& -& \frac{\mathrm{\Delta }t}{2h}[(s^{y|r}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-(s^{y|r}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}]\\ & -& \frac{\mathrm{\Delta }t}{2h}[(s^{r|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-(s^{r|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}]\\ & -& \frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}}(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}-u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}})+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}}.\end{array}\end{array}$$

  else

  (216)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{y|r}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{y|r})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{r|y}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{r|y})+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}-{\mathbf{e}}_x},\\ & & \\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{y|r}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|r})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{r|y}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{r|y})+\frac{\mathrm{\Delta }t}{2}{f}_{\mathbf{i}}.\end{array}\end{array}$$

  end if

• Account for the $\partial w_0/\partial r$ term:

  if is_vel and comp = 2 then

  (217)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+{\tilde{w}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}-\frac{\mathrm{\Delta }t}{4h}({\tilde{w}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}}+{\tilde{w}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{M}\mathrm{A}\mathrm{C}})(w_{0,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}-w_{0,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_r})\end{array}\end{array}$$

  end if

• Upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

  (218)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

• Predict $s$ to y-faces using a full-dimensional extrapolation.

• Predict $s$ to r-faces using a full-dimensional extrapolation.

3D Spherical Case

The spherical case is the same as the plane-parallel 3D Cartesian case, except the $\partial w_0/\partial r$ term enters in the full dimensional extrapolation for each direction when predicting velocity to faces. As in the plane-parallel case, make sure upwind based on the full velocity.

Computing ${\mathbf{U}}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }$ in VARDEN

2D Cartesian Case

• We do a 1D Taylor series extrapolation to get both components of velocity at the x-face:

  (219)$$u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=u_{\mathbf{i}-{\mathbf{e}}_x}+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}\mathrm{m}\mathrm{a}\mathrm{x}(0,u_{\mathbf{i}-{\mathbf{e}}_x})]{\mathrm{\Delta }}_x{u}_{\mathbf{i}-{\mathbf{e}}_x},$$
  (220)$$u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=u_{\mathbf{i}}+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}\mathrm{m}\mathrm{i}\mathrm{n}(0,u_{\mathbf{i}})]{\mathrm{\Delta }}_x{u}_{\mathbf{i}}.$$
  (221)$$v_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=v_{\mathbf{i}-{\mathbf{e}}_x}+[\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}\mathrm{m}\mathrm{a}\mathrm{x}(0,v_{\mathbf{i}-{\mathbf{e}}_x})]{\mathrm{\Delta }}_x{v}_{\mathbf{i}-{\mathbf{e}}_x},$$
  (222)

  We obtain the normal velocity using the Riemann problem:

  (223)$$\begin{array}{r}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=\{\begin{array}{ll}0,& (u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}\le 0\mathrm{A}\mathrm{N}\mathrm{D}{u}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}\ge 0)\mathrm{O}\mathrm{R}|u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}|<ϵ,\\ u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D},& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}>0,\\ u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D},& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}<0.\end{array}\end{array}$$

  We obtain the transverse velocity by upwinding based on $u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}$:

  (224)$$\begin{array}{r}{v}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=\{\begin{array}{ll}\frac{1}{2}(v_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}+v_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}\right|<ϵ\\ v_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}>0,\\ v_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}<0.\end{array}\end{array}$$

• We perform analogous operations to compute both components of velocity at the y-faces, ${\mathbf{U}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}$.

• Now we do a full-dimensional extrapolation to get the MAC velocity at the x-faces (note that we only compute the normal components):

  (225)$$\begin{array}{r}\begin{array}{rlr}{u}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }& =& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{1D}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{1D})(u_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{1D}-u_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{1D})+\frac{\mathrm{\Delta }t}{2}{f}_{u,\mathbf{i}-{\mathbf{e}}_x},\\ u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }& =& u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{1D}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D})(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{1D}-u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D})+\frac{\mathrm{\Delta }t}{2}{f}_{u,\mathbf{i}}.\end{array}\end{array}$$

  Then we solve a Riemann problem:

  (226)$$\begin{array}{r}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }=\{\begin{array}{ll}0,& (u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\le 0\mathrm{A}\mathrm{N}\mathrm{D}{u}_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }\ge 0)\mathrm{O}\mathrm{R}|u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }|<ϵ,\\ u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }>0,\\ u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast },& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }<0.\end{array}\end{array}$$

• We perform analogous operations to compute the normal velocity at the y-faces, $v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }$.

3D Cartesian Case

This is more complicated than the 2D case because we include corner coupling. We compute ${\mathbf{U}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D},{\mathbf{U}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}$, and ${\mathbf{U}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{1D}$ in an analogous manner as varden U_L^1D]-Eq.224. Then we compute an intermediate state, $u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z}$, which is described as “state $u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}$ that has been updated to account for the transverse derivatives in the z direction”, using:

(227)$$\begin{array}{r}\begin{array}{rlr}{u}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z}& =& u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}-\frac{\mathrm{\Delta }t}{6h}(w_{\mathbf{i}-{\mathbf{e}}_y+\frac{1}{2}{\mathbf{e}}_z}^{1D}+w_{\mathbf{i}-{\mathbf{e}}_y-\frac{1}{2}{\mathbf{e}}_z}^{1D})(u_{\mathbf{i}-{\mathbf{e}}_y+\frac{1}{2}{\mathbf{e}}_z}^{1D}-u_{\mathbf{i}-{\mathbf{e}}_y-\frac{1}{2}{\mathbf{e}}_z}^{1D}),\\ u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z}& =& u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}-\frac{\mathrm{\Delta }t}{6h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{1D}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{1D})(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{1D}-u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{1D}).\end{array}\end{array}$$

Then upwind based on $v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}$:

(228)$$\begin{array}{r}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z}=\{\begin{array}{ll}\frac{1}{2}(u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z}+u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z}),& \left|v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}\right|<ϵ\\ u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z},& v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}>0,\\ u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z},& v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}<0.\end{array}\end{array}$$

We use an analogous procedure to compute five more intermediate states, $u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{z|y},v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|z},v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{z|x},w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}$, and $w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|x}$. Then we do a full-dimensional extrapolation to get the MAC velocities at normal faces:

(229)$$\begin{array}{r}\begin{array}{rlr}{u}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }=u_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{1D}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{1D})(u_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{y|z}-u_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{y|z})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_z}^{1D}+w_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_z}^{1D})(u_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_z}^{z|y}-u_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_z}^{z|y})+\frac{\mathrm{\Delta }t}{2}{f}_{u,\mathbf{i}-{\mathbf{e}}_x},\\ u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }=u_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{1D}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D})(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{y|z}-u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{1D}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{1D})(u_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{z|y}-u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{z|y})+\frac{\mathrm{\Delta }t}{2}{f}_{u,\mathbf{i}}.\end{array}\end{array}$$

Then we use the Riemann solver given above for the 2D case (Eq.226) to compute $u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }$. We use an analogous procedure to obtain $v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }$ and $w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C},\ast }$.

Computing ${\mathbf{U}}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ and ${\rho }^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ in VARDEN

To compute ${\mathbf{U}}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$, VARDEN uses the exact same algorithm as the $s^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ case in MAESTROeX. The algorithm for ${\rho }^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ in VARDEN is slightly different than the $s^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ case in MAESTROeX since it uses a “conservative” formulation. Here, $s$ is used in place of either $\rho ,u,v$, or $w$ (in 3D).

2D Cartesian Case

The 1D extrapolation is:

(230)$$s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=s_{\mathbf{i}-{\mathbf{e}}_x}^n+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{2h}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_x{s}_{\mathbf{i}-{\mathbf{e}}_x}^n,$$

(231)$$s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=s_{\mathbf{i}}-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{2h}{u}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}){\mathrm{\Delta }}_x{s}_{\mathbf{i}}^n.$$

Then we upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

(232)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

We use an analogous procedure to obtain $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}$. Now we do a full-dimensional extrapolation of $s$ to each face. The extrapolation of a “non-conserved” $s$ to x-faces is:

(233)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{1D}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{1D})+\frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}-{\mathbf{e}}_x},\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{1D}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D})+\frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}}.\end{array}\end{array}$$

The extrapolation of a “conserved” $s$ to x-faces is:

(234)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{2h}[(s^{1D}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-(s^{1D}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}]\\ & -& \frac{\mathrm{\Delta }t}{2}{s}_{\mathbf{i}-{\mathbf{e}}_x}(\mathrm{\nabla }\cdot {\mathbf{U}}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x}+\frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}-{\mathbf{e}}_x}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}-v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})+\frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}-{\mathbf{e}}_x},\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{2h}[(s^{1D}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-(s^{1D}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}]\\ & -& \frac{\mathrm{\Delta }t}{2}{s}_{\mathbf{i}}(\mathrm{\nabla }\cdot {\mathbf{U}}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}}+\frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}-v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})+\frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}}.\end{array}\end{array}$$

Then we upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$.

(235)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

We use an analogous procedure to compute $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$.

3D Cartesian Case

This is more complicated than the 2D case because we include corner coupling. We first compute $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}$, $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}$, and $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{1D}$ in an analogous manner to Eq.230 and Eq.231. Then we compute six intermediate states, $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y},s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|z},s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|x},s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z},s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{z|x}$, and $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{z|y}$. For the “non-conservative case”, we use, for example:

(236)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{6h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{1D}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{1D}),\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{6h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{1D}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{1D}).\end{array}\end{array}$$

For the “conservative” case, we use, for example:

(237)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{3h}[(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}],\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}& =& s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}-\frac{\mathrm{\Delta }t}{3h}[(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-(s{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}].\end{array}\end{array}$$

Then we upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$:

(238)$$\begin{array}{r}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}=\{\begin{array}{ll}\frac{1}{2}(s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}+s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y}),& \left|u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}\right|<ϵ\\ s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0,\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{x|y},& u_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0.\end{array}\end{array}$$

We use an analogous procedure to compute the other five intermediate states. Now we do a full-dimensional extrapolation of $s$ to each face. The extrapolation of a “non-conserved” $s$ to x-faces is:

(239)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{y|z}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{y|z})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_z}^{z|y}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_z}^{z|y})\\ & +& \frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}-{\mathbf{e}}_x},\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{4h}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{y|z}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{y|z})\\ & -& \frac{\mathrm{\Delta }t}{4h}(w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{z|y}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{z|y})\\ & +& \frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}}.\end{array}\end{array}$$

The extrapolation of a “conserved” $s$ to x-faces is:

(240)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{L,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{2h}[(s^{y|z}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-(s^{y|z}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}]\\ & -& \frac{\mathrm{\Delta }t}{2h}[(s^{z|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_z}-(s^{z|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_z}]\\ & -& \frac{\mathrm{\Delta }t}{2}{s}_{\mathbf{i}-{\mathbf{e}}_x}(\mathrm{\nabla }\cdot {\mathbf{U}}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-{\mathbf{e}}_x}\\ & +& \frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}-{\mathbf{e}}_x}(v_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}-v_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}}-w_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}})\\ & +& \frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}-{\mathbf{e}}_x},\\ s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}=s_{R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^{1D}& -& \frac{\mathrm{\Delta }t}{2h}[(s^{y|z}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-(s^{y|z}{v}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}]\\ & -& \frac{\mathrm{\Delta }t}{2h}[(s^{z|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}-(s^{z|y}{w}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}]\\ & -& \frac{\mathrm{\Delta }t}{2}{s}_{\mathbf{i}}(\mathrm{\nabla }\cdot {\mathbf{U}}^{\mathrm{M}\mathrm{A}\mathrm{C}}{)}_{\mathbf{i}}\\ & +& \frac{\mathrm{\Delta }t}{2h}{s}_{\mathbf{i}}(v_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}-v_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{M}\mathrm{A}\mathrm{C}}+w_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}}-w_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}^{\mathrm{M}\mathrm{A}\mathrm{C}})\\ & +& \frac{\mathrm{\Delta }t}{2}{f}_{s,\mathbf{i}}.\end{array}\end{array}$$

Then we upwind based on $u^{\mathrm{M}\mathrm{A}\mathrm{C}}$, as in Eq.235. We use an analogous procedure to compute both $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^{\mathrm{E}\mathrm{D}\mathrm{G}\mathrm{E}}$ and $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_z}$.

ESTATE_FPU in GODUNOV_2D/3D.f

• First, the normal predictor.

  (241)$$\begin{array}{r}\begin{array}{rlr}{s}_L^x& =& s_{\mathbf{i}-{\mathbf{e}}_x}+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}-{\mathbf{e}}_x}+\underset{\text{IF USE\_MINION}}{\underbrace{\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}-{\mathbf{e}}_x}}}\\ s_R^x& =& s_{\mathbf{i}}-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}+\underset{\text{IF USE\_MINION}}{\underbrace{\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}}}}\end{array}\end{array}$$

  If USE_MINION and ICONSERVE then:

  (242)$$\begin{array}{r}\begin{array}{rlr}{s}_L^x& =& s_L^x-\frac{\mathrm{\Delta }t}{2}{s}_{\mathbf{i}-{\mathbf{e}}_x}{\text{DIVU}}_{\mathbf{i}-{\mathbf{e}}_x}\\ s_R^x& =& s_R^x-\frac{\mathrm{\Delta }t}{2}{s}_{\mathbf{i}}{\text{DIVU}}_{\mathbf{i}}\end{array}\end{array}$$

  Apply boundary conditions on $s_L^x$ and $s_R^x$. Then,

  (243)$$\begin{array}{r}{\text{s}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^x=\{\begin{array}{ll}{s}_L^x,& {\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}>0,\\ s_R^x,& \text{else}.\end{array}\end{array}$$

• Then, if $|{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}|\le ϵ$, we set $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^x=(s_L^x+s_R^x)/2$. The procedure to obtain $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^y$ is analogous.

• Now, the transverse terms.

  If ICONSERVE then:

  (244)$$\begin{array}{r}\begin{array}{rlr}{\text{sedge}}_L^x& =& s_{\mathbf{i}-{\mathbf{e}}_x}+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}-{\mathbf{e}}_x}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}-{\mathbf{e}}_x}\\ & & -\frac{\mathrm{\Delta }t}{2}[\frac{{\text{VEDGE}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}{s}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}^y-{\text{VEDGE}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}{s}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y}^y}{h_y}\\ & & -\frac{s_{\mathbf{i}-{\mathbf{e}}_x}({\text{VEDGE}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-{\text{VEDGE}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y})}{h_y}+s_{\mathbf{i}-{\mathbf{e}}_x}{\text{DIVU}}_{\mathbf{i}-{\mathbf{e}}_x}]\\ {\text{sedge}}_R^x& =& s_{\mathbf{i}}-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}}\\ & & -\frac{\mathrm{\Delta }t}{2}[\frac{{\text{VEDGE}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}{s}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}^y-{\text{VEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}{s}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}^y}{h_y}\\ & & -\frac{s_{\mathbf{i}}({\text{VEDGE}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-{\text{VEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y})}{h_y}+s_{\mathbf{i}}{\text{DIVU}}_{\mathbf{i}}]\end{array}\end{array}$$

• Now, define ${\text{VBAR}}_{\mathbf{i}}=({\text{VEDGE}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}+{\text{VEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y})/2$.

  If NOT ICONSERVE and ${\text{VEDGE}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}\cdot {\text{VEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}<0$ and ${\text{VBAR}}_{\mathbf{i}}<0$ then:

  (245)$$\begin{array}{r}\begin{array}{rl}{\text{sedge}}_L^x=s_{\mathbf{i}-{\mathbf{e}}_x}& +(\frac{1}{2}-\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}-{\mathbf{e}}_x}\\ & -\frac{\mathrm{\Delta }t}{2}\left[\frac{{\text{VBAR}}_{\mathbf{i}-{\mathbf{e}}_x}(s_{\mathbf{i}-{\mathbf{e}}_x+{\mathbf{e}}_y}-s_{\mathbf{i}-{\mathbf{e}}_x})}{h_y}\right]\end{array}\end{array}$$
  (246)$$\begin{array}{r}\begin{array}{rl}{\text{sedge}}_R^x=s_{\mathbf{i}}& -(\frac{1}{2}+\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}}\\ & -\frac{\mathrm{\Delta }t}{2}\left[\frac{{\text{VBAR}}_{\mathbf{i}}(s_{\mathbf{i}+{\mathbf{e}}_y}-s_{\mathbf{i}})}{h_y}\right]\end{array}\end{array}$$

  Else If NOT ICONSERVE and ${\text{VEDGE}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}\cdot {\text{VEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}<0$ and ${\text{VBAR}}_{\mathbf{i}}\ge 0$ then:

  (247)$$\begin{array}{r}\begin{array}{rlr}{\text{sedge}}_L^x=s_{\mathbf{i}-{\mathbf{e}}_x}& +& (\frac{1}{2}-\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}-{\mathbf{e}}_x}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}-{\mathbf{e}}_x}\\ & & -\frac{\mathrm{\Delta }t}{2}\left[\frac{{\text{VBAR}}_{\mathbf{i}-{\mathbf{e}}_x}(s_{\mathbf{i}-{\mathbf{e}}_x}-s_{\mathbf{i}-{\mathbf{e}}_x-{\mathbf{e}}_y})}{h_y}\right]\\ {\text{sedge}}_R^x=s_{\mathbf{i}}& -& (\frac{1}{2}+\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}}\\ & & -\frac{\mathrm{\Delta }t}{2}\left[\frac{{\text{VBAR}}_{\mathbf{i}}(s_{\mathbf{i}}-s_{\mathbf{i}-{\mathbf{e}}_y})}{h_y}\right]\end{array}\end{array}$$

  Else If NOT ICONSERVE and ${\text{VEDGE}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}\cdot {\text{VEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y}\ge 0$ then:

  (248)$$\begin{array}{r}\begin{array}{rl}{\text{sedge}}_L^x& =s_{\mathbf{i}-{\mathbf{e}}_x}+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}-{\mathbf{e}}_x}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}-{\mathbf{e}}_x}\\ & & -\frac{\mathrm{\Delta }t}{2}\left[\frac{({\text{VEDGE}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}+{\text{VEDGE}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y})}{2h_y}\right]\\ {\text{sedge}}_R^x& =& s_{\mathbf{i}}-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{h_x}{\text{UEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}}\\ & & -\frac{\mathrm{\Delta }t}{2}\left[\frac{({\text{VEDGE}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}+{\text{VEDGE}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y})}{2h_y}\right]\end{array}\end{array}$$

• Finally, upwind analogous to Eq.243 to get ${\text{sedge}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}$.

ESTATE in GODUNOV_2D/3D.f

First, the normal predictor.

(249)$$\begin{array}{r}\begin{array}{rlr}{s}_L^x& =& s_{\mathbf{i}-{\mathbf{e}}_x}+(\frac{1}{2}-\frac{\mathrm{\Delta }t}{h_x}{u}_{\mathbf{i}-{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}-{\mathbf{e}}_x}\\ s_R^x& =& s_{\mathbf{i}}-(\frac{1}{2}+\frac{\mathrm{\Delta }t}{h_x}{u}_{\mathbf{i}}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}\end{array}\end{array}$$

If USE_MINION then:

(250)$$\begin{array}{r}\begin{array}{rlr}{s}_L^x& =& s_L^x+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}-{\mathbf{e}}_x}\\ s_R^x& =& s_R^x+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}}\end{array}\end{array}$$

Apply boundary conditions on $s_L^x$ and $s_R^x$. Then,

(251)$$\begin{array}{r}{\text{s}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^x=\{\begin{array}{ll}{s}_L^x,& {\text{UAD}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}>0,\\ s_R^x,& \text{else}.\end{array}\end{array}$$

Then, if $|{\text{UAD}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}|\le ϵ$, we set $s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}^x=(s_L^x+s_R^x)/2$.

(252)$$\begin{array}{r}\begin{array}{rlr}{\text{sedge}}_L^x=s_{\mathbf{i}-{\mathbf{e}}_x}& +& (\frac{1}{2}-\frac{\mathrm{\Delta }t}{h_x}{u}_{\mathbf{i}-{\mathbf{e}}_x}){\mathrm{\Delta }}^x{s}_{\mathbf{i}-{\mathbf{e}}_x}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}-{\mathbf{e}}_x}\\ & & -\frac{\mathrm{\Delta }t}{2}\left[\frac{({\text{VAD}}_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}+{\text{VAD}}_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y})(s_{\mathbf{i}-{\mathbf{e}}_x+\frac{1}{2}{\mathbf{e}}_y}-s_{\mathbf{i}-{\mathbf{e}}_x-\frac{1}{2}{\mathbf{e}}_y})}{2h_y}\right]\\ {\text{sedge}}_R^x=s_{\mathbf{i}}& -& (\frac{1}{2}+\frac{\mathrm{\Delta }t}{h_x}{u}_{\mathbf{i}}){\mathrm{\Delta }}^x{s}_{\mathbf{i}}+\frac{\mathrm{\Delta }t}{2}{\text{TFORCES}}_{\mathbf{i}}\\ & & -\frac{\mathrm{\Delta }t}{2}\left[\frac{({\text{VAD}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}+{\text{VAD}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y})(s_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_y}-s_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_y})}{2h_y}\right]\end{array}\end{array}$$

Note that the 2D and 3D algorithms are different - in 3D the transverse terms use upwinding analogous to Eq.245-Eq.248, using UAD instead of UEDGE. Finally, upwind analogous to Eq.251 to get ${\text{sedge}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}$, but use UEDGE instead of UAD.

Piecewise Parabolic Method (PPM)

Consider a scalar, $s$, which we wish to predict to time-centered edges. The PPM method is an improvement over the piecewise-linear method. Using our notation, we modify equations (Eq.209 and Eq.210) in Section 4 to obtain better estimates for the time-centered 1D edge states, $s_{L/R,\mathbf{i}-\frac{1}{2}{\mathbf{e}}_x}$, etc.. Once these states are obtained, we continue with the full-dimensional extrapolations as described before.

The PPM method is described in a series of papers:

• Colella and Woodward 1984 - describes the basic method.

• Miller and Colella 2002 - describes how to apply PPM to a multidimensional unsplit Godunov method and generalizes the characteristic tracing for more complicated systems. Note that we only need to upwind based on the fluid velocity, so we don’t need to use fancy characteristic tracing.

• Colella and Sekora 2008 - describes new fancy quadratic limiters. There are several errors in the printed text, so we have implemented a corrected version from Phil Colella.

Here are the steps for the $x$-direction. For simplicity, we replace the vector index notation with a simple scalar notation ($\mathbf{i}+{\mathbf{e}}_x\to i+1$, etc.).

• Step 1: Compute $s_{i,+}$ and $s_{i,-}$, which are spatial interpolations of $s$ to the hi and lo faces of cell $i$, respectively. See Sections 9.1 and 9.2 for the two options.

• Step 2: Construct a quadratic profile within each cell.

  (253)$$s_i^I(x)=s_{i,-}+\xi [s_{i,+}-s_{i,-}+s_{6,i}(1-\xi )],$$
  (254)$$s_{6,i}=6s_i-3(s_{i,-}+s_{i,+}),$$
  (255)$$\xi =\frac{x-(i-\frac{1}{2})h}{h},0\le \xi \le 1.$$

• Step 3: Integrate quadratic profiles to get the average value swept over the face over time. Define the following integrals, where $\sigma =|u|\mathrm{\Delta }t/h$:

  (256)$$\begin{array}{r}\begin{array}{rlr}{\mathcal{I}}_{i,+}(\sigma )& =& \frac{1}{\sigma h}{\int }_{(i+\frac{1}{2})h-\sigma h}^{(i+\frac{1}{2})h}{s}_i^I(x)dx\\ {\mathcal{I}}_{i,-}(\sigma )& =& \frac{1}{\sigma h}{\int }_{(i-\frac{1}{2})h}^{(i-\frac{1}{2})h+\sigma h}{s}_i^I(x)dx\end{array}\end{array}$$

  Plugging in ([Quadratic Interp]) gives:

  (257)$$\begin{array}{r}\begin{array}{rlr}{\mathcal{I}}_{i,+}(\sigma )& =& s_{j,+}-\frac{\sigma }{2}[s_{j,+}-s_{j,-}-(1-\frac{2}{3}\sigma )s_{6,i}],\\ {\mathcal{I}}_{i,-}(\sigma )& =& s_{j,-}+\frac{\sigma }{2}[s_{j,+}-s_{j,-}+(1-\frac{2}{3}\sigma )s_{6,i}].\end{array}\end{array}$$

• Step 4: Obtain 1D edge states. Perform a 1D extrapolation, without source terms, to get left and right edge states. Add the source terms later if desired/necessary.

  (258)$$\begin{array}{r}\begin{array}{rlr}{s}_{L,i-\frac{1}{2}}& =& \{\begin{array}{ll}{\mathcal{I}}_{i-1,+}(\sigma ),& u_{i-1}\text{or}{u}_{i-\frac{1}{2}}^{\mathrm{M}\mathrm{A}\mathrm{C}}>0\\ s_{i-1},& \text{else}.\end{array}\\ s_{R,i-\frac{1}{2}}& =& \{\begin{array}{ll}{\mathcal{I}}_{i,-}(\sigma ),& u_i\text{or}{u}_{i-\frac{1}{2}}^{\mathrm{M}\mathrm{A}\mathrm{C}}<0\\ s_i,& \text{else}.\end{array}\end{array}\end{array}$$

Colella and Woodward Based Approach

Spatially interpolate $s$ to edges. Use a 4th-order interpolation in space with van Leer limiting to obtain edge values:

(259)$$s_{i+\frac{1}{2}}=\frac{1}{2}(s_i+s_{i+1})-\frac{1}{6}(\delta s_{i+1}^{vL}-\delta s_i^{vL}),$$

(260)$$\delta s_i=\frac{1}{2}(s_{i+1}-s_{i-1}),$$

(261)$$\begin{array}{r}\delta s_i^{vL}=\{\begin{array}{ll}\text{sign}(\delta s_i)min(|\delta s_i|,2|s_{i+1}-s_i|,2|s_i-s_{i-1}|),& \mathrm{i}\mathrm{f}(s_{i+1}-s_i)(s_i-s_{i-1})>0,\\ 0,& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\end{array}\end{array}$$

A more compact way of writing this is

(262)$$s=\text{sign}(\delta s_i),$$

(263)$$\delta s_i^{vL}=smax\{min[s\delta s_i,2s(s_{i+1}-s_i),2s(s_i-s_{i-1})],0\}$$

Without the limiters, Eq.259 is the familiar 4th-order spatial interpolation formula:

(264)$$s_{i+\frac{1}{2}}=\frac{7}{12}(s_{i+1}+s_i)-\frac{1}{12}(s_{i+2}+s_{i-1}).$$

Next, we must ensure that $s_{i+\frac{1}{2}}$ lies between the adjacent cell-centered values:

(265)$$min(s_i,s_{i+1})\le s_{i+\frac{1}{2}}\le max(s_i,s_{i+1}).$$

In anticipation of further limiting, we set double-valued face-centered values:

(266)$$s_{i,+}=s_{i+1,-}=s_{i+\frac{1}{2}}.$$

Modify $s_{i,\pm }$ using a quadratic limiter. First, we test whether $s_i$ is a local extreumum with the condition:

(267)$$(s_{i,+}-s_i)(s_i-s_{i,-})\le 0,$$

If this condition is true, we constrain $s_{i,\pm }$ by setting $s_{i,+}=s_{i,-}=s_i$. If not, we then apply a second test to determine whether $s_i$ is sufficiently close to $s_{i,\pm }$ so that a quadratic interpolate would contain a local extremum. We define ${\alpha }_{i,\pm }=s_{i,\pm }-s_i$. If one of $|{\alpha }_{i,\pm }|\ge 2|{\alpha }_{i,\mp }|$ holds, then for that choice of $\pm =+,-$ we set:

(268)$$s_{i,\pm }=3s_i-2s_{i,\mp }.$$

Colella and Sekora Based Approach

• Spatially interpolate $s$ to edges.

  Use a 4th-order interpolation in space to obtain edge values:

  (269)$$s_{i+\frac{1}{2}}=\frac{7}{12}(s_{i+1}+s_i)-\frac{1}{12}(s_{i+2}+s_{i-1}).$$

  Then, if $(s_{i+\frac{1}{2}}-s_i)(s_{i+1}-s_{i+\frac{1}{2}})<0$, we limit $s_{i+\frac{1}{2}}$ using a nonlinear combination of approximations to the second derivative. First, define:

  (270)$$\begin{array}{r}\begin{array}{rlr}(D^{2}s{)}_{i+\frac{1}{2}}& =& 3(s_i-2s_{i+\frac{1}{2}}+s_{i+1})\\ (D^{2}s{)}_{i+\frac{1}{2},L}& =& s_{i-1}-2s_i+s_{i+1}\\ (D^{2}s{)}_{i+\frac{1}{2},R}& =& s_i-2s_{i+1}+s_{i+2}\end{array}\end{array}$$

  Then, define

  (271)$$s=\text{sign}[(D^{2}s{)}_{i+\frac{1}{2}}],$$
  (272)$$(D^{2}s{)}_{i+\frac{1}{2},\text{lim}}=smax\{min[Cs(D^{2}s{)}_{i+\frac{1}{2},L},Cs(D^{2}s{)}_{i+\frac{1}{2},R},s(D^{2}s{)}_{i+\frac{1}{2}}],0\},$$

  where $C=1.25$ was used in Colella and Sekora. Then,

  (273)$$s_{i+\frac{1}{2}}=\frac{1}{2}(s_i+s_{i+1})-\frac{1}{6}(D^{2}s{)}_{i+\frac{1}{2},\text{lim}}.$$

  Now we implement Phil’s new version of the algorithm to eliminate sensitivity to roundoff. First we need to detect whether a particular cell corresponds to an “extremum”. There are two tests. For the first test, define

  (274)$${\alpha }_{i,\pm }=s_{i\pm \frac{1}{2}}-s_i.$$

  If ${\alpha }_{i,+}{\alpha }_{i,-}\ge 0$, then we are at an extremum. We apply the second test if either $|{\alpha }_{i,\pm }|>2|{\alpha }_{i,\mp }|$. Then, we define:

  (275)$$\begin{array}{r}\begin{array}{rlr}(Ds{)}_{i,\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e},-}& =& s_{i-\frac{1}{2}}-s_{i-\text{sfrac}32}\\ (Ds{)}_{i,\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e},+}& =& s_{i+\text{sfrac}32}-s_{i+\frac{1}{2}}\end{array}\end{array}$$
  (276)$$(Ds{)}_{i,\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e},\mathrm{m}\mathrm{i}\mathrm{n}}=min[|(Ds{)}_{i,\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e},-}|,|(Ds{)}_{i,\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e},+}|].$$
  (277)$$\begin{array}{r}\begin{array}{rlr}(Ds{)}_{i,\mathrm{c}\mathrm{c},-}& =& s_i-s_{i-1}\\ (Ds{)}_{i,\mathrm{c}\mathrm{c},+}& =& s_{i+1}-s_i\end{array}\end{array}$$
  (278)$$(Ds{)}_{i,\mathrm{c}\mathrm{c},\mathrm{m}\mathrm{i}\mathrm{n}}=min[|(Ds{)}_{i,\mathrm{c}\mathrm{c},-}|,|(Ds{)}_{i,\mathrm{c}\mathrm{c},+}|].$$

  If $(Ds{)}_{i,\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e},\mathrm{m}\mathrm{i}\mathrm{n}}\ge (Ds{)}_{i,\mathrm{c}\mathrm{c},\mathrm{m}\mathrm{i}\mathrm{n}}$, set $(Ds{)}_{i,\pm }=(Ds{)}_{i,\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e},\pm }$. Otherwise, set $(Ds{)}_{i,\pm }=(Ds{)}_{i,\mathrm{c}\mathrm{c},\pm }$. Finally, we are at an extreumum if $(Ds{)}_{i,+}(Ds{)}_{i,-}\le 0$.

• Now that we have finished the extremum tests, if we are at an extremum, we scale ${\alpha }_{i,\pm }$. First, we define

  (279)$$\begin{array}{r}\begin{array}{rlr}(D^{2}s{)}_i& =& 6({\alpha }_{i,+}+{\alpha }_{i,-})\\ (D^{2}s{)}_{i,L}& =& s_{i-2}-2s_{i-1}+s_i\\ (D^{2}s{)}_{i,R}& =& s_i-2s_{i+1}+s_{i+2}\\ (D^{2}s{)}_{i,C}& =& s_{i-1}-2s_i+s_{i+1}\end{array}\end{array}$$

  Then, define

  (280)$$s=\text{sign}[(D^{2}s{)}_i],$$
  (281)$$(D^{2}s{)}_{i,\text{lim}}=max\{min[s(D^{2}s{)}_i,Cs(D^{2}s{)}_{i,L},Cs(D^{2}s{)}_{i,R},Cs(D^{2}s{)}_{i,C}],0\}.$$

  Then,

  (282)$${\alpha }_{i,\pm }=\frac{{\alpha }_{i,\pm }(D^{2}s{)}_{i,\text{lim}}}{max[|(D^{2}s{)}_i|,1\times {10}^{-10}]}$$

  Otherwise, if we are not at an extremum and $|{\alpha }_{i,\pm }|>2|{\alpha }_{i,\mp }|$, then define

  (283)$$s=\text{sign}({\alpha }_{i,\mp })$$
  (284)$$\delta {\mathcal{I}}_{\text{ext}}=\frac{-{\alpha }_{i,\pm }^2}{4({\alpha }_{j,+}+{\alpha }_{j,-})},$$
  (285)$$\delta s=s_{i\mp 1}-s_i,$$

  If $s\delta {\mathcal{I}}_{\text{ext}}\ge s\delta s$, then we perform the following test. If $s\delta s-{\alpha }_{i,\mp }\ge 1\times {10}^{-10}$, then

  (286)$${\alpha }_{i,\pm }=-2\delta s-2s{[(\delta s{)}^2-\delta s{\alpha }_{i,\mp }]}^{\frac{1}{2}}$$

  otherwise,

  (287)$${\alpha }_{i,\pm }=-2{\alpha }_{i,\mp }$$

  Finally, $s_{i,\pm }=s_i+{\alpha }_{i,\pm }$.