Godunov Interface States
These are working notes for the Godunov step in MAESTROeX and VARDEN.
MAESTROeX Notation
For 2D, \(\Ub = (u,w)\) and \(\Ubt = (\ut,\wt)\). Note that \(u = \ut\). We will use the shorthand \(\ib = (x,r)\).
For 3D plane parallel, \(\Ub = (u,v,w)\) and \(\Ubt = (\ut,\vt,\wt)\). Note that \(u = \ut\) and \(v = \vt\). We will use the shorthand \(\ib = (x,y,r)\).
For 3D spherical, \(\Ub = (u,v,w)\) and \(\Ubt = (\ut,\vt,\wt)\). We will use the shorthand \(\ib = (x,y,z)\).
Computing \(\Ub\) From \(\Ubt\)
For plane-parallel problems, in order to compute \(w\) from \(\wt\), we use simple averaging
For spherial problems, in order to compute \(\Ub\) from \(\Ubt\), we first map \(w_0\) onto \(w_0^{\mac}\) using put_w0_on_edges, where \(w_0^{\mac}\) only contains normal velocities at each face. Then we construct \(\Ub\) by using
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:
For spherical problems, we compute the radial bin centered gradient using
Then we put \(\partial w_0/\partial r\) onto a Cartesian grid using put_1d_array_on_cart_3d_sphr.
Computing \(\Ubt^{\trans}\) in MAESTROeX
(169)\[\frac{\partial\Ubt}{\partial t} = -\Ub\cdot\nabla\Ubt \underbrace{- (\Ubt\cdot\eb_r)\frac{\partial w_0}{\partial r}\eb_r - \frac{1}{\rho}\nabla\pi + \frac{1}{\rho_0}\frac{\partial\pi_0}{\partial r}\eb_r - \frac{(\rho-\rho_0)}{\rho}g\eb_r}_{\hbox{forcing terms}}.\label{Perturbational Velocity Equation}\]We extrapolate each velocity component to edge-centered, time-centered locations. For example,
(170)\[\begin{split}\begin{aligned} \ut_{R,\ib-\half\eb_x} &=& \ut_{\ib}^n + \frac{h}{2}\frac{\partial\ut_{\ib}^n}{\partial x} + \frac{\dt}{2}\frac{\partial\ut_{\ib}^n}{\partial t} \nonumber \\ &=& \ut_{\ib}^n + \frac{h}{2}\frac{\partial\ut_{\ib}^n}{\partial x} + \frac{\dt}{2} \left(-\ut_{\ib}^n\frac{\partial\ut_{\ib}^n}{\partial x} - \wt_{\ib}^n\frac{\partial\ut_{\ib}^n}{\partial r} + \text{forcing terms}\right)\end{aligned}\end{split}\]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)\[\ut_{R,\ib-\half\eb_x} = \ut_{\ib}^n + \left[\frac{1}{2} - \frac{\dt}{2h}\min(0,\ut_{\ib}^n)\right]\partial\ut_{\ib}^n\]
2D Cartesian Case
We predict \(\ut\) to x-faces using a 1D extrapolation:
We pick the final trans states using a Riemann solver:
We predict \(\wt\) to r-faces using a 1D extrapolation:
We pick the final \(\trans\) states using a Riemann solver, noting that we upwind based on the full velocity.
3D Cartesian Case
We use the exact same procedure in 2D and 3D to compute \(\ut^{\trans}\) and \(\wt^{\trans}\). The procedure for computing \(\vt^{\trans}\) is analogous to computing \(\ut^{\trans}\). We predict \(\vt\) to y-faces using the 1D extrapolation:
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 \(\Ubt^{\mac,*}\) in MAESTROeX
(178)\[\frac{\partial\Ubt}{\partial t} = -\Ub\cdot\nabla\Ubt \underbrace{- (\Ubt\cdot\eb_r)\frac{\partial w_0}{\partial r}\eb_r \underbrace{- \frac{1}{\rho}\nabla\pi + \frac{1}{\rho_0}\frac{\partial\pi_0}{\partial r}\eb_r - \frac{(\rho-\rho_0)}{\rho}g\eb_r}_{\hbox{terms included in $\fb_{\Ubt}$}}}_{\hbox{forcing terms}}.\]Note that the \(\partial w_0/\partial r\) term is treated like a forcing term, but it is not actually part of \(\fb_{\Ubt}\). We make use of the 1D extrapolations used to compute \(\Ubt^{\trans}\) (\(\ut_{L/R,\ib-\half\eb_x}\), \(\vt_{L/R,\ib-\half\eb_y}\), and \(\wt_{L/R,\ib-\half\eb_r}\)), as well as the “\(\trans\)” states (\(\ut_{\ib-\half\eb_x}^{\trans}\), \(\vt_{\ib-\half\eb_y}^{\trans}\), and \(\wt_{\ib-\half\eb_r}^{\trans}\))
2D Cartesian Case
Predict \(\ut\) to r-faces using a 1D extrapolation.
Predict \(\ut\) to x-faces using a full-dimensional extrapolation.
Predict \(\wt\) to x-faces using a 1D extrapolation.
Predict \(\wt\) to r-faces using a full-dimensional extrapolation.
Predict \(\ut\) to r-faces using a 1D extrapolation:
Upwind based on \(w^{\trans}\):
Predict \(\ut\) to x-faces using a full-dimensional extrapolation,
Solve a Riemann problem:
Predict \(\wt\) to x-faces using a 1D extrapolation:
Upwind based on \(u^{\trans}\):
Predict \(\wt\) to r-faces using a full-dimensional extrapolation:
Solve a Riemann problem:
3D Cartesian Case
This algorithm is more complicated than the 2D case since we include the effects of corner coupling.
Predict \(\ut\) to y-faces using a 1D extrapolation.
Predict \(\ut\) to r-faces using a 1D extrapolation.
Predict \(\vt\) to x-faces using a 1D extrapolation.
Predict \(\vt\) to r-faces using a 1D extrapolation.
Predict \(\wt\) to x-faces using a 1D extrapolation.
Predict \(\wt\) to y-faces using a 1D extrapolation.
Update prediction of \(\ut\) to y-faces by accounting for \(r\)-derivatives.
Update prediction of \(\ut\) to r-faces by accounting for \(y\)-derivatives.
Update prediction of \(\vt\) to x-faces by accounting for \(r\)-derivatives.
Update prediction of \(\vt\) to r-faces by accounting for \(x\)-derivatives.
Update prediction of \(\wt\) to x-faces by accounting for \(y\)-derivatives.
Update prediction of \(\wt\) to y-faces by accounting for \(x\)-derivatives.
Predict \(\ut\) to x-faces using a full-dimensional extrapolation.
Predict \(\vt\) to y-faces using a full-dimensional extrapolation.
Predict \(\wt\) to r-faces using a full-dimensional extrapolation.
Predict \(\ut\) to y-faces using a 1D extrapolation.
(187)\[\begin{split}\begin{aligned} \ut_{L,\ib-\half\eb_y} &=& \ut_{\ib-\eb_y}^n + \left[\half - \frac{\dt}{2h}\max(0,v_{\ib-\eb_y}^n)\right]\Delta_y \ut_{\ib-\eb_y}^n, \\ \ut_{R,\ib-\half\eb_y} &=& \ut_{\ib} - \left[\half + \frac{\dt}{2h}\min(0,v_{\ib}^n)\right]\Delta_y \ut_{\ib}^n.\end{aligned}\end{split}\]Upwind based on \(v^{\trans}\):
(188)\[\begin{split}\ut_{\ib-\half\eb_y} = \begin{cases} \half\left(\ut_{L,\ib-\half\eb_y} + \ut_{R,\ib-\half\eb_y}\right), & \left|v^{\trans}_{\ib-\half\eb_y}\right| < \epsilon \\ \ut_{L,\ib-\half\eb_y}, & v^{\trans}_{\ib-\half\eb_y} > 0, \\ \ut_{R,\ib-\half\eb_y}, & v^{\trans}_{\ib-\half\eb_y} < 0. \\ \end{cases}\end{split}\]Predict \(\ut\) to r-faces using a 1D extrapolation.
Predict \(\vt\) to x-faces using a 1D extrapolation.
Predict \(\vt\) to r-faces using a 1D extrapolation.
Predict \(\wt\) to x-faces using a 1D extrapolation.
Predict \(\wt\) to y-faces using a 1D extrapolation.
Update prediction of \(\ut\) to y-faces by accounting for \(r\)-derivatives. The notation \(\ut_{\ib-\half\eb_y}^{y|r}\) means state \(\ut_{\ib-\half\eb_y}\) that has been updated to account for transverse derives in the r-direction.
(189)\[\begin{split}\begin{aligned} \ut_{L,\ib-\half\eb_y}^{y|r} &=& \ut_{L,\ib-\half\eb_y} - \frac{\dt}{6h}\left(w_{\ib-\eb_y+\half\eb_r}^{\trans}+w_{\ib-\eb_y-\half\eb_r}^{\trans}\right)\left(\ut_{\ib-\eb_y+\half\eb_r}-\ut_{\ib-\eb_y-\half\eb_r}\right), \\ \ut_{R,\ib-\half\eb_y}^{y|r} &=& \ut_{R,\ib-\half\eb_y} - \frac{\dt}{6h}\left(w_{\ib+\half\eb_r}^{\trans}+w_{\ib-\half\eb_r}^{\trans}\right)\left(\ut_{\ib+\half\eb_r}-\ut_{\ib-\half\eb_r}\right).\end{aligned}\end{split}\]Upwind based on \(v^{\trans}\):
(190)\[\begin{split}\ut_{\ib-\half\eb_y}^{y|r} = \begin{cases} \half\left(\ut_{L,\ib-\half\eb_y}^{y|r} + \ut_{R,\ib-\half\eb_y}^{y|r}\right), & \left|v_{\ib-\half\eb_y}^{\trans}\right| < \epsilon \\ \ut_{L,\ib-\half\eb_y}^{y|r}, & v_{\ib-\half\eb_y}^{\trans} > 0, \\ \ut_{R,\ib-\half\eb_y}^{y|r}, & v_{\ib-\half\eb_y}^{\trans} < 0. \end{cases}\end{split}\]Update prediction of \(\ut\) to r-faces by accounting for \(y\)-derivatives.
Update prediction of \(\vt\) to x-faces by accounting for \(r\)-derivatives.
Update prediction of \(\vt\) to r-faces by accounting for \(x\)-derivatives.
Update prediction of \(\wt\) to x-faces by accounting for \(y\)-derivatives.
Update prediction of \(\wt\) to y-faces by accounting for \(x\)-derivatives.
Predict \(\ut\) to x-faces using a full-dimensional extrapolation.
(191)\[\begin{split}\begin{aligned} \ut_{L,\ib-\half\eb_x}^{\mac,*} = \ut_{L,\ib-\half\eb_x} &-& \frac{\dt}{4h}\left(v_{\ib-\eb_x+\half\eb_y}^{\trans}+v_{\ib-\eb_x-\half\eb_y}^{\trans}\right)\left(\ut_{\ib-\eb_x+\half\eb_y}^{y|r}-\ut_{\ib-\eb_x-\half\eb_y}^{y|r}\right) \nonumber \\ &-& \frac{\dt}{4h}\left(w_{\ib-\eb_x+\half\eb_r}^{\trans}+w_{\ib-\eb_x-\half\eb_r}^{\trans}\right)\left(\ut_{\ib-\eb_x+\half\eb_r}^{r|y}-\ut_{\ib-\eb_x-\half\eb_r}^{r|y}\right) + \frac{\dt}{2}f_{u,\ib-\eb_x}, \nonumber \\ && \\ \ut_{R,\ib-\half\eb_x}^{\mac,*} = \ut_{R,\ib-\half\eb_x} &-& \frac{\dt}{4h}\left(v_{\ib+\half\eb_y}^{\trans}+v_{\ib-\half\eb_y}^{\trans}\right)\left(\ut_{\ib+\half\eb_y}^{y|r}-\ut_{\ib-\half\eb_y}^{y|r}\right) \nonumber \\ &-& \frac{\dt}{4h}\left(w_{\ib+\half\eb_r}^{\trans}+w_{\ib-\half\eb_r}^{\trans}\right)\left(\ut_{\ib+\half\eb_r}^{r|y}-\ut_{\ib-\half\eb_r}^{r|y}\right) + \frac{\dt}{2}f_{u,\ib}.\end{aligned}\end{split}\]Solve a Riemann problem:
(192)\[\begin{split}\ut_{\ib-\half\eb_x}^{\mac,*} = \begin{cases} 0, & \left(u_{L,\ib-\half\eb_x}^{\mac,*} \le 0 ~~ {\rm AND} ~~ u_{R,\ib-\half\eb_x}^{\mac,*} \ge 0\right) ~~ {\rm OR} ~~ \left|u_{L,\ib-\half\eb_x}^{\mac,*} + u_{R,\ib-\half\eb_x}^{\mac,*}\right| < \epsilon, \\ \ut_{L,\ib-\half\eb_x}^{\mac,*}, & u_{L,\ib-\half\eb_x}^{\mac,*} + u_{R,\ib-\half\eb_x}^{\mac,*} > 0, \\ \ut_{R,\ib-\half\eb_x}^{\mac,*}, & u_{L,\ib-\half\eb_x}^{\mac,*} + u_{R,\ib-\half\eb_x}^{\mac,*} < 0. \end{cases}\end{split}\]Predict \(\vt\) to y-faces using a full-dimensional extrapolation.
Predict \(\wt\) 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{split}\begin{aligned} \wt_{L,\ib-\half\eb_r}^{\mac,*} &=& \wt_{L,\ib-\half\eb_r}^{\mac,*} - \frac{\dt}{4h}\left(\wt^{\trans}_{\ib+\half\eb_r} + \wt^{\trans}_{\ib-\half\eb_r}\right)\left(w_{0,\ib+\half\eb_r}-w_{0,\ib-\half\eb_r}\right) \\ \wt_{R,\ib-\half\eb_r}^{\mac,*} &=& \wt_{R,\ib-\half\eb_r}^{\mac,*} - \frac{\dt}{4h}\left(\wt^{\trans}_{\ib-\half\eb_r} + \wt^{\trans}_{\ib-\frac{3}{2}\eb_r}\right)\left(w_{0,\ib-\half\eb_r}-w_{0,\ib-\frac{3}{2}\eb_r}\right)\end{aligned}\end{split}\]
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^{'\edge}, X_k^{\edge},(\rho h)^{'\edge}\), and \(\Ubt^{\edge}\) in MAESTROeX
We call make_edge_scal to compute \(\rho^{'\edge}, X_k^{\edge}, (\rho h)^{'\edge}\), and \(\Ubt^{\edge}\) 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'\) and \(X_k\) to faces, and the choice of energy prediction is as follows:
For enthalpy_pred_type = 1, we predict \((\rho h)'\) 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'\) to faces.
For enthalpy_pred_type = 6, we predict \(T'\) to faces.
We are using enthalpy_pred_type = 1 for now. The equations of motion are:
2D Cartesian Case
Predict \(s\) to r-faces using a 1D extrapolation.
Predict \(s\) to x-faces using a full-dimensional extrapolation.
Predict \(s\) to x-faces using a 1D extrapolation.
Predict \(s\) to r-faces using a full-dimensional extrapolation.
Predict \(s\) to r-faces using a 1D extrapolation:
(195)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_r} &=& s_{\ib-\eb_r}^n + \left(\half - \frac{\dt}{2h}w_{\ib-\half\eb_r}^{\mac}\right)\Delta_r s_{\ib-\eb_r}^n, \\ s_{R,\ib-\half\eb_r} &=& s_{\ib} - \left(\half + \frac{\dt}{2h}w_{\ib-\half\eb_r}^{\mac}\right)\Delta_r s_{\ib}^n.\end{aligned}\end{split}\]Upwind based on \(w^{\mac}\):
(196)\[\begin{split}s_{\ib-\half\eb_r} = \begin{cases} \half\left(s_{L,\ib-\half\eb_r} + s_{R,\ib-\half\eb_r}\right), & \left|w^{\mac}_{\ib-\half\eb_r}\right| < \epsilon \\ s_{L,\ib-\half\eb_r}, & w^{\mac}_{\ib-\half\eb_r} > 0, \\ s_{R,\ib-\half\eb_r}, & w^{\mac}_{\ib-\half\eb_r} < 0. \\ \end{cases}\end{split}\]Predict \(s\) to x-faces using a full-dimensional extrapolation. First, the normal derivative and forcing terms:
(197)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{\edge} &=& s_{\ib-\eb_x}^n + \left(\half - \frac{\dt}{2h}u_{\ib-\half\eb_x}^{\mac}\right)\Delta_x s_{\ib-\eb_x}^n + \frac{\dt}{2}f_{\ib-\eb_x}^n \\ s_{R,\ib-\half\eb_x}^{\edge} &=& s_{\ib}^n - \left(\half + \frac{\dt}{2h}u_{\ib-\half\eb_x}^{\mac}\right)\Delta_x s_{\ib}^n + \frac{\dt}{2}f_{\ib}^n \end{aligned}\end{split}\]Account for the transverse terms:
if is_conservative then
(198)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{\edge} &=& s_{L,\ib-\half\eb_x}^{\edge} - \frac{\dt}{2h}\left[\left(w^{\mac}s\right)_{\ib-\eb_x+\half\eb_r} - \left(w^{\mac}s\right)_{\ib-\eb_x-\half\eb_r}\right] - \frac{\dt}{2h}s_{\ib-\eb_x}^{n}\left(u_{\ib-\half\eb_x}^{\mac}-u_{\ib-\frac{3}{2}\eb_x}^{\mac}\right)\nonumber \\ &&\\ s_{R,\ib-\half\eb_x}^{\edge} &=& s_{R,\ib-\half\eb_x}^{\edge} - \frac{\dt}{2h}\left[\left(w^{\mac}s\right)_{\ib+\half\eb_r} - \left(w^{\mac}s\right)_{\ib-\half\eb_r}\right] - \frac{\dt}{2h}s_{\ib}^{n}\left(u_{\ib+\half\eb_x}^{\mac}-u_{\ib-\half\eb_x}^{\mac}\right)\end{aligned}\end{split}\]else
(199)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{\edge} &=& s_{L,\ib-\half\eb_x}^{\edge} - \frac{\dt}{4h}\left(w^{\mac}_{\ib-\eb_x+\half\eb_r} + w^{\mac}_{\ib-\eb_x-\half\eb_r}\right)\left(s_{\ib-\eb_x+\half\eb_r} - s_{\ib-\eb_x-\half\eb_r}\right)\\ s_{R,\ib-\half\eb_x}^{\edge} &=& s_{R,\ib-\half\eb_x}^{\edge} - \frac{\dt}{4h}\left(w^{\mac}_{\ib+\half\eb_r} + w^{\mac}_{\ib-\half\eb_r}\right)\left(s_{\ib+\half\eb_r} - s_{\ib-\half\eb_r}\right)\end{aligned}\end{split}\]end if
Account for the \(\partial w_0/\partial r\) term:
if is_vel and comp = 2 then
(200)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{\edge} &=& s_{L,\ib-\half\eb_x}^{\edge} - \frac{\dt}{4h}\left(\wt^{\mac}_{\ib-\eb_x+\half\eb_r} + \wt^{\mac}_{\ib-\eb_x-\half\eb_r}\right)\left(w_{0,\ib+\half\eb_r}-w_{0,\ib-\half\eb_r}\right) \\ s_{R,\ib-\half\eb_x}^{\edge} &=& s_{R,\ib-\half\eb_x}^{\edge} - \frac{\dt}{4h}\left(\wt^{\mac}_{\ib+\half\eb_r} + \wt^{\mac}_{\ib-\half\eb_r}\right)\left(w_{0,\ib+\half\eb_r}-w_{0,\ib-\half\eb_r}\right) \\\end{aligned}\end{split}\]end if
Upwind based on \(u^{\mac}\).
(201)\[\begin{split}s_{\ib-\half\eb_x}^{\edge} = \begin{cases} \half\left(s_{L,\ib-\half\eb_x}^{\edge} + s_{R,\ib-\half\eb_x}^{\edge}\right), & \left|u^{\mac}_{\ib-\half\eb_x}\right| < \epsilon \\ s_{L,\ib-\half\eb_x}^{\edge}, & u^{\mac}_{\ib-\half\eb_x} > 0, \\ s_{R,\ib-\half\eb_x}^{\edge}, & u^{\mac}_{\ib-\half\eb_x} < 0. \end{cases}\end{split}\]Predict \(s\) to x-faces using a 1D extrapolation:
(202)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x} &=& s_{\ib-\eb_x}^n + \left(\half - \frac{\dt}{2h}u_{\ib-\half\eb_x}^{\mac}\right)\Delta_x s_{\ib-\eb_x}^n, \\ s_{R,\ib-\half\eb_x} &=& s_{\ib} - \left(\half + \frac{\dt}{2h}u_{\ib-\half\eb_x}^{\mac}\right)\Delta_x s_{\ib}^n.\end{aligned}\end{split}\]Upwind based on \(u^{\mac}\):
(203)\[\begin{split}s_{\ib-\half\eb_x} = \begin{cases} \half\left(s_{L,\ib-\half\eb_x} + s_{R,\ib-\half\eb_x}\right), & \left|u^{\mac}_{\ib-\half\eb_x}\right| < \epsilon \\ s_{L,\ib-\half\eb_x}, & u^{\mac}_{\ib-\half\eb_x} > 0, \\ s_{R,\ib-\half\eb_x}, & u^{\mac}_{\ib-\half\eb_x} < 0. \\ \end{cases}\end{split}\]Predict \(s\) to r-faces using a full-dimensional extrapolation. First, the normal derivative and forcing terms:
(204)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_r}^{\edge} &=& s_{\ib-\eb_r}^n + \left(\half - \frac{\dt}{2h}w_{\ib-\half\eb_r}^{\mac}\right)\Delta_r s_{\ib-\eb_r}^n + \frac{\dt}{2}f_{\ib-\eb_r}^n \\ s_{R,\ib-\half\eb_r}^{\edge} &=& s_{\ib}^n - \left(\half + \frac{\dt}{2h}w_{\ib-\half\eb_r}^{\mac}\right)\Delta_r s_{\ib}^n + \frac{\dt}{2}f_{\ib}^n \end{aligned}\end{split}\]Account for the transverse terms: if is_conservative then
(205)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_r}^{\edge} &=& s_{L,\ib-\half\eb_r}^{\edge} - \frac{\dt}{2h}\left[\left(u^{\mac}s\right)_{\ib-\eb_r+\half\eb_x} - \left(u^{\mac}s\right)_{\ib-\eb_r-\half\eb_x}\right] - \frac{\dt}{2h}s_{\ib-\eb_r}^{n}\left(w_{\ib-\half\eb_r}^{\mac}-w_{\ib-\frac{3}{2}\eb_r}^{\mac}\right)\nonumber\\ && \\ s_{R,\ib-\half\eb_r}^{\edge} &=& s_{R,\ib-\half\eb_r}^{\edge} - \frac{\dt}{2h}\left[\left(u^{\mac}s\right)_{\ib+\half\eb_x} - \left(u^{\mac}s\right)_{\ib-\half\eb_x}\right] - \frac{\dt}{2h}s_{\ib}^{n}\left(w_{\ib+\half\eb_r}^{\mac}-w_{\ib-\half\eb_r}^{\mac}\right)\end{aligned}\end{split}\]else
(206)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_r}^{\edge} &=& s_{L,\ib-\half\eb_r}^{\edge} - \frac{\dt}{4h}\left(u^{\mac}_{\ib-\eb_r+\half\eb_x} + u^{\mac}_{\ib-\eb_r-\half\eb_x}\right)\left(s_{\ib-\eb_r+\half\eb_x} - s_{\ib-\eb_r-\half\eb_x}\right)\\ s_{R,\ib-\half\eb_r}^{\edge} &=& s_{R,\ib-\half\eb_r}^{\edge} - \frac{\dt}{4h}\left(u^{\mac}_{\ib+\half\eb_x} + u^{\mac}_{\ib-\half\eb_x}\right)\left(s_{\ib+\half\eb_x} - s_{\ib-\half\eb_x}\right)\end{aligned}\end{split}\]end if
Account for the \(\partial w_0/\partial r\) term: if is_vel and comp = 2 then
(207)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_r}^{\edge} &=& s_{L,\ib-\half\eb_r}^{\edge} - \frac{\dt}{4h}\left(\wt^{\mac}_{\ib-\half\eb_r} + \wt^{\mac}_{\ib-\frac{3}{2}\eb_r}\right)\left(w_{0,\ib-\half\eb_r}-w_{0,\ib-\frac{3}{2}\eb_r}\right) \\ s_{R,\ib-\half\eb_r}^{\edge} &=& s_{R,\ib-\half\eb_r}^{\edge} - \frac{\dt}{4h}\left(\wt^{\mac}_{\ib+\half\eb_r} + \wt^{\mac}_{\ib-\half\eb_r}\right)\left(w_{0,\ib+\half\eb_r}-w_{0,\ib-\half\eb_r}\right) \\\end{aligned}\end{split}\]end if
Upwind based on \(w^{\mac}\):
(208)\[\begin{split}s_{\ib-\half\eb_r} = \begin{cases} \half\left(s_{L,\ib-\half\eb_r} + s_{R,\ib-\half\eb_r}\right), & \left|w^{\mac}_{\ib-\half\eb_r}\right| < \epsilon \\ u_{L,\ib-\half\eb_r}, & w^{\mac}_{\ib-\half\eb_r} > 0, \\ u_{R,\ib-\half\eb_r}, & w^{\mac}_{\ib-\half\eb_r} < 0. \\ \end{cases}\end{split}\]
3D Cartesian Case
This algorithm is more complicated than the 2D case since we include the effects of corner coupling.
Predict \(s\) to x-faces using a 1D extrapolation.
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.
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.
Predict \(s\) to y-faces using a full-dimensional extrapolation.
Predict \(s\) to r-faces using a full-dimensional extrapolation.
Predict \(s\) to x-faces using a 1D extrapolation.
(209)\[s_{L,\ib-\half\eb_x} &=& s_{\ib-\eb_x}^n + \left(\half - \frac{\dt}{2h}u_{\ib-\half\eb_x}^{\mac}\right)\Delta_x s_{\ib-\eb_x}^n,\](210)\[s_{R,\ib-\half\eb_x} &=& s_{\ib} - \left(\half + \frac{\dt}{2h}u_{\ib-\half\eb_x}^{\mac}\right)\Delta_x s_{\ib}^n.\]Upwind based on \(u^{\mac}\):
(211)\[\begin{split}s_{\ib-\half\eb_x} = \begin{cases} \half\left(s_{L,\ib-\half\eb_x} + s_{R,\ib-\half\eb_x}\right), & \left|u^{\mac}_{\ib-\half\eb_x}\right| < \epsilon \\ s_{L,\ib-\half\eb_x}, & u^{\mac}_{\ib-\half\eb_x} > 0, \\ s_{R,\ib-\half\eb_x}, & u^{\mac}_{\ib-\half\eb_x} < 0. \\ \end{cases}\end{split}\]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_{\ib-\half\eb_x}^{x|y}\) means “state \(s_{\ib-\half\eb_x}\) that has been updated to account for the transverse derivatives in the \(y\)-direction”.
if is_conservative then
(212)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{x|y} &=& s_{L,\ib-\half\eb_x} - \frac{\dt}{3h}\left[(sv^{\mac})_{\ib-\eb_x+\half\eb_y}-(sv^{\mac})_{\ib-\eb_x-\half\eb_y}\right], \\ s_{R,\ib-\half\eb_x}^{x|y} &=& s_{R,\ib-\half\eb_x} - \frac{\dt}{3h}\left[(sv^{\mac})_{\ib+\half\eb_y}-(sv^{\mac})_{\ib-\half\eb_y}\right].\end{aligned}\end{split}\]else
(213)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{x|y} &=& s_{L,\ib-\half\eb_x} - \frac{\dt}{6h}\left(v_{\ib-\eb_x+\half\eb_y}^{\mac} + v_{\ib-\eb_x-\half\eb_y}^{\mac}\right)\left(s_{\ib-\eb_x+\half\eb_y} - s_{\ib-\eb_x-\half\eb_y}\right), \\ s_{R,\ib-\half\eb_x}^{x|y} &=& s_{R,\ib-\half\eb_x} - \frac{\dt}{6h}\left(v_{\ib+\half\eb_y}^{\mac} + v_{\ib-\half\eb_y}^{\mac}\right)\left(s_{\ib+\half\eb_y} - s_{\ib-\half\eb_y}\right).\end{aligned}\end{split}\]end if
Upwind based on \(u^{\mac}\):
(214)\[\begin{split}s_{\ib-\half\eb_x}^{x|y} = \begin{cases} \half\left(s_{L,\ib-\half\eb_x}^{x|y} + s_{R,\ib-\half\eb_x}^{x|y}\right), & \left|u^{\mac}_{\ib-\half\eb_x}\right| < \epsilon \\ s_{L,\ib-\half\eb_x}^{x|y}, & u^{\mac}_{\ib-\half\eb_x} > 0, \\ s_{R,\ib-\half\eb_x}^{x|y}, & u^{\mac}_{\ib-\half\eb_x} < 0. \end{cases}\end{split}\]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{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{\edge} = s_{L,\ib-\half\eb_x} &-& \frac{\dt}{2h}\left[(s^{y|r}v^{\mac})_{\ib-\eb_x+\half\eb_y}-({s^{y|r}v^{\mac})_{\ib-\eb_x-\half\eb_y}}\right] \nonumber \\ &-& \frac{\dt}{2h}\left[(s^{r|y}w^{\mac})_{\ib-\eb_x+\half\eb_r}-({s^{r|y}w^{\mac})_{\ib-\eb_x-\half\eb_r}}\right] \nonumber \\ &-& \frac{\dt}{2h}s_{\ib-\eb_x}\left(u_{\ib-\half\eb_x}^{\mac}-u_{\ib-\frac{3}{2}\eb_x}^{\mac}\right) + \frac{\dt}{2}f_{\ib-\eb_x}, \\ s_{R,\ib-\half\eb_x}^{\edge} = s_{R,\ib-\half\eb_x} &-& \frac{\dt}{2h}\left[(s^{y|r}v^{\mac})_{\ib+\half\eb_y}-({s^{y|r}v^{\mac})_{\ib-\half\eb_y}}\right] \nonumber \\ &-& \frac{\dt}{2h}\left[(s^{r|y}w^{\mac})_{\ib+\half\eb_r}-({s^{r|y}w^{\mac})_{\ib-\half\eb_r}}\right] \nonumber \\ &-& \frac{\dt}{2h}s_{\ib}\left(u_{\ib+\half\eb_x}^{\mac}-u_{\ib-\half\eb_x}^{\mac}\right) + \frac{\dt}{2}f_{\ib}.\end{aligned}\end{split}\]else
(216)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{\edge} = s_{L,\ib-\half\eb_x} &-& \frac{\dt}{4h}\left(v_{\ib-\eb_x+\half\eb_y}^{\mac}+v_{\ib-\eb_x-\half\eb_y}^{\mac}\right)\left(s_{\ib-\eb_x+\half\eb_y}^{y|r}-s_{\ib-\eb_x-\half\eb_y}^{y|r}\right) \nonumber \\ &-& \frac{\dt}{4h}\left(w_{\ib-\eb_x+\half\eb_r}^{\mac}+w_{\ib-\eb_x-\half\eb_r}^{\mac}\right)\left(s_{\ib-\eb_x+\half\eb_r}^{r|y}-s_{\ib-\eb_x-\half\eb_r}^{r|y}\right) + \frac{\dt}{2}f_{\ib-\eb_x}, \nonumber \\ && \\ s_{R,\ib-\half\eb_x}^{\edge} = s_{R,\ib-\half\eb_x} &-& \frac{\dt}{4h}\left(v_{\ib+\half\eb_y}^{\mac}+v_{\ib-\half\eb_y}^{\mac}\right)\left(s_{\ib+\half\eb_y}^{y|r}-s_{\ib-\half\eb_y}^{y|r}\right) \nonumber \\ &-& \frac{\dt}{4h}\left(w_{\ib+\half\eb_r}^{\mac}+w_{\ib-\half\eb_r}^{\mac}\right)\left(s_{\ib+\half\eb_r}^{r|y}-s_{\ib-\half\eb_r}^{r|y}\right) + \frac{\dt}{2}f_{\ib}.\end{aligned}\end{split}\]end if
Account for the \(\partial w_0/\partial r\) term:
if is_vel and comp = 2 then
(217)\[\begin{split}\begin{aligned} s_{L,\ib-\half\eb_x}^{\edge} &=& s_{L,\ib-\half\eb_x}^{\edge} - \frac{\dt}{4h}\left(\wt^{\mac}_{\ib-\eb_x+\half\eb_r} + \wt^{\mac}_{\ib-\eb_x-\half\eb_r}\right)\left(w_{0,\ib+\half\eb_r}-w_{0,\ib-\half\eb_r}\right) \\ s_{R,\ib-\half\eb_x}^{\edge} &=& s_{R,\ib-\half\eb_x}^{\edge} - \frac{\dt}{4h}\left(\wt^{\mac}_{\ib+\half\eb_r} + \wt^{\mac}_{\ib-\half\eb_r}\right)\left(w_{0,\ib+\half\eb_r}-w_{0,\ib-\half\eb_r}\right) \\\end{aligned}\end{split}\]end if
Upwind based on \(u^{\mac}\):
(218)\[\begin{split}s_{\ib-\half\eb_x}^{\edge} = \begin{cases} \half\left(s_{L,\ib-\half\eb_x}^{\edge} + s_{R,\ib-\half\eb_x}^{\edge}\right), & \left|u^{\mac}_{\ib-\half\eb_x}\right| < \epsilon \\ s_{L,\ib-\half\eb_x}^{\edge}, & u^{\mac}_{\ib-\half\eb_x} > 0, \\ s_{R,\ib-\half\eb_x}^{\edge}, & u^{\mac}_{\ib-\half\eb_x} < 0. \end{cases}\end{split}\]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 \(\Ub^{\mac,*}\) 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,\ib-\half\eb_x}^{1D} = u_{\ib-\eb_x} + \left[\half - \frac{\dt}{2h}{\rm max}(0,u_{\ib-\eb_x})\right]\Delta_xu_{\ib-\eb_x},\](220)\[u_{R,\ib-\half\eb_x}^{1D} = u_{\ib} + \left[\half - \frac{\dt}{2h}{\rm min}(0,u_{\ib})\right]\Delta_xu_{\ib} .\](221)\[v_{L,\ib-\half\eb_x}^{1D} = v_{\ib-\eb_x} + \left[\half - \frac{\dt}{2h}{\rm max}(0,v_{\ib-\eb_x})\right]\Delta_xv_{\ib-\eb_x},\](222)\[v_{R,\ib-\half\eb_x}^{1D} &=& v_{\ib} + \left[\half - \frac{\dt}{2h}{\rm min}(0,v_{\ib})\right]\Delta_xv_{\ib}.\]We obtain the normal velocity using the Riemann problem:
(223)\[\begin{split}u_{\ib-\half\eb_x}^{1D} = \begin{cases} 0, & \left(u_{L,\ib-\half\eb_x}^{1D} \le 0 ~~ {\rm AND} ~~ u_{R,\ib-\half\eb_x}^{1D} \ge 0\right) ~~ {\rm OR} ~~ \left|u_{L,\ib-\half\eb_x}^{1D} + u_{R,\ib-\half\eb_x}^{1D}\right| < \epsilon, \\ u_{L,\ib-\half\eb_x}^{1D}, & u_{L,\ib-\half\eb_x}^{1D} + u_{R,\ib-\half\eb_x}^{1D} > 0, \\ u_{R,\ib-\half\eb_x}^{1D}, & u_{L,\ib-\half\eb_x}^{1D} + u_{R,\ib-\half\eb_x}^{1D} < 0. \end{cases}\end{split}\]We obtain the transverse velocity by upwinding based on \(u_{\ib-\half\eb_x}^{1D}\):
(224)\[\begin{split}v_{\ib-\half\eb_x}^{1D} = \begin{cases} \half\left(v_{L,\ib-\half\eb_x}^{1D} + v_{R,\ib-\half\eb_x}^{1D}\right), & \left|u_{\ib-\half\eb_x}^{1D}\right| < \epsilon \\ v_{L,\ib-\half\eb_x}^{1D}, & u_{\ib-\half\eb_x}^{1D} > 0, \\ v_{R,\ib-\half\eb_x}^{1D}, & u_{\ib-\half\eb_x}^{1D} < 0. \end{cases}\end{split}\]We perform analogous operations to compute both components of velocity at the y-faces, \(\Ub_{\ib-\half\eb_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{split}\begin{aligned} u_{L,\ib-\half\eb_x}^{\mac,*} &=& u_{L,\ib-\half\eb_x}^{1D} - \frac{\dt}{4h}\left(v_{\ib-\eb_x+\half\eb_y}^{1D}+v_{\ib-\eb_x-\half\eb_y}^{1D}\right)\left(u_{\ib-\eb_x+\half\eb_y}^{1D} - u_{\ib-\eb_x-\half\eb_y}^{1D}\right) + \frac{\dt}{2}f_{u,\ib-\eb_x}, \\ u_{R,\ib-\half\eb_x}^{\mac,*} &=& u_{R,\ib-\half\eb_x}^{1D} - \frac{\dt}{4h}\left(v_{\ib+\half\eb_y}^{1D}+v_{\ib-\half\eb_y}^{1D}\right)\left(u_{\ib+\half\eb_y}^{1D} - u_{\ib-\half\eb_y}^{1D}\right) + \frac{\dt}{2}f_{u,\ib}.\end{aligned}\end{split}\]Then we solve a Riemann problem:
(226)\[\begin{split}u_{\ib-\half\eb_x}^{\mac,*} = \begin{cases} 0, & \left(u_{L,\ib-\half\eb_x}^{\mac,*} \le 0 ~~ {\rm AND} ~~ u_{R,\ib-\half\eb_x}^{\mac,*} \ge 0\right) ~~ {\rm OR} ~~ \left|u_{L,\ib-\half\eb_x}^{\mac,*} + u_{R,\ib-\half\eb_x}^{\mac,*}\right| < \epsilon, \\ u_{L,\ib-\half\eb_x}^{\mac,*}, & u_{L,\ib-\half\eb_x}^{\mac,*} + u_{R,\ib-\half\eb_x}^{\mac,*} > 0, \\ u_{R,\ib-\half\eb_x}^{\mac,*}, & u_{L,\ib-\half\eb_x}^{\mac,*} + u_{R,\ib-\half\eb_x}^{\mac,*} < 0. \end{cases}\end{split}\]We perform analogous operations to compute the normal velocity at the y-faces, \(v^{\mac,*}_{\ib-\half\eb_y}\).
3D Cartesian Case
This is more complicated than the 2D case because we include corner
coupling. We compute \(\Ub_{\ib-\half\eb_x}^{1D},
\Ub_{\ib-\half\eb_y}^{1D}\), and \(\Ub_{\ib-\half\eb_z}^{1D}\) in
an analogous manner as varden U_L^1D]-Eq.224. Then we compute an intermediate state,
\(u_{\ib-\half\eb_y}^{y|z}\), which is described as “state
\(u_{\ib-\half\eb_y}^{1D}\) that has been updated to account for
the transverse derivatives in the z direction”, using:
Then upwind based on \(v_{\ib-\half\eb_y}^{1D}\):
We use an analogous procedure to compute five more intermediate states, \(u_{\ib-\half\eb_z}^{z|y}, v_{\ib-\half\eb_x}^{x|z}, v_{\ib-\half\eb_z}^{z|x}, w_{\ib-\half\eb_x}^{x|y}\), and \(w_{\ib-\half\eb_y}^{y|x}\). Then we do a full-dimensional extrapolation to get the MAC velocities at normal faces:
Then we use the Riemann solver given above for the 2D case (Eq.226) to compute \(u_{\ib-\half\eb_x}^{\mac,*}\). We use an analogous procedure to obtain \(v_{\ib-\half\eb_y}^{\mac,*}\) and \(w_{\ib-\half\eb_z}^{\mac,*}\).
Computing \(\Ub^{\edge}\) and \(\rho^{\edge}\) in VARDEN
To compute \(\Ub^{\edge}\), VARDEN uses the exact same algorithm as the \(s^{\edge}\) case in MAESTROeX. The algorithm for \(\rho^{\edge}\) in VARDEN is slightly different than the \(s^{\edge}\) 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:
Then we upwind based on \(u^{\mac}\):
We use an analogous procedure to obtain \(s_{\ib-\half\eb_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:
The extrapolation of a “conserved” \(s\) to x-faces is:
Then we upwind based on \(u^{\mac}\).
We use an analogous procedure to compute \(s_{\ib-\half\eb_y}^{\edge}\).
3D Cartesian Case
This is more complicated than the 2D case because we include corner coupling. We first compute \(s_{\ib-\half\eb_x}^{1D}\), \(s_{\ib-\half\eb_y}^{1D}\), and \(s_{\ib-\half\eb_z}^{1D}\) in an analogous manner to Eq.230 and Eq.231. Then we compute six intermediate states, \(s_{\ib-\half\eb_x}^{x|y}, s_{\ib-\half\eb_x}^{x|z}, s_{\ib-\half\eb_y}^{y|x}, s_{\ib-\half\eb_y}^{y|z}, s_{\ib-\half\eb_z}^{z|x}\), and \(s_{\ib-\half\eb_z}^{z|y}\). For the “non-conservative case”, we use, for example:
For the “conservative” case, we use, for example:
Then we upwind based on \(u^{\mac}\):
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:
The extrapolation of a “conserved” \(s\) to x-faces is:
Then we upwind based on \(u^{\mac}\), as in Eq.235. We use an analogous procedure to compute both \(s_{\ib-\half\eb_y}^{\edge}\) and \(s_{\ib-\half\eb_z}\).
ESTATE_FPU in GODUNOV_2D/3D.f
First, the normal predictor.
(241)\[\begin{split}\begin{aligned} s_L^x &=& s_{\ib-\eb_x} + \left(\half - \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib-\eb_x} + \underbrace{\frac{\dt}{2}\text{TFORCES}_{\ib-\eb_x}}_{\text{IF USE\_MINION}} \\ s_R^x &=& s_{\ib} - \left(\half + \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib} + \underbrace{\frac{\dt}{2}\text{TFORCES}_{\ib}}_{\text{IF USE\_MINION}}\end{aligned}\end{split}\]If USE_MINION and ICONSERVE then:
(242)\[\begin{split}\begin{aligned} s_L^x &=& s_L^x - \frac{\dt}{2}s_{\ib-\eb_x}\text{DIVU}_{\ib-\eb_x} \\ s_R^x &=& s_R^x - \frac{\dt}{2}s_{\ib}\text{DIVU}_{\ib}\end{aligned}\end{split}\]Apply boundary conditions on \(s_L^x\) and \(s_R^x\). Then,
(243)\[\begin{split}\text{s}_{\ib-\half\eb_x}^x = \begin{cases} s_L^x, & \text{UEDGE}_{\ib-\half\eb_x} > 0, \\ s_R^x, & \text{else}. \\ \end{cases}\end{split}\]Then, if \(|\text{UEDGE}_{\ib-\half\eb_x}| \le \epsilon\), we set \(s_{\ib-\half\eb_x}^x = (s_L^x+s_R^x)/2\). The procedure to obtain \(s_{\ib-\half\eb_y}^y\) is analogous.
Now, the transverse terms.
If ICONSERVE then:
(244)\[\begin{split}\begin{aligned} \text{sedge}_L^x &=& s_{\ib-\eb_x} + \left(\half - \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib-\eb_x} + \frac{\dt}{2}\text{TFORCES}_{\ib-\eb_x} \nonumber\\ && - \frac{\dt}{2}\left[\frac{\text{VEDGE}_{\ib-\eb_x+\half\eb_y}s_{\ib-\eb_x+\half\eb_y}^y - \text{VEDGE}_{\ib-\eb_x-\half\eb_y}s_{\ib-\eb_x-\half\eb_y}^y}{h_y}\right.\nonumber\\ && ~~~~~~~~~~ \left. - \frac{s_{\ib-\eb_x}(\text{VEDGE}_{\ib-\eb_x+\half\eb_y}-\text{VEDGE}_{\ib-\eb_x-\half\eb_y})}{h_y}+s_{\ib-\eb_x}\text{DIVU}_{\ib-\eb_x}\right]\\ \text{sedge}_R^x &=& s_{\ib} - \left(\half + \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib} + \frac{\dt}{2}\text{TFORCES}_{\ib} \nonumber\\ && - \frac{\dt}{2}\left[\frac{\text{VEDGE}_{\ib+\half\eb_y}s_{\ib+\half\eb_y}^y - \text{VEDGE}_{\ib-\half\eb_y}s_{\ib-\half\eb_y}^y}{h_y}\right.\nonumber\\ && ~~~~~~~~~~ \left. - \frac{s_{\ib}(\text{VEDGE}_{\ib+\half\eb_y}-\text{VEDGE}_{\ib-\half\eb_y})}{h_y}+s_{\ib}\text{DIVU}_{\ib}\right]\end{aligned}\end{split}\]Now, define \(\text{VBAR}_{\ib} = (\text{VEDGE}_{\ib+\half\eb_y}+\text{VEDGE}_{\ib-\half\eb_y})/2\).
If NOT ICONSERVE and \(\text{VEDGE}_{\ib+\half\eb_y}\cdot\text{VEDGE}_{\ib-\half\eb_y} < 0\) and \(\text{VBAR}_{\ib} < 0\) then:
(245)\[\begin{split}\begin{align} \text{sedge}_L^x = s_{\ib-\eb_x} &+ \left(\half - \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib} + \frac{\dt}{2}\text{TFORCES}_{\ib-\eb_x} \nonumber\\ & - \frac{\dt}{2}\left[\frac{\text{VBAR}_{\ib-\eb_x}(s_{\ib-\eb_x+\eb_y}-s_{\ib-\eb_x})}{h_y}\right] \end{align}\end{split}\](246)\[\begin{split}\begin{align} \text{sedge}_R^x = s_{\ib} &- \left(\half + \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib} + \frac{\dt}{2}\text{TFORCES}_{\ib} \nonumber\\ & - \frac{\dt}{2}\left[\frac{\text{VBAR}_{\ib}(s_{\ib+\eb_y}-s_{\ib})}{h_y}\right] \end{align}\end{split}\]Else If NOT ICONSERVE and \(\text{VEDGE}_{\ib+\half\eb_y}\cdot\text{VEDGE}_{\ib-\half\eb_y} < 0\) and \(\text{VBAR}_{\ib} \ge 0\) then:
(247)\[\begin{split}\begin{aligned} \text{sedge}_L^x = s_{\ib-\eb_x} &+& \left(\half - \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib-\eb_x} + \frac{\dt}{2}\text{TFORCES}_{\ib-\eb_x} \nonumber\\ && - \frac{\dt}{2}\left[\frac{\text{VBAR}_{\ib-\eb_x}(s_{\ib-\eb_x}-s_{\ib-\eb_x-\eb_y})}{h_y}\right] \\ \text{sedge}_R^x = s_{\ib} &-& \left(\half + \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib} + \frac{\dt}{2}\text{TFORCES}_{\ib} \nonumber\\ && - \frac{\dt}{2}\left[\frac{\text{VBAR}_{\ib}(s_{\ib}-s_{\ib-\eb_y})}{h_y}\right]\end{aligned}\end{split}\]Else If NOT ICONSERVE and \(\text{VEDGE}_{\ib+\half\eb_y}\cdot\text{VEDGE}_{\ib-\half\eb_y} \ge 0\) then:
(248)\[\begin{split}\begin{align} \text{sedge}_L^x &= s_{\ib-\eb_x} + \left(\half - \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib-\eb_x} + \frac{\dt}{2}\text{TFORCES}_{\ib-\eb_x} \nonumber\\ && - \frac{\dt}{2}\left[\frac{(\text{VEDGE}_{\ib-\eb_x+\half\eb_y}+\text{VEDGE}_{\ib-\eb_x-\half\eb_y})(s_{\ib-\eb_x+\half\eb_y}-s_{\ib-\eb_x-\half\eb_y})}{2h_y}\right] \\ \text{sedge}_R^x &=& s_{\ib} - \left(\half + \frac{\dt}{h_x}\text{UEDGE}_{\ib-\half\eb_x}\right)\Delta^x s_{\ib} + \frac{\dt}{2}\text{TFORCES}_{\ib} \nonumber\\ && - \frac{\dt}{2}\left[\frac{(\text{VEDGE}_{\ib+\half\eb_y}+\text{VEDGE}_{\ib-\half\eb_y})(s_{\ib+\half\eb_y}-s_{\ib-\half\eb_y})}{2h_y}\right] \end{align}\end{split}\]Finally, upwind analogous to Eq.243 to get \(\text{sedge}_{\ib-\half\eb_x}\).
ESTATE in GODUNOV_2D/3D.f
First, the normal predictor.
If USE_MINION then:
Apply boundary conditions on \(s_L^x\) and \(s_R^x\). Then,
Then, if \(|\text{UAD}_{\ib-\half\eb_x}| \le \epsilon\), we set \(s_{\ib-\half\eb_x}^x = (s_L^x+s_R^x)/2\).
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}_{\ib-\half\eb_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,\ib-\myhalf\eb_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 (\(\ib+\eb_x \rightarrow 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\left[s_{i,+} - s_{i,-} + s_{6,i}(1-\xi)\right],\label{Quadratic Interp}\](254)\[s_{6,i}= 6s_{i} - 3\left(s_{i,-}+s_{i,+}\right),\](255)\[\xi = \frac{x - (i-\myhalf)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|\Delta t/h\):
(256)\[\begin{split}\begin{aligned} \mathcal{I}_{i,+}(\sigma) &=& \frac{1}{\sigma h}\int_{(i+\myhalf)h-\sigma h}^{(i+\myhalf)h}s_i^I(x)dx \\ \mathcal{I}_{i,-}(\sigma) &=& \frac{1}{\sigma h}\int_{(i-\myhalf)h}^{(i-\myhalf)h+\sigma h}s_i^I(x)dx\end{aligned}\end{split}\]Plugging in ([Quadratic Interp]) gives:
(257)\[\begin{split}\begin{aligned} \mathcal{I}_{i,+}(\sigma) &=& s_{j,+} - \frac{\sigma}{2}\left[s_{j,+}-s_{j,-}-\left(1-\frac{2}{3}\sigma\right)s_{6,i}\right], \\ \mathcal{I}_{i,-}(\sigma) &=& s_{j,-} + \frac{\sigma}{2}\left[s_{j,+}-s_{j,-}+\left(1-\frac{2}{3}\sigma\right)s_{6,i}\right].\end{aligned}\end{split}\]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{split}\begin{aligned} s_{L,i-\myhalf} &=& \begin{cases} \mathcal{I}_{i-1,+}(\sigma), & u_{i-1} ~ \text{or} ~ u_{i-\myhalf}^{\mac} > 0 \\ s_{i-1}, & \text{else}. \end{cases}\\ s_{R,i-\myhalf} &=& \begin{cases} \mathcal{I}_{i,-}(\sigma), & u_{i} ~ \text{or} ~ u_{i-\myhalf}^{\mac} < 0 \\ s_{i}, & \text{else}. \end{cases}\end{aligned}\end{split}\]
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:
A more compact way of writing this is
Without the limiters, Eq.259 is the familiar 4th-order spatial interpolation formula:
Next, we must ensure that \(s_{i+\myhalf}\) lies between the adjacent cell-centered values:
In anticipation of further limiting, we set double-valued face-centered values:
Modify \(s_{i,\pm}\) using a quadratic limiter. First, we test whether \(s_i\) is a local extreumum with the condition:
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:
Colella and Sekora Based Approach
Spatially interpolate \(s\) to edges.
Use a 4th-order interpolation in space to obtain edge values:
(269)\[s_{i+\myhalf} = \frac{7}{12}\left(s_{i+1}+s_i\right) - \frac{1}{12}\left(s_{i+2}+s_{i-1}\right).\]Then, if \((s_{i+\myhalf}-s_i)(s_{i+1}-s_{i+\myhalf}) < 0\), we limit \(s_{i+\myhalf}\) using a nonlinear combination of approximations to the second derivative. First, define:
(270)\[\begin{split}\begin{aligned} (D^2s)_{i+\myhalf} &=& 3\left(s_{i}-2s_{i+\myhalf}+s_{i+1}\right) \\ (D^2s)_{i+\myhalf,L} &=& s_{i-1}-2s_{i}+s_{i+1} \\ (D^2s)_{i+\myhalf,R} &=& s_{i}-2s_{i+1}+s_{i+2}\end{aligned}\end{split}\]Then, define
(271)\[s = \text{sign}\left[(D^2s)_{i+\myhalf}\right],\](272)\[(D^2s)_{i+\myhalf,\text{lim}} = s\max\left\{\min\left[Cs(D^2s)_{i+\myhalf,L},Cs(D^2s)_{i+\myhalf,R},s(D^2s)_{i+\myhalf}\right],0\right\},\]where \(C=1.25\) was used in Colella and Sekora. Then,
(273)\[s_{i+\myhalf} = \frac{1}{2}\left(s_{i}+s_{i+1}\right) - \frac{1}{6}(D^2s)_{i+\myhalf,\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\myhalf} - 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{split}\begin{aligned} (Ds)_{i,{\rm face},-} &=& s_{i-\myhalf} - s_{i-\sfrac{3}{2}} \\ (Ds)_{i,{\rm face},+} &=& s_{i+\sfrac{3}{2}} - s_{i+\myhalf}\end{aligned}\end{split}\](276)\[(Ds)_{i,{\rm face,min}} = \min\left[\left|(Ds)_{i,{\rm face},-}\right|,\left|(Ds)_{i,{\rm face},+}\right|\right].\](277)\[\begin{split}\begin{aligned} (Ds)_{i,{\rm cc},-} &=& s_{i} - s_{i-1} \\ (Ds)_{i,{\rm cc},+} &=& s_{i+1} - s_{i}\end{aligned}\end{split}\](278)\[(Ds)_{i,{\rm cc,min}} = \min\left[\left|(Ds)_{i,{\rm cc},-}\right|,\left|(Ds)_{i,{\rm cc},+}\right|\right].\]If \((Ds)_{i,{\rm face,min}} \ge (Ds)_{i,{\rm cc,min}}\), set \((Ds)_{i,\pm} = (Ds)_{i,{\rm face},\pm}\). Otherwise, set \((Ds)_{i,\pm} = (Ds)_{i,{\rm cc},\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{split}\begin{aligned} (D^2s)_{i} &=& 6(\alpha_{i,+}+\alpha_{i,-}) \\ (D^2s)_{i,L} &=& s_{i-2}-2s_{i-1}+s_{i} \\ (D^2s)_{i,R} &=& s_{i}-2s_{i+1}+s_{i+2} \\ (D^2s)_{i,C} &=& s_{i-1}-2s_{i}+s_{i+1}\end{aligned}\end{split}\]Then, define
(280)\[s = \text{sign}\left[(D^2s)_{i}\right],\](281)\[(D^2s)_{i,\text{lim}} = \max\left\{\min\left[s(D^2s)_{i},Cs(D^2s)_{i,L},Cs(D^2s)_{i,R},Cs(D^2s)_{i,C}\right],0\right\}.\]Then,
(282)\[\alpha_{i,\pm} = \frac{\alpha_{i,\pm}(D^2s)_{i,\text{lim}}}{\max\left[\left|(D^2s)_{i}\right|,1\times 10^{-10}\right]}\]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\left(\alpha_{j,+}+\alpha_{j,-}\right)},\](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\left[(\delta s)^2 - \delta s \alpha_{i,\mp}\right]^{\myhalf}\]otherwise,
(287)\[\alpha_{i,\pm} = -2\alpha_{i,\mp}\]Finally, \(s_{i,\pm} = s_i + \alpha_{i,\pm}\).