We carry around three different 1D radial quantities: \(\etarhoec\) (edge-centered), \(\etarhocc\) (cell-centered), and \(\divetarho\) (cell-centered). These notes discuss when each of these is used, and how they are computed, in both plane-parallel and spherical.

The Mixing Term, \(\etarho\)

The base state evolves in response to heating and mixing in the star. The density evolution is governed by

(116)\[\frac{\partial \rho_0}{\partial t} = - \nablab \cdotb \left( \rho_0 w_0 \er \right) - \nablab \cdotb \left( \etarho \er \right) ,\]

with

(117)\[\etarho(r) = \overline{\left(\rhop \Ubt \cdot \er \right)} = \frac{1}{A(\Omega_H)} \int_{\Omega_H} (\rhop \Ubt \cdot \er ) \; dA ,\]

designed to keep the average value of the full density, \(\rho\), over a layer of constant radius in the star equal to \(\rho_0\). To complete the update of the base state, we need evolution equations for the pressure, \(p_0\), and velocity, \(w_0\). For spherical geometry, the derivation of \(w_0\) constraint equation is shown in the multilevel paper, resulting in the following system

(118)\[\begin{split}\begin{aligned} w_0 &= \ow + \delta w_0 \\ \frac{1}{r^2} \frac{\partial}{\partial r} \left (r^2 \ow \right ) &= \Sbar \\ \frac{\partial}{\partial r} \left[ \frac{\gammabar p_0}{r^2} \frac{\partial}{\partial r} (r^2 \dw) \right] &= - \frac{g}{r^2} \frac{\partial (r^2 \etarho)}{\partial r} - \frac{4 (\ow + \dw) \rho_0 g}{r} - 4 \pi G \rho_0 \etarho \label{eq:dw0constraint}\end{aligned}\end{split}\]

In paper III, we introduced a mixing term, \(\etarho\), to the density evolution equation (Eq.116), with the objective of keeping the base state density equal to the average of the density over a layer. For a spherical base state, it is best to define the average in terms of spherical coordinates,

(119)\[\overline{q} = \frac{1}{4\pi} \int_{\Omega_H} q(r,\theta,\phi) \; d\Omega\]

where \(\int_{\Omega_H} d\Omega = 4\pi\) represents the integral over the spherical \(\theta\) and \(\phi\) angles at constant radius.

Recall from Paper III that if initially \(\overline{\rho'} = 0\), there is no guarantee that \(\overline{\rho'} = 0\) will hold at later time. To see this, start with the equation for the perturbational density

