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

1. Enter advance_timestep with $$[\etarhoec, \etarhocc]^{n-\myhalf}$$.

2. Call make_w0. The spherical version uses uses $$\etarho^{{\rm ec},n-\myhalf}$$ and $$\etarho^{{\rm cc},n-\myhalf}$$.

3. Call density_advance. The plane-parallel version computes $$\etarho^{{\rm flux},n+\myhalf,*}$$.

4. Call make_etarho to compute $$[\etarhoec, \etarhocc]^{n+\myhalf,*}$$. The plane-parallel version uses $$\etarho^{{\rm flux},n+\myhalf,*}$$.

5. Call make_psi. The plane-parallel version uses $$\etarho^{{\rm cc},n+\myhalf,*}$$.

6. Call make_w0. The spherical version uses uses $$\etarho^{{\rm ec},n+\myhalf,*}$$ and $$\etarho^{{\rm cc},n+\myhalf,*}$$.

7. Call density_advance. The plane-parallel version computes $$\etarho^{{\rm flux},n+\myhalf}$$.

8. Call make_etarho to compute $$[\etarhoec, \etarhocc]^{n+\myhalf}$$. The plane-parallel version uses $$\etarho^{{\rm flux},n+\myhalf}$$.

9. 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$$

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.