(120)\[\frac{\partial\rho'}{\partial t} = -\Ub\cdot\nabla\rho' - \rho'\nabla\cdot\Ub - \nabla\cdot\left(\rho_0\Ubt\right),\]

written in a slightly different form:

(121)\[\frac{\partial\rho'}{\partial t} + \nabla\cdot\left(\rho'\Ub\right) = -\nabla\cdot\left(\rho_0\Ubt\right).\]

We integrate this over a spherical shell of thickness \(2h\) at radius \(r_0\), i.e., \(\Omega_H \times (r_0-h, r_0+h)\), and normalize by the integration volume, which is \(\sim 4\pi r_0^2 2h\) for small \(h\), to obtain:

(122)\[\begin{split}\begin{aligned} \frac{1}{4\pi r_0^2 2h}\int_{r_0-h}^{r_0+h}r^2 dr\int_{\Omega_H}\left[\frac{\partial\rho'}{\partial t} + \nabla\cdot\left(\rho'\Ub\right)\right]d\Omega =& - \frac{1}{4\pi r_0^2 2h}\int_{r_0-h}^{r_0+h}r^2 dr\int_{\Omega_H}\nabla\cdot\left(\rho_0\Ubt\right)d\Omega \nonumber \\ =& \left. -\frac{1}{4\pi r_0^2 2h}\int_{\Omega_H}\left[\rho_0\left(\Ubt\cdot\eb_r\right)\right]r^2 d\Omega\right|_{r_0-h}^{r_0+h} \nonumber \\ =& 0,\end{aligned}\end{split}\]

where we have used the divergence theorem in spherical coordinates to transform the volume integral on the right hand side into an area integral over \(\Omega_H\). We see that the right hand side disappears since \(\int_{\Omega_H}\rho_0(\Ubt\cdot\eb_r)d\Omega=0\), which follows from the definition of \(\Ubt\). Now, expanding the remaining terms and taking the limit as \(h\rightarrow 0,\) we can write

(123)\[\begin{split}\begin{aligned} 0 =& \lim_{h\rightarrow 0} \frac{1}{4\pi r_0^2 2h} \int_{r_0-h}^{r_0+h} r^2 dr \int_{\Omega_H} \left[\frac{\partial \rho'}{\partial t} + \nabla\cdot(\rho'\Ub)\right] d\Omega \nonumber \\ =& \frac{\partial}{\partial t} \left( \lim_{h\rightarrow 0} \frac{1}{4\pi r_0^2} \frac{1}{2h} \int_{r_0-h}^{r_0+h} r^2 dr \int_{\Omega_H} \rho' d\Omega \right) + \lim_{h\rightarrow 0} \left[\frac{1}{4\pi r_0^2 2h} \int_{r_0-h}^{r_0+h} r^2 dr \int_{\Omega_H} \nabla \cdot ( \rho' \Ub ) d\Omega\right] \nonumber \\ =& \frac{\partial}{\partial t} \left(\frac{1}{4\pi} \int_{\Omega_H} \rho' d\Omega \right) + \lim_{h\rightarrow 0} \left\{\left.\frac{1}{r_0^2 2h}\left[\frac{1}{4\pi}\int_{\Omega_H} \rho' (\Ub \cdot \eb_r) d\Omega \right] r^2 \right |_{r_0-h}^{r_0+h} \right\} \nonumber \\ =& \frac{\partial}{\partial t} \overline{\rho'} + \lim_{h\rightarrow 0} \left\{ \left . \frac{1}{r_0^2 2h} \left[ \overline{\rho' (\Ub \cdot \eb_r)} \right] r^2 \right |_{r_0-h}^{r_0+h} \right\} \nonumber \\ =& \frac{\partial}{\partial t} \overline{\rho'} + \lim_{h\rightarrow 0} \frac{1}{r_0^2 2h} \int_{r_0-h}^{r_0+h} \nabla \cdot \left[\overline{\rho' (\Ub \cdot \eb_r)} \eb_r \right] r^2 dr \nonumber \\ =& \frac{\partial}{\partial t} \overline{\rho'} + \nabla \cdot \left[ \overline{\rho' (\Ub \cdot \eb_r)} \eb_r \right]\end{aligned}\end{split}\]

again using the divergence theorem, extracting the time derivative from the spatial integral, and switching the order of operations as appropriate.

In short,

(124)\[\frac{\partial}{\partial t} \overline{\rho'} = - \nabla\cdot\left[\overline{\rho'\left(\Ub\cdot\eb_r\right)}\eb_r\right] = -\nabla\cdot\left(\etarho\eb_r\right),\]

and thus, \(\etarho = \overline{\left(\rho'\Ub\cdot\eb_r\right)}\).

We need both \(\etarho\) alone and its divergence for the various terms in the construction of \(w_0\) and the correction to \(\rho_0\). The quantity \(\etarho\) is edge-centered on our grid, and for Cartesian geometries, we constructed it by averaging the appropriate fluxes through the grid boundaries. For a spherical base state, this does not work, since the spherical shells do not line up with the Cartesian grid boundaries.

Therefore, we take a different approach. We compute \(\etarho\) by constructing the quantity \(\rhop \Ubt \cdot \er\) in each cell, and then use our average routine to construct a 1-d, cell-centered \(\eta_{\rho,r}\) (this is essentially numerically solving the integral in Eq.117. The edge-centered values of \(\etarho\), \(\eta_{\rho,r+1/2}\) are then constructed by simple averaging:

(125)\[\eta_{\rho,r+1/2} = \frac{\eta_{\rho,r} + \eta_{\rho,r+1}}{2} .\]

Instead of differencing \(\eta_{\rho,r+1/2}\) to construct the divergence, we instead use Eq.124 directly, by writing:

(126)\[\left [ \nabla \cdot (\etarho \er ) \right ]^{n+1/2} = - \frac{\overline{\rhop^{n+1}} - \overline{\rhop^n}}{\Delta t} = - \frac{\overline{\rhop^{n+1}}}{\Delta t} ,\]

where we have made use of the fact that \(\overline{\rhop^n} = 0\) by construction.

\(\eta\) Flow Chart

Enter advance_timestep with \([\etarhoec, \etarhocc]^{n-\myhalf}\).
Call make_w0. The spherical version uses uses \(\etarho^{{\rm ec},n-\myhalf}\) and \(\etarho^{{\rm cc},n-\myhalf}\).
Call density_advance. The plane-parallel version computes \(\etarho^{{\rm flux},n+\myhalf,*}\).
Call make_etarho to compute \([\etarhoec, \etarhocc]^{n+\myhalf,*}\). The plane-parallel version uses \(\etarho^{{\rm flux},n+\myhalf,*}\).
Call make_psi. The plane-parallel version uses \(\etarho^{{\rm cc},n+\myhalf,*}\).
Call make_w0. The spherical version uses uses \(\etarho^{{\rm ec},n+\myhalf,*}\) and \(\etarho^{{\rm cc},n+\myhalf,*}\).
Call density_advance. The plane-parallel version computes \(\etarho^{{\rm flux},n+\myhalf}\).
Call make_etarho to compute \([\etarhoec, \etarhocc]^{n+\myhalf}\). The plane-parallel version uses \(\etarho^{{\rm flux},n+\myhalf}\).
Call make_psi. The plane-parallel version uses \(\etarho^{{\rm cc},n+\myhalf}\).

Computing \(\etarhoec\) and \(\etarhocc\)

This is done in make_eta.f90.

Plane-Parallel

We first compute a radial edge-centered multifab, \(\etarho^{\rm flux}\), using

(127)\[\eta_{\rho,\ib+\myhalf\eb_r}^{\rm flux} = \left[\left(\Ubt_{\ib+\myhalf\eb_r}^{n+\myhalf}\cdot\eb_r\right) + w_{0,r+\myhalf}^{n+\myhalf}\right] \rho_{\ib+\myhalf\eb_r}^{n+\myhalf} - w_{0,r+\myhalf}^{n+\myhalf}\rho_{0,r+\myhalf}^{n+\myhalf, {\rm pred}}\]

\(\etarhoec\) is the edge-centered “average” value of \(\eta_{\rho}^{\rm flux}\),

(128)\[\eta_{\rho,r+\myhalf}^{\rm ec} = \overline{\eta_{\rho,\ib+\myhalf\eb_r}^{\rm flux}}\]

\(\etarhocc\) is a cell-centered average of \(\etarhoec\),

(129)\[\eta_{\rho,r}^{\rm cc} = \frac{\eta_{\rho,r+\myhalf}^{\rm ec} + \eta_{\rho,r-\myhalf}^{\rm ec}}{2}.\]

Spherical

First, construct \(\eta_{\rho}^{\rm cart} = [\rho'(\Ubt\cdot\eb_r)]^{n+\myhalf}\) using:

(130)\[\left[\frac{\rho^n+\rho^{n+1}}{2}-\left(\frac{\rho_0^n+\rho_0^{n+1}}{2}\right)^{\rm cart}\right] \sum_d\left(\frac{\Ubt_{\ib+\myhalf\eb_d}^{n+\myhalf}\cdot\eb_d+\Ubt_{\ib-\myhalf\eb_d}^{n+\myhalf}\cdot\eb_d}{2}\right)n_d.\]

Then, \(\etarhocc\) is the cell-centered average of \(\eta_{\rho}^{\rm cart}\),

(131)\[\etarhocc = \overline{\eta_{\rho}^{\rm cart}}.\]

On interior faces, \(\etarhoec\) is the average of \(\etarhocc\),

(132)\[\eta_{\rho,r-\myhalf}^{\rm ec} = \frac{\eta_{\rho,r-1}^{\rm cc} + \eta_{\rho,r}^{\rm cc}}{2}.\]

At the upper and lower boundaries, we use

(133)\[\begin{split}\begin{aligned} \eta_{\rho,-\myhalf}^{\rm ec} &=& 0, \\ \eta_{\rho,{\rm nr}-\myhalf}^{\rm ec} &=& \eta_{\rho,{\rm nr}-1}^{\rm cc}.\end{aligned}\end{split}\]

Using \(\etarhoec\)

Plane-Parallel

NOT USED.

Spherical

In make_w0, \(\etarhoec\) is used in the construction of the RHS for the \(\delta w_0\) equation.

Using \(\etarhocc\)

Plane-Parallel

In make_psi, \(\psi = \etarhocc g\).

Spherical

In make_w0, \(\etarhocc\) is used in the construction of the RHS for the \(\delta w_0\) equation.