MAESTRO: AN ADAPTIVE LOW MACH NUMBER HYDRODYNAMICS ALGORITHM FOR STELLAR FLOWS

For any Cartesian cell-centered field, ϕ, we define $\overline{\phi } =$ Avg(ϕ) as the lateral average over a layer at constant radius r, Ω_H, as

$$\begin{equation} \overline{\phi (r)} = \frac{1}{A(\Omega _H)}\int _{\Omega _H}\phi (r,{\bf x})d\Omega ; \qquad A(\Omega _H) = \int _{\Omega _H}d\Omega. \end{equation} \tag{ 20 }$$

Planar. This is a straightforward arithmetic average of cells at a particular height since the radial cell centers are in alignment with the Cartesian grid cell centers.

Spherical. It can be shown that any Cartesian cell center is a radius $\hat{r}_m = \Delta x\sqrt{3/4 + 2m}$ from the center of the star, where m ⩾ 0 is an integer. For example, the Cartesian cell with coordinates (i, j, k) = (1, 1, 1) relative to the center of the star lies at a distance of $\Delta x\sqrt{3/4+6}$ from the center of the star, corresponding to m = 3. The Cartesian cells with coordinates (i, j, k) = (2, 0, 0), (0, 2, 0), or (0, 0, 2) relative to the center of the star also lie at that same distance. For the 384³ resolution examples in this paper, we have verified that a non-zero set of Cartesian cell centers map into each radius $\hat{r}_m$ until m is large enough to correspond to a radius larger than half the width of the computational domain (i.e., the edge of the domain, not the corner of the domain). Figure 3 shows the number of Cartesian cells that map into each radius $\hat{r}_m$, which we refer to as the "hit count," for a 384³ domain. We use this mapping to help construct the lateral average, using the following steps.

• 1.  
  Create an itemized list, $\hat{\phi }_m$, where each element is associated with a radius $\hat{r}_m = \Delta x\sqrt{3/4 + 2m}$ from the center of the star.
• 2.  
  For each $\hat{\phi }_m$, compute the arithmetic average value of the Cartesian cells whose centers lie at the associated radius. As an additional element in the itemized list, include the center of the star (corresponding to a radius of r = 0). Compute this additional value of $\hat{\phi }$ at this location using quadratic interpolation with $\hat{\phi }_0$, $\hat{\phi }_1$, and a homogeneous Neumann condition at r = 0 as the stencil points. Note that for very large values of m, it is possible that no Cartesian cell centers exist at a radius $\hat{r}_m$ (i.e., the hit count is zero). If so, we say that $\hat{\phi }_m$ has an undefined/invalid value, and we ignore such values for the rest of this procedure.
• 3.  
  To compute the lateral average, use quadratic interpolation using the value in the itemized list with the closest associated radius, $\hat{\phi }_k$, and the nearest values above and below, $\hat{\phi }_{k_+}$ and $\hat{\phi }_{k_-}$, using divided differences:

  $$\begin{eqnarray*} \overline{\phi (r)} &=& \hat{\phi }_{k_-} + \frac{\hat{\phi }_k-\hat{\phi }_{k_-}}{\hat{r}_k-\hat{r}_{k_-}}(r-\hat{r}_{k_-})\nonumber\\ && + \frac{\frac{\hat{\phi }_{k_+}-\hat{\phi }_k}{\hat{r}_{k_+}-\hat{r}_k}-\frac{\hat{\phi }_k-\hat{\phi }_{k_-}}{\hat{r}_k-\hat{r}_{k_-}}}{\hat{r}_{k_+}-\hat{r}_{k_-}}(r-\hat{r}_{k_-})(r-\hat{r}_k), \end{eqnarray*}$$

  where $\hat{r}_{k_-}$, $\hat{r}_k$, and $\hat{r}_{k_+}$ are the three radii associated with $\hat{\phi }_{k_-}$, $\hat{\phi }_k$, and $\hat{\phi }_{k_+}$, respectively. Finally, constrain $\overline{\phi (r)}$ to lie within the range of $\hat{\phi }_{k_-}$, $\hat{\phi }_k$, and $\hat{\phi }_{k_+}$ so as to not introduce any new maxima or minima.

Zoom In Zoom Out Reset image size

In Section 6.1, we show the improvement of this averaging procedure over the Paper IV procedure.

There are four different mappings from a one-dimensional radial array to the three-dimensional Cartesian grid; below we describe the procedures for planar and spherical geometries separately.

4.2.1. Planar

4.2.2. Spherical

• 1.  
  To map a cell-centered radial array onto Cartesian cell centers, we use quadratic interpolation from the nearest three radial cell centers (see Figure 4(a)). This is a departure from Paper IV, in which we used piecewise constant interpolation.
• 2.  
  To map an edge-centered radial array onto Cartesian cell centers, we use linear interpolation from the nearest two points (see Figure 4(b)).
• 3.  
  To map a cell-centered radial array onto Cartesian edges, we first map the radial array onto Cartesian cell centers (see (1)), then average the two neighboring centers to obtain the Cartesian edge values (see Figure 4(c)).
• 4.  
  To map an edge-centered radial array onto Cartesian edges, we first map the radial array onto Cartesian cell centers (see (2)), then average the two neighboring centers to obtain the Cartesian edge values (see Figure 4(c)).

Zoom In Zoom Out Reset image size

At each time step the state data are defined on a nested hierarchy of grids, ranging from the base level (ℓ = 1), which covers the entire computational domain, to the finest level (ℓ = ℓ_max). At each level there is a union of non-intersecting rectangular grids with the same spatial resolution. For simplicity, we require that the cells composing the grids be square (Δx = Δy = Δz), and that the refinement ratio between levels be 2. The grids in the interior of the computational domain are required to be properly nested, i.e., the union of grids at level ℓ + 1 is completely contained in the union of grids at level ℓ. Additionally, in the interior, we require that each grid at level ℓ + 1 be a distance of at least two level ℓ cells from the boundary between level ℓ and level ℓ − 1 grids; this allows us to always fill "ghost cells" at level ℓ + 1 from the level ℓ data (or the physical boundary conditions, if appropriate). We initialize the grid hierarchy and regrid following the procedure outlined in Bell et al. (1994). A user-specified error estimation routine is used to tag cells where more resolution is desired. The tagged cells are grouped into rectangular patches following Berger & Rigoutsos (1991), and subsequently refined to create new grids at next level. Refinement continues until the maximum level is reached.

During Step 0,³ grids at all levels are filled directly from the initial data. As the simulation progresses, we periodically check our refinement criteria and regrid as necessary. This regridding takes places during Step 12, before computing the next time step. Newly created grids are filled by using data from previous grids at the same refinement level (if available) or by interpolating from underlying coarser grids.

Since we use the same time step to advance the solution at all levels, much of the complication associated with synchronization of data between levels in a subcycling algorithm (see Almgren et al. 1998) is eliminated. The MAC projections in Steps 3 and 7 enforce that $\widetilde{{\bf U}}^{{\rm ADV},\star }$ and $\widetilde{{\bf U}}^{\rm ADV},$ respectively, on any coarse edge underlying fine edges are the average of the values on the fine edges. Similarly, the nodal projection in Step 12 enforces that at any coarse node underlying a fine node, the value of ϕ on the coarse node is identical to the value on the fine node above it. The additional communication of data between levels occurs as follows.

• 1.  
  Before any explicit operation at level ℓ>1, data in ghost cells at that level are filled by interpolating from level ℓ − 1, or imposing physical boundary conditions, as appropriate.
• 2.  
  Edge-based fluxes at level ℓ < ℓ_max that underlie edges at level ℓ + 1 are defined to be the average of the fluxes on level ℓ + 1 at that edge. This enforces conservation.
• 3.  
  After any update to the solution, data at finer levels are conservatively averaged onto the underlying coarse grid cells, starting at the finest level.

Our specific treatment of AMR is guided by our initial scientific applications, including Type I X-ray bursts and the convective phase of SNe Ia, as well as numerical concerns, most notably the presence of the one-dimensional base state. Our treatment of the base state in an AMR framework differs for planar versus spherical problems.

For planar problems, our approach is to define a radial base state array with variable mesh spacing. A general localized fine Cartesian grid would require either a base state that exists at multiple resolutions at a particular height, or an interpolation algorithm to obtain the base state value at a particular height if Δr ≠ Δx. Both of these methods pose problems, as they generate oscillations in perturbational quantities (such as ρ', (ρh)' and the $S-\overline{S}$ term on the right-hand side of the divergence constraint) since the lateral average routine is only defined when the base state is aligned with the Cartesian grid across the width of the domain. Any attempt at interpolation will cause oscillations in the perturbational quantities directly related to the interpolation error. We have found that such oscillations can be detrimental to the results. With these issues in mind, we choose to only allow fine grids to exist that span the width of the domain. This way, the base state exists as a single seamless entity with multiple resolutions depending on height (see Figure 5). We will take advantage of this type of grid structure in our studies of Type I X-ray bursts.

Zoom In Zoom Out Reset image size

Next, we define ghost cell values for the finer base state levels, and fill these values by interpolating coarser data. This makes the algorithm directly compatible with the one-dimensional time-centered edge state calculation used in Advect Base Density⁴ and Advect Base Enthalpy. In particular, the slopes can be used with a consistent stencil at each level, that is not dependent on the data from any other level once the ghost cells are set.

Finally, whenever we regrid the Cartesian grid data, we regrid the base state to match the grid structure of the Cartesian grid. Then, we set $\rho _0=\overline{\rho }$ and compute p₀ using Enforce HSE. To compute ψ and w₀ on the new base state array, we use piecewise linear interpolation of the coarser data to fill any new fine radial cells/edges.

For spherical geometry, we first note that even in the single-level case the radial base state is not aligned with the Cartesian grid. Therefore, we use a base state with a fixed Δr for all levels (see Figure 5). As in the single-level algorithm, we choose Δr = Δx/5, but here, Δx corresponds to resolution of the Cartesian grid at the finest level.

Our next consideration is defining the radial average. First, we first create an itemized list associated with each level of refinement using only Cartesian cells that are not covered by cells at a finer level. At this point one option would be to merge the lists and proceed as in the single-level algorithm; this was tested and found to be problematic. Instead, we choose the list from a chosen particular level and define the average using quadratic interpolation with only this list, as in the single-level case. To decide which list to use, we first examine the three points that would be used by quadratic interpolation at each level. The guiding principle is to avoid using interpolation points with low hit counts. Thus, at each level we find the minimum hit count of the three points; the level which has the largest minimum hit count is the level whose list we use for interpolation. We note that this multi-level average works particularly well when the center of the star is fully refined. This is the case since near the center of the star, there are relatively few Cartesian cells that contribute to each radial bin, so by fully refining the center of the star, we ensure that the multi-level averaging procedure retains the accuracy of a single-level spherical average near the center. For our studies of SNe Ia, we will take advantage of this fact by always refining the center of the star, which is our region of interest (see Figure 12 in Section 6.6, as an example). In Section 6.1, we present numerical tests of the new multi-level averaging procedure for spherical geometry.

For both planar and spherical problems, after regridding the Cartesian grid, we make the state thermodynamically consistent by computing T = T(ρ, h, X_k) (for planar problems) or T = T(ρ, p₀, X_k) (for spherical problems). Then, we recompute $\overline{\Gamma _1}$ and β₀ as described in Steps 10 and 11.

We parallelize the algorithm by distributing the grids on each level across processors. Each grid carries a perimeter of ghost cells that are filled from neighboring grids at the same level or interpolated from coarser grids as needed. This allows the data on each grid to be updated independently of the other grids. A typical grid is large enough (e.g., 32³ cells) that for explicit operations the cost of computation within each grid greatly exceeds the cost of communication between grids. The linear solves necessary for the MAC projection and approximation nodal projection have higher communication costs, but we still obtain good parallel efficiency for the overall algorithm. A scaling study for MAESTRO can be found in Almgren et al. (2007).

Since the one-dimensional base state arrays are so much smaller than the three-dimensional arrays holding the full solution, each processor owns a copy of the entire one-dimensional base state arrays. Operations such as averaging to define base state quantities require a collection operation among grids, followed by a distribution of the average state to each processor.

To test the new spherical fill and lateral average routines from Section 4, we first create a unit cube with 384³ resolution and no refinement. We create a radial array with Δr = Δx/5, and initialize the radial array cell centers with the Gaussian profile $\phi _{\rm exact}(r) = e^{-10r^2}$. We map ϕ_exact to Cartesian cell centers with the fill operation, then compute the lateral average of this Cartesian field, Avg(ϕ). We repeat this process by choosing the grid with two levels of refinement used in the full-star simulation test in Section 6.6, shown in Figure 12.

Figure 6 (left) shows the relative error between ϕ_exact and Avg(ϕ) for the new mapping procedures and the mapping procedures from Paper IV for the single-level test. The new mapping procedure greatly decreases the relative error. Figure 6 (right) is a zoom-in of the relative error for the new mapping. For the single-level grid, the relative error is $\mathcal {O}(10^{-8})$ for r ∈ [0, 0.036] and is $\mathcal {O}(10^{-13})$ for r ∈ [0.036, 0.7]. For the test with two levels of refinement, the relative error is $\mathcal {O}(10^{-8})$ for r ∈ [0, 0.7], which is still a vast improvement when compared to the Paper IV mapping applied to a single-level simulation.

Zoom In Zoom Out Reset image size

To test the base state expansion for spherical geometry, we perform a series of tests similar to those in Paper II which tested the base state expansion in planar problems. We run the same test using three codes—a one-dimensional version of the compressible code, CASTRO (Almgren et al. 2010), in spherical coordinates; a one-dimensional version of MAESTRO in spherical coordinates, and a full three-dimensional spherical star in MAESTRO.

Our initial model is generated by specifying a core density (2.6 × 10⁹ g cm^-3), temperature (6 × 10⁸ K), and a uniform composition (X(¹²C) = 0.3, X(¹⁶O) = 0.7) and integrating the equation of hydrostatic equilibrium outward while constraining the specific entropy, s, to be constant. In discrete form, we solve

$$\begin{equation} p_{0,j+1} = p_{0,j} + \case{1}{2} \Delta r(\rho _{0,j} + \rho _{0,j+1}) g_{j+1/2}. \end{equation} \tag{ 22 }$$

$$\begin{equation} s_{0,j+1} = s_{0,j}. \end{equation} \tag{ 23 }$$

We begin with a guess of ρ_{0,j + 1} and T_{0,j + 1} and use the equation of state and Newton–Raphson iterations to find the values that satisfy our system. Since this is a spherical, self-gravitating star, the gravitation acceleration, g_{j + 1/2}, is updated each iteration based on the current value of the density. Once the temperature falls below 10⁷ K, we keep the temperature constant, and continue determining the density via hydrostatic equilibrium. This uniquely determines the initial model.

For the one-dimensional simulations, we map the inner 5 × 10⁸ cm of the model onto a one-dimensional array with 1280 elements with Δr = 3.90625 × 10⁵ cm. For the full-star three-dimensional simulation, we map the model onto a 5 × 10⁸ cm³ domain with 256³ Cartesian grid cells with Δx = 5Δr = 19.53125 × 10⁵ cm. For the one- and three-dimensional MAESTRO calculations, we use cutoff densities (see Appendix A.5) of ρ_cutoff = 10⁵ g cm^-3 and ρ_anelastic = 10⁶ g cm^-3, corresponding to radii of approximately 1.8 × 10⁸ cm and 1.9 × 10⁸ cm, so the star easily fits within the computational domain for each problem. For the full-star three-dimensional simulation, we use an inner sponge with r_sp equal to the radius where ρ₀ = 10⁷ g cm^-3, r_md equal to the radius where ρ₀ = 3 × 10⁶ g cm^-3, and κ = 10 s⁻¹. We use the same outer sponge as in Paper IV. All boundary conditions are outflow, except for the center of the one-dimensional simulations, which uses a symmetry boundary condition. We run each simulation using a CFL number of 0.5.

We heat the center of the star for 0.5 s and look at the solution at 2.0 s (after the compressible solution has had time to re-equilibrate). We use $H_{\rm ext} = H_0 e^{-(r/10^7\; {\rm cm})^2}$, with H₀ = 10¹⁶ erg g^-1 s⁻¹ chosen to be much larger than the nuclear energy generation rate during the convective phase of SN Ia, in order to see a measurable effect over a few seconds. A comparison of ρ₀, p₀, and $\overline{T}$ at t = 0 and t = 2 s for each code is shown in Figure 7. There is excellent agreement between each of the simulations.

Zoom In Zoom Out Reset image size

In Paper III, we demonstrated that our single-level algorithm is second order in space and time by tracking a hot bubble rising in a white-dwarf environment using two-dimensional planar geometry. Here, we perform the same test to show that our algorithm, with a level of refinement, is second order in space and time.

We choose a domain size of 7.2 × 10⁷ cm by 2.88 × 10⁸ cm and generate a high resolution initial model with Δr = 7.03125 × 10⁴ cm (which is equal to Δx for our highest resolution, "exact" solution) using the method described in Appendix C with r_base = 0. For each of the remaining, lower resolution simulations for this test, we generate an initial model using Δr equal to the effective Δx of each simulation by linearly interpolating values from the high resolution model. Next, we add a temperature perturbation of the form

$$\begin{equation} T_{ij} = T_{0,j} + 0.3\left[1 + \tanh \left(2-\frac{d_{ij}}{\sigma }\right)\right], \end{equation} \tag{ 24 }$$

where σ = 2.5 × 10⁶ cm and d_ij is the physical distance between the cell center corresponding to cell (i, j) and the location (3.6 × 10⁷ cm, 3.2 × 10⁷ cm). Then, we call the equation of state to compute a consistent ρ, h = ρ, h(p₀, T, X_k) everywhere. We use the reaction network described in Section 4.2 in Paper III. We use cutoff densities of ρ_cutoff = ρ_anelastic = 3 × 10⁶ g cm^-3, and a sponge with r_sp equal to the radius where ρ₀ = 10ρ_cutoff, r_md equal to the radius where ρ₀ = ρ_cutoff, and κ = 10 s⁻¹. We specify periodic boundary conditions on the side walls, outflow at the top, and a solid wall at the bottom of the domain.

Since we do not have an exact analytical solution, we consider a single-level simulation run with 1024 × 4096 cells and Δt = 3.125 × 10⁻³ s to be the exact solution. We perform three single-level simulations using resolutions of 64 × 256, 128 × 512, and 256 × 1024 grid cells using fixed time steps of Δt = 0.05 s, 0.025 s, and 0.0125 s, respectively. We also perform two simulations with a single level of refinement with effective resolutions of 128 × 512 and 256 × 1024 grid cells with fixed time steps of Δt = 0.025 s and 0.0125 s, respectively. These fixed time steps correspond to a CFL of 0.9. We refine all cells in the range r ∈ [1.8 × 10⁷, 5.4 × 10⁷] cm, so effectively we have refined 1/8th of the domain, making sure the hot spot is contained within the refined region. We run each simulation to t = 1 s. For this test, whenever we call Enforce HSE to compute p₀ from ρ₀, we use r = 1.8 × 10⁷ cm as the starting point for integration rather than the location of the cutoff density to ensure that numerical errors due to integrating the equation of hydrostatic equilibrium across simulations with different resolutions are minimized.

In order to compute the L₁ error norm for each simulation, we average the data from the exact solution onto a grid with corresponding resolution. We measure the L₁ error norm in the physical space corresponding to the refined region using

$$\begin{equation} L_1 = \frac{1}{n_{\rm cell}} \sum _{i,j} |\phi _{ij} - \phi _{ij,{\rm exact}}|, \end{equation} \tag{ 25 }$$

where n_cell is the number of cells we sum over. This form of the L₁ error norm gives us the average error per cell. We compute the convergence rate, p, between a coarser and finer simulation using

$$\begin{equation} p = \log _2\left(\frac{L_{1,{\rm coarser}}}{L_{1,{\rm finer}}}\right). \end{equation} \tag{ 26 }$$

Tables 1 and 2 show the L₁ error norms and convergence rates for the single-level and multi-level solutions, respectively. The convergence rates correspond to the two columns on either side of the reported value. We note second-order convergence in each variable. Additionally, the magnitude of the L₁ error norms for the multi-level simulations is comparable to the corresponding resolution error norms for the single-level simulations. This means that the multi-level simulations are accurately capturing the finer-scale features, as compared to the single-level simulations, i.e., the presence of coarse grid data and/or coarse-fine interfaces is not harming the solution in the refined region.

To test the ability for an adaptive, three-dimensional planar simulation to track a localized feature, we examine a hot bubble rising in a white-dwarf environment. The problem setup is exactly the same as in Section 6.3, except that we now compute in three dimensions and allow the grid structure to change with time. We choose a domain size of 7.2 × 10⁷ cm by 7.2 × 10⁷ cm 2.88 × 10⁸ cm and for each simulation, we generate an initial model with Δr = 5.625 × 10⁵ cm (which is equal to the effective Δx for both of the simulations in this test) using the method described in Appendix C with r_base = 0. We add a temperature perturbation of the form given in Equation (24), but in three dimensions with the hot spot centered at location (3.6 × 10⁷ cm, 3.6 × 10⁷ cm, 3.2 × 10⁷ cm). We will show that the adaptive simulation captures the same dynamics as the single-level simulation in a more computationally efficient manner.

We compute a single-level simulation with 128 × 128 × 512 grid cells, and an adaptive simulation with two levels of refinement at the same effective resolution. For each cell, if the $T-\overline{T} > 3\times 10^7$ K, we tag all cells at that height to ensure they are computed at the finest refinement level. We run each simulation to t = 2.5 s using a CFL number of 0.9.

Figure 8 shows the initial profile of $T-\overline{T}$ and the initial grid structure of the multi-level run. The single-level simulation has 8,388,608 grid cells and takes approximately 32 s per time step. The adaptive simulation initially has 131,072 grid cells at the coarsest level, 114,688 cells at the first level of refinement, and 688,128 cells at the finest level of refinement (the number of grid cells at the finer levels changes slightly with time as the grid structure changes to track the bubble) and takes approximately 12 s per time step, for a factor of 2.7 speedup. Both computations were performed using 32 processors on the Franklin XT4 machine at NERSC. Figure 9 shows a series of planar slices of the simulations at 0.5 s intervals in order to show that the adaptive simulation captures the same dynamics as the single-level simulation. The vertical distance shown is from z = 0 to 9.2 × 10⁷ cm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To test the expansion of the base state in a multi-level, two-dimensional planar simulation, we simulate a white-dwarf environment with an external heating layer. We also show that the multi-level simulation captures the same fine-scale structure as a single-level simulation at the same effective resolution for a short time. We performed a similar test in Paper III, but without refinement.

We choose a domain size of 2.5 × 10⁸ cm by 4 × 10⁸ cm and for each simulation, we generate an initial model with Δr equal to the effective Δx for that simulation using the method described in Appendix C with r_base = 5 × 10⁷ cm. The low entropy layer beneath our model atmosphere prevents the convective flow from reaching our lower boundary.

We apply an external heating layer of the form

$$\begin{equation} H_{{\rm ext},ij} = H_0 e^{(r_j-r_{\rm layer})^2}\left[1+\sum _{m=1}^3 b_m\sin \left(\frac{k_m\pi x_i}{L_x} + \Psi _m\right)\right] \end{equation} \tag{ 27 }$$

with r_layer = 1.25 × 10⁸ cm, H₀ = 2.5 × 10¹⁶ erg g⁻¹ s⁻¹, L_x = 2.5 × 10⁸ cm, and r_j and x_i being the radial and horizontal physical coordinates of cell (i, j). The perturbation parameters are b = (0.00625, 0.01875, 0.0125), k = (2, 6, 8), and Ψ = (0, π/3, π/5). We disable reactions for this test, since the heating layer was chosen to expand the base state over a very short time period, rather than accurately model reactions.

We use cutoff densities of ρ_cutoff = ρ_anelastic = 3 × 10⁶ g cm^-3 and a sponge with r_sp equal to the radius where ρ₀ = 10⁸ g cm^-3, r_md equal to the radius where ρ₀ = ρ_cutoff, and κ = 100 s⁻¹. We specify periodic boundary conditions on the side walls, outflow at the top, and a solid wall at the bottom of the domain.

In this test, we use a CFL number of 0.9. We perform two single-level simulations using 80 × 128 and 320 × 512 cells, and a simulation with two levels of refinement and an effective resolution of 320 × 512 cells. For the multi-level simulation, we fix the refined grids, ensuring that r ∈ [9.375 × 10⁷, 1.5626 × 10⁸] cm (which contains the external heating layer) is at the finest level of refinement. We run each simulation to t = 30 s.

Figure 10 shows ρ₀, p₀, and $\overline{T}$ after 30 s of convection for the single-level fine grid simulation and the multi-level simulation. There is an excellent agreement between these two simulations, except in the temperature profiles at the top of the domain. However, this corresponds to a region with density below the cutoff densities where the temperature is extremely sensitive to small density perturbations, and furthermore, is not fully refined in this test, so this is an acceptable difference. Both simulations were performed using four processors; the single-level fine grid run required approximately 1.9 s per time step and the multi-level run required approximately 0.6 s per time step, for a factor of 3 speedup.

Zoom In Zoom Out Reset image size

Figure 11 shows the temperature profile after 3 and 4 s for each of the three simulations. The vertical distance shown is from z = 5 × 10⁷ cm to 2.2 × 10⁸ cm. At t = 3 s, the multi-level simulation is able to capture the finer-scale structure visible in the single-level fine grid simulation, which is not captured in the single-level coarse grid simulation. At t = 4 s, a finer-scale structure is still visible in the multi-level simulation, but the solution begins to drift from the single-level fine grid simulation, which is expected since we are deliberately refining only a part of the convective region.

Zoom In Zoom Out Reset image size

We now compute three-dimensional, full-star calculations of convection in a white dwarf. This problem models the convection and energetics of a white dwarf that is a few hours from reaching ignition. We performed similar simulations using an earlier version of the algorithm in Paper IV.

We begin with the one-dimensional white-dwarf model described in Section 2.4 in Paper IV. We map this one-dimensional model into the center of a computational domain of 5 × 10⁸ cm on a side. The first simulation is a single-level simulation with 384³ cells. The second simulation is adaptive with two levels of refinement and an effective resolution of 384³ cells. We use the reaction network strategy from Chamulak et al. (2008) to compute the energetics from the carbon burning. This modified network differs from the one used in Paper IV in the ash composition (we now burn to an ash with A = 18 and Z = 8.8) and the energy release (we use a quadratic fit to the q-values tabulated on page 164 of Chamulak et al. 2008). Finally, instead of destroying two carbon nuclei for each reaction, we use the M₁₂ value of 2.93 described in that paper. We initialize the simulation with a velocity perturbation described exactly as in Section 2.4 in Paper IV. We use the same cutoff densities, sponge parameters, and boundary conditions as in Section 6.2.

Figure 12 shows the initial grid structure of the adaptive simulation. Based on our work in Paper IV, we choose to fully refine all cells where ρ>10⁸ g cm^-3, since at early times, the dynamics of the star are driven by the reactions and convection in this inner region. We wish to examine whether the adaptive simulation can give the same result as the single-level simulation, and the computational efficiency of each simulation.

Zoom In Zoom Out Reset image size

We use a CFL number of 0.7 and compute to t = 900 s. We choose two diagnostics used in Paper IV to compare the simulations. Peak temperature is a useful diagnostic since the reaction rates are extremely sensitive to temperature, and thus peak temperature serves as a good guide for observing the progression toward ignition. Peak radial velocity is another useful diagnostic as it is a simple measure of the strength of the convection within the star. Since the solution of our equation is highly nonlinear, and the reaction rates scale with temperature as ∼T²³ (Woosley et al. 2004), we expect that errors from the coarse grid will perturb the solution at the finest level, eventually causing significant differences in the exact flow field. However, as shown in Paper IV, when we run our simulation to the point of ignition, we require upward of hundreds of thousands of time steps. Therefore, in the comparison diagnostics in this test, it is sufficient to compare the overall qualitative solution. An exact quantitative comparison is impossible over long times.

Figure 13 shows the evolution of the peak temperature and peak radial velocity over the first 900 s for both simulations. The adaptive simulation gives the same qualitative result as the single-level simulation. As mentioned before, the curves do not match up more closely because the equations are highly nonlinear, and slight differences in the solution caused by solver tolerance and discretization error change details of the results, but not the overall qualitative results. The single-level simulation has 56,623,760 grid cells and takes approximately 36 s per time step. The adaptive simulation initially has 884,736 grid cells at the coarsest level, 3,511,808 cells at the first level of refinement, and 4,282,048 cells at the finest level of refinement (the number of grid cells at the finer levels changes slightly with time) and takes approximately 18 s per time step, for a factor of 2 speedup. Also note that the overall memory requirements are significantly less for the adaptive simulation, as can be seen by the reduced total cell count. Each simulation was run using the Franklin XT4 machine at NERSC with 216 processors.

Zoom In Zoom Out Reset image size

We note that in the future, when we perform longer calculations up to the point of ignition, we may have to refine a greater portion of the star in order to properly capture the overall dynamics as the convective region expands. In a forthcoming paper, we plan on performing higher resolution studies, where AMR will save us both time and computational resources.

The single-level algorithm has gone through numerous changes since Papers III and IV. The current single-level algorithm is presented in Sections A.2–A.5; we summarize here the changes since Papers III and IV:

• 1.  
  We have extended the base state evolution to spherical problems by defining Advect Base Density Spherical (Section A.3.2), Advect Base Enthalpy Spherical (Section A.3.4), and Compute w₀ Spherical (Section A.3.5). We have also defined spherical discretizations for ψ and η_ρ (Steps 4C, 4F, 8C, and 8F in Section A.4).
• 2.  
  For spherical problems, we have improved the accuracy of the mapping of data between the one-dimensional radial array and the three-dimensional Cartesian grid (Section 4).
• 3.  
  We have upgraded our unsplit piecewise-linear Godunov method for time-centered edge state prediction for the base state and Cartesian grid data to use the unsplit piecewise-parabolic method (PPM; Colella & Woodward 1984) with full corner coupling in three dimensions (Miller & Colella 2002; Saltzman 1994). We shall refer to this procedure as "computing the time-centered edge states" (Section A.3.2, Section A.3.4, and Steps 3, 4, 7, 8, and 11 in Section A.4).
• 4.  
  As introduced in Paper IV, we update T after the call to React State (Section A.3.1).
• 5.  
  Rather than evolve p₀ using an evolution equation, we simply update p₀ using the hydrostatic equilibrium equation (Section A.3.3). These two methods are analytically equivalent, but in our experience, the numerical drift associated with evolving p₀ using an evolution equation causes the entire solution to drift from thermodynamic equilibrium over time.
• 6.  
  As explained in Section 2, in the advection step, we predict (ρh)' to edges instead of T (Steps 4Hi and 8Hi in Section A.4). Thus, a base state enthalpy, (ρh)₀ is required in order to construct the enthalpy fluxes for the conservative update. Unlike ρ₀, we do not need to carry the complete evolution of (ρh)₀. In practice, we set (ρh)₀ equal to the lateral average of the full state enthalpy after the first call to React State, i.e., $(\rho h)_0^n = \overline{(\rho h)^{(1)}}$. We then advect (ρh)₀ without reactions or heating to mirror the treatment of (ρh) in the advection step (Steps 4G and 8G in Section A.4).
• 7.  
  In the advection step, rather than simultaneously updating each of the base state quantities, followed by an update of all the full state quantities, we have interwoven these updates in order to obtain better accuracy. In order: we advect ρ₀, advect ρ, correct ρ₀, advance p₀, advect (ρh)₀, and advect (ρh) (Steps 4 and 8 in Section A.4). This enables us to use a time-centered ψ rather than a time-lagged ψ for the enthalpy update.
• 8.  
  Rather than compute ∇ · η_ρ and use it to adjust ρ₀ to account for convecting overturning, as was done in Correct Base in Paper III, we simply use the lateral average operator to enforce $\rho _0 = \overline{\rho }$, which is simpler and analytically equivalent (Steps 4D and 8D in Section A.4). We use the improved averaging procedure described in Section 4.1, which greatly improves the accuracy of this mapping for spherical problems.
• 9.  
  We have moved the first reaction step (formerly Step 3 in Section A.4) to occur before Steps 1 and 2 in Section A.4. This is a non-functional change; it was only made so that Steps 2–5 mirror Steps 6–9 in the overall predictor–corrector scheme.
• 10.  
  For planar problems, we evolve the enthalpy for the sole purpose of updating the temperature, which is subsequently fed into the reaction network and used in computing thermodynamic derivatives. For spherical problems, as described in Paper IV, we instead define T from ρ, p₀, and X_k. This completely decouples the enthalpy from the rest of the solution. Our experience has shown that, with spherical geometry, the discretization errors are minimized by using the hydrostatic, radial base state pressure to define temperature. We still evolve the enthalpy, but we use it as a diagnostic to examine to what extent our numerical method keeps the solution in thermodynamic equilibrium.

For planar problems, we treat gravity as constant in space and time. For spherical problems, gravity is computed directly from ρ₀ in the following manner. First, we define the enclosed mass at radial edges and cell centers:

$$\begin{eqnarray} &&m_{{\rm encl},j-1/2} = \sum _{k=1}^j \frac{4}{3}\pi \big(r_{k-1/2}^3 - r_{k-3/2}^3\big)\rho _{0,k-1}; \qquad m_{{\rm encl},j} = m_{{\rm encl},j-1/2} + \frac{4}{3}\pi \big(r_j^3 - r_{j-1/2}^3\big)\rho _{0,j}. \end{eqnarray} \tag{ A1 }$$

In practice, we compute (r³_{k − 1/2} − r³_{k − 3/2}) as Δr(r²_{k − 1/2} + r_{k − 1/2}r_{k − 3/2} + r²_{k − 3/2}), in order to minimize roundoff error at large radii and use a similar formula for the radial difference term in the equation for m_encl,j. Next, we define the gravity at both radial edges and cell centers:

$$\begin{equation} g_{j-1/2} = \frac{Gm_{{\rm encl},j-1/2}}{r_{j-1/2}^2}; \qquad g_j = \frac{Gm_{{\rm encl},j}}{r_j^2}. \end{equation} \tag{ A2 }$$

We indicate which base state density is used to compute gravity by using a shorthand notation with superscripts on g, e.g., gⁿ ≡ g(ρⁿ₀).

We make use of the following shorthand notations in outlining the algorithm:

A.3.1. Reactions

React State $[\rho ^{\rm in},(\rho h)^{\rm in},X_k^{\rm in},T^{\rm in}, (\rho H_{\rm ext})^{\rm in}, p_0^{\rm in}] \rightarrow [\rho ^{\rm out}, (\rho h)^{\rm out}, X_k^{\rm out}, T^{\rm out}, (\rho \dot{\omega }_k)^{\rm out}, (\rho H_{\rm nuc})^{\rm out}]$ evolves the species and enthalpy due to reactions through Δt/2 according to

$$\begin{equation} \frac{dX_k}{dt} = \dot{\omega }_k(\rho,X_k,T) ; \qquad \frac{dT}{dt} = \frac{1}{c_p} \left(-\!\sum _k \xi _k \dot{\omega }_k + H_{\rm nuc}\right). \end{equation} \tag{ A3 }$$

Full details of the solution procedure can be found in Paper III. We then define

$$\begin{eqnarray} (\rho \dot{\omega }_k)^{\rm out} &=& \frac{\rho ^{\rm out}\big(X_k^{\rm out} - X_k^{\rm in}\big)}{\Delta t/2}, \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} (\rho h)^{\rm out} &=& (\rho h)^{\rm in} + \frac{\Delta t}{2} (\rho H_{\rm nuc})^{\rm out} + \frac{\Delta t}{2} (\rho H_{\rm ext})^{\rm in}, \end{eqnarray} \tag{ A5 }$$

where the enthalpy update includes external heat sources (ρH_ext)ⁱⁿ. As introduced in Paper IV, we update the temperature using T^out = T(ρ^out, h^out, X^out_k) for planar geometry or T^out = T(ρ^out, pⁱⁿ₀, X^out_k) for spherical geometry. Note that the density remains unchanged within React State, i.e., ρ^out = ρⁱⁿ.

A.3.2. ρ₀ Advection

Advect Base Density [ρⁱⁿ₀, wⁱⁿ₀] → [ρ^out₀, ρ^{out,n + 1/2}₀] is the process by which we update the base state density through Δt in time. We keep the time-centered edge states, ρ^{out,n + 1/2}₀, since they are used later in discretization of η_ρ for planar problems.

Planar. We discretize Equation (14) to compute the new base state density,

$$\begin{eqnarray} \rho _{0,j}^{\rm out} &=& \rho _{0,j}^{\rm in} - \frac{\Delta t}{\Delta r} \big[ \big(\rho _0^{{\rm out},{n+1/2}} w_0^{\rm in}\big)_{j+1/2} - \big(\rho _0^{{\rm out},{n+1/2}} w_0^{\rm in}\big)_{j-1/2} \big]. \end{eqnarray} \tag{ A6 }$$

We compute the time-centered edge states, ρ^{out,n + 1/2}₀, by discretizing an expanded form of Equation (14):

$$\begin{equation} \frac{\partial \rho _0}{\partial t} + w_0 \frac{\partial \rho _0}{\partial r} = - \rho _0 \frac{\partial w_0}{\partial r}, \end{equation} \tag{ A7 }$$

where the right-hand side is used as the force term.

Spherical. The base state density update now includes the area factors in the divergences:

$$\begin{eqnarray} \rho _{0,j}^{\rm out} &=& \rho _{0,j}^{\rm in} - \frac{1}{r_j^2} \frac{\Delta t}{\Delta r} \big[ \big(r^2 \rho _0^{{\rm out},{n+1/2}} w_0^{\rm in}\big)_{j+1/2} - \big(r^2 \rho _0^{{\rm out},{n+1/2}} w_0^{\rm in}\big)_{j-1/2} \big]. \end{eqnarray} \tag{ A8 }$$

In order to compute the time-centered edge states, an additional geometric term is added to the forcing, due to the spherical discretization of (14):

$$\begin{equation} \frac{\partial \rho _0}{\partial t} + w_0 \frac{\partial \rho _0}{\partial r} = - \rho _0 \frac{\partial w_0}{\partial r} - \frac{2 \rho _0 w_0}{r}. \end{equation} \tag{ A9 }$$

A.3.3. p₀ Update

Enforce HSE [pⁱⁿ₀, ρⁱⁿ₀] → [p^out₀] has replaced Advect Base Pressure from Paper III as the process by which we update the base state pressure. Rather than discretizing the evolution equation for p₀, we enforce hydrostatic equilibrium directly, which is numerically simpler and analytically equivalent. We first set p^out_{0,j = 0} = pⁱⁿ_{0,j = 0} and then update p^out₀ using

$$\begin{equation} p_{0,j+1}^{\rm out} = p_{0,j}^{\rm out} + \Delta r g_{j+1/2}\frac{\big(\rho _{0,j+1}^{\rm in}+\rho _{0,j}^{\rm in}\big)}{2}, \end{equation} \tag{ A10 }$$

where g = g(ρⁱⁿ₀). As soon as ρⁱⁿ_0,j < ρ_cutoff, we set p^out_{0,j + 1} = p^out_0,j for all remaining values of j. Then we compare $p_{0,j_{\rm max}}^{\rm out}$ with $p_{0,j_{\rm max}}^{\rm in}$ and offset every element in p^out₀ so that $p_{0,j_{\rm max}}^{\rm out}= p_{0,j_{\rm max}}^{\rm in}$. We are effectively using the location where the ρⁱⁿ₀ drops below ρ_cutoff as the starting point for integration.

A.3.4. (ρh)₀ Advection

Advect Base Enthalpy [(ρh)ⁱⁿ₀, wⁱⁿ₀, ψⁱⁿ] → [(ρh)^out₀] is the process by which we update the base state enthalpy through Δt in time

Planar. We discretize Equation (15), neglecting reaction source terms, to compute the new base state enthalpy,

$$\begin{eqnarray} (\rho h)_{0,j}^{\rm out} &=& (\rho h)_{0,j}^{\rm in} - \frac{\Delta t}{\Delta r} \big\lbrace \big[ (\rho h)_0^{{n+1/2}} w_0^{\rm in}\big]_{j+1/2} - \big[ (\rho h)_0^{{n+1/2}} w_0^{\rm in}\big]_{j-1/2} \big\rbrace + \Delta t\psi _j^{\rm in}. \end{eqnarray} \tag{ A11 }$$

We compute the time-centered edge states, (ρh)^{n + 1/2}₀, by discretizing an expanded form of Equation (15):

$$\begin{equation} \frac{\partial (\rho h)_0}{\partial t} + w_0 \frac{\partial (\rho h)_0}{\partial r} = -(\rho h)_0 \frac{\partial w_0}{\partial r} + \psi. \end{equation} \tag{ A12 }$$

Spherical. The base state enthalpy update now includes the area factors in the divergences:

$$\begin{eqnarray} (\rho h)_{0,j}^{\rm out} &=& (\rho h)_{0,j}^{\rm in} - \frac{1}{r_j^2} \frac{\Delta t}{\Delta r} \big\lbrace \big[ r^2 (\rho h)_0^{{n+1/2}} w_0^{\rm in}\big]_{j+1/2} - \big[ r^2 (\rho h)_0^{{n+1/2}} w_0^{\rm in}\big]_{j-1/2} \big\rbrace +\Delta t\psi ^{{\rm in},{n+1/2}}. \end{eqnarray} \tag{ A13 }$$

In order to compute the time-centered edge states, an additional geometric term is added to the forcing, due to the spherical discretization of (15):

$$\begin{equation} \frac{\partial (\rho h)_0}{\partial t} + w_0 \frac{\partial (\rho h)_0}{\partial r} = -(\rho h)_0 \frac{\partial w_0}{\partial r} - \frac{2 (\rho h)_0 w_0}{r} + \psi. \end{equation} \tag{ A14 }$$

A.3.5. Computing w₀

Here we describe the process by which we compute w₀. The arguments are different for planar and spherical geometries.

Compute w₀ Planar$[\overline{S}^{\rm in},\overline{\Gamma _1}^{\rm in}, p_0^{\rm in},\psi ^{\rm in}]\rightarrow [w_0^{\rm out}]$:

In Paper III, we showed that ψ = η_ρg for planar geometries, and derived an alternate expression for Equation (11). We compute w₀ using Equation (35) in Paper III using the following discretization:

$$\begin{equation} \frac{w_{0,j+1/2}^{\rm out}-w_{0,j-1/2}^{\rm out}}{\Delta r} = \left(\overline{S}^{\rm in} - \frac{1}{\overline{\Gamma _1}^{\rm in} p_0^{\rm in}}\psi ^{\rm in}\right)_j, \end{equation} \tag{ A15 }$$

with w_0,−1/2 = 0.

Compute w₀ Spherical $[\overline{S}^{\rm in},\overline{\Gamma _1}^{\rm in},\rho _0^{\rm in},p_0^{\rm in},\eta _\rho ^{\rm in}] \rightarrow [w_0^{\rm out}]$:

We begin with Equation (11) written in spherical coordinates:

$$\begin{equation} \frac{1}{r^2}\frac{\partial }{\partial r} \left(r^2 \beta _0 w_0 \right) = \beta _0 \left(\overline{S}- \frac{1}{\overline{\Gamma _1}p_0} \frac{\partial p_0}{\partial t} \right). \end{equation} \tag{ A16 }$$

We expand the spatial derivative and recall from Paper I that

$$\begin{equation} \frac{1}{\overline{\Gamma _1}p_0} \frac{\partial p_0}{\partial r} = \frac{1}{\beta _0} \frac{\partial \beta _0}{\partial r}, \end{equation} \tag{ A17 }$$

giving

$$\begin{equation} \frac{1}{r^2} \frac{\partial }{\partial r} (r^2 w_0) = \overline{S}- \frac{1}{\overline{\Gamma _1}p_0} \underbrace{\left(\frac{\partial p_0}{\partial t} + w_0 \frac{\partial p_0}{\partial r} \right)}_{\psi }. \end{equation} \tag{ A18 }$$

We solve this equation for w₀ as described in Appendix B.

We now describe the main algorithm, making frequent use of the shorthand developed above. In summary, in the predictor step (Steps 2–5) we use an estimate of the expansion term, S, to compute a preliminary solution at the new time level, denoted with an "n + 1, ⋆" superscript. In the corrector step (Steps 6–9), we use the results from the predictor step to compute a more accurate expansion term, and compute the final solution at the new time level, denoted with an "n + 1" superscript. We use Strang splitting to achieve second-order accuracy in time. See Figure 14 for a flow chart of the algorithm, including the notation used as we advance the solution by Δt. Figure 15 is a flow chart of the advection steps (Steps 4 and 8), which includes the notation we use as we advect the solution through a time interval of Δt.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The discussion that follows mirrors closely that in Paper III, but has been updated to reflect all the changes throughout the algorithm. The advance of the state through a single time step appears as follows:

Step 0. Initialization

This step remains unchanged from Paper III. The initialization step only occurs at the beginning of the simulation. The initial values for U⁰, ρ⁰, (ρh)⁰, X⁰_k, T⁰, ρ⁰₀, p⁰₀, and $\overline{\Gamma _1^0}$ are specified from the problem-dependent initial conditions. The initial time step, Δt⁰, is computed as in Paper III. Finally, initial values for w^−1/2₀, η^−1/2_ρ, ψ^−1/2, π^−1/2, S⁰, and S¹ come from a preliminary pass through the algorithm.

Step 1. React the full state through the first time interval of Δt/2.

Call React State$[\rho ^n, (\rho h)^n, X_k^n, T^n, (\rho H_{\rm ext})^n, p_0^n] \rightarrow [\rho ^{(1)},(\rho h)^{(1)},X_k^{(1)},T^{(1)},(\rho \dot{\omega }_k)^{(1)},(\rho H_{\rm nuc})^{(1)}]$.

Step 2. Compute the provisional time-centered expansion, S^{n + 1/2,⋆⋆}, provisional base state velocity, w^{n + 1/2,⋆}₀, and provisional base state velocity forcing.

• A.  
  Compute S^{n + 1/2,⋆⋆}. We compute an estimate for the time-centered expansion term in the velocity divergence constraint (Equation (12)). For the first time step (n = 0), we set

  $$\begin{equation} S^{n+1/2,\star \star } = \frac{S^0 + S^1}{2}, \end{equation} \tag{ A19 }$$

  where S¹ is found during initialization. For other time steps (n ≠ 0), following Bell et al. (2004), we extrapolate to the half-time using S at the previous and current time levels

  $$\begin{equation} S^{{n+1/2},\star \star } = S^n + \frac{\Delta t^n}{2} \frac{S^n - S^{n-1}}{\Delta t^{n-1}}. \end{equation} \tag{ A20 }$$

  Next, compute

  $$\begin{equation} \overline{S^{{n+1/2},\star \star }} = {\rm \bf Avg} (S^{{n+1/2},\star \star }). \end{equation} \tag{ A21 }$$
• B.  
  Compute w^{n + 1/2,⋆}₀.For planar geometry, callCompute w₀ Planar$[\overline{S^{{n+1/2},\star \star }},\overline{\Gamma _1^n},p_0^n,\psi ^{n-1/2}] \rightarrow [w_0^{{n+1/2},\star }]$.For spherical geometry, callCompute w₀ Spherical$[\overline{S^{{n+1/2},\star \star }},\overline{\Gamma _1^n},\rho _0^n,p_0^n,\eta _\rho ^{n-1/2}] \rightarrow [w_0^{{n+1/2},\star }]$.
• C.  
  Compute the provisional base state velocity forcing. Rearrange Equation (9),

  $$\begin{equation} -\frac{1}{\rho _0} \frac{\partial \pi _0}{\partial r} = \frac{\partial w_0}{\partial t} + w_0 \frac{\partial w_0}{\partial r}, \end{equation} \tag{ A22 }$$

  with the following discretization:

  $$\begin{equation} \left(\frac{1}{\rho _0} \frac{\partial \pi _0}{\partial r} \right)^{n,\star } = -\frac{w_0^{{n+1/2},\star } - w_0^{n-1/2}}{(\Delta t^n+\Delta t^{n-1})/2} - w_0^{n,\star } \left(\frac{\partial w_0}{\partial r}\right)^{n,\star }, \end{equation} \tag{ A23 }$$

  where w^n,⋆₀ and (∂w₀/∂r)^n,⋆ are defined as

  $$\begin{eqnarray} w_0^{n,\star } &=& \frac{\Delta t^{n} w_0^{{n-1/2}} + \Delta t^{n-1} w_0^{{n+1/2},\star }}{\Delta t^n+\Delta t^{n-1}}, \end{eqnarray} \tag{ A24 }$$

  $$\begin{eqnarray} \left(\frac{\partial w_0}{\partial r}\right)^{n,\star } &=& \frac{1}{\Delta t^n+\Delta t^{n-1}}\left[ \Delta t^{n} \left(\frac{\partial w_0 }{ \partial r}\right)^{{n-1/2}} + \Delta t^{n-1} \left(\frac{\partial w_0 }{ \partial r}\right)^{{n+1/2},\star } \right]. \end{eqnarray} \tag{ A25 }$$

  If n = 0, we use Δt⁻¹ = Δt⁰.

Step 3. Construct the provisional time-centered advective velocity on edges, $\widetilde{{\bf U}}^{{\rm ADV},\star }$.

Using Equation (10), we compute time-centered edge velocities, $\widetilde{{\bf U}}^{{\rm ADV},\dagger,\star }$, using ${\bf U}= \widetilde{{\bf U}}^n + w_0^{{n+1/2},\star }$. The † superscript refers to the fact that the predicted velocity field does not satisfy the divergence constraint. We then construct $\widetilde{{\bf U}}^{{\rm ADV},\star }$ from $\widetilde{{\bf U}}^{{\rm ADV},\dagger,\star }$ using a MAC projection, as described in detail in Appendix B of Paper III. We note that $\widetilde{{\bf U}}^{{\rm ADV},\star }$ satisfies the discrete version of $\overline{(\widetilde{{\bf U}}^{{\rm ADV},\star }\cdot {\bf e}_r)}=0$ as well as

$$\begin{equation} \nabla \cdot \big(\beta _0^n \widetilde{{\bf U}}^{{\rm ADV},\star }\big) = \beta _0^n \big(S^{{n+1/2},\star \star } - \overline{S^{{n+1/2},\star \star }}\big), \end{equation} \tag{ A26 }$$

$$\begin{equation} \beta _0^n = \beta _0 \big(\rho _0^n, p_0^n, \overline{\Gamma _1^n}\big), \end{equation} \tag{ A27 }$$

where β₀ is computed as described in Appendix C of Paper III.

Step. 4. Advect the base state and full state through a time interval of Δt.

• A.  
  Update ρ₀, saving the time-centered density at radial edges by callingAdvect Base Density [ρⁿ₀, w^{n + 1/2,⋆}₀] → [ρ^(2a),⋆₀, ρ^{n + 1/2,⋆,pred}₀].
• B.  
  Update (ρX_k) using a discretized version of Equation (1) omitting the reaction terms, which were already accounted for in React State. The update consists of two steps:

  • i.  
    Compute the time-centered species edge states, (ρX_k)^{n + 1/2,⋆,pred}, for the conservative update of (ρX_k)⁽¹⁾. We use Equations (17) and (13) to predict $\rho ^{^{\prime }(1)} = \rho ^{(1)} - \rho _0^n$ and X⁽¹⁾_k = (ρX_k)⁽¹⁾/ρ⁽¹⁾ to time-centered edges using ${\bf U}= \widetilde{{\bf U}}^{{\rm ADV},\star }+w_0^{{n+1/2},\star }{\bf e}_r$, yielding $\rho ^{^{\prime }{n+1/2},\star,{\rm pred}}$ and X^{n + 1/2,⋆,pred}_k. We convert the perturbational density full state density using

    $$\begin{equation} \rho ^{{n+1/2},\star,{\rm pred}} = \rho ^{^{\prime }{n+1/2},\star,{\rm pred}} + \frac{\rho _0^n + \rho _0^{(2a),\star }}{2}, \end{equation} \tag{ A28 }$$

    where the base state density terms are mapped to Cartesian edges. Then,(ρX_k)^{n + 1/2,⋆,pred} = (ρ^{n + 1/2,⋆,pred}X^{n + 1/2,⋆,pred}_k).
  • ii.  
    Evolve (ρX_k)⁽¹⁾ → (ρX_k)^(2),⋆ using

    $$\begin{eqnarray} (\rho X_k)^{(2),\star } &=& (\rho X_k)^{(1)} - \Delta t\big\lbrace \nabla \cdot \big[ \big(\widetilde{{\bf U}}^{{\rm ADV},\star }+w_0^{{n+1/2},\star } {\bf e}_r\big) (\rho X_k)^{{n+1/2},\star,{\rm pred}} \big] \big\rbrace, \end{eqnarray} \tag{ A29 }$$

    $$\begin{equation} \rho ^{(2),\star } = \sum _k (\rho X_k)^{(2),\star }, \qquad X_k^{(2),\star } = (\rho X_k)^{(2),\star } / \rho ^{(2),\star }. \end{equation} \tag{ A30 }$$
• C.  
  Define a radial edge-centered η^{n + 1/2,⋆}_ρ.For planar geometry, since $\eta _\rho = \overline{\rho ^{\prime }({\bf U}\cdot {\bf e}_r)} = \overline{\rho ({\bf U}\cdot {\bf e}_r)}-\overline{\rho _0({\bf U}\cdot {\bf e}_r}) = \overline{\rho ({\bf U}\cdot {\bf e}_r)} - \rho _0w_0$,

  $$\begin{eqnarray} \eta _\rho ^{{n+1/2},\star } &=& {\rm {\bf Avg}} \sum _k \big[ \big(\widetilde{{\bf U}}^{{\rm ADV},\star }\cdot {\bf e}_r + w_0^{{n+1/2},\star }\big) (\rho X_k)^{{n+1/2},\star,{\rm pred}} \big]- w_0^{{n+1/2},\star } \rho _0^{{n+1/2},\star,{\rm pred}}. \end{eqnarray} \tag{ A31 }$$

  For spherical geometry, first construct η^{cart,n + 1/2,⋆}_ρ = [ρ'(U · e_r)]^{n + 1/2,⋆} on Cartesian cell centers using

  $$\begin{eqnarray} \eta _\rho ^{{\rm cart},{n+1/2},\star } &=& \left[\left(\frac{\rho ^{(1)}+\rho ^{(2),\star }}{2}\right)-\left(\frac{\rho _0^n+\rho _0^{(2a),\star }}{2}\right)\right] \cdot \big(\widetilde{{\bf U}}^{{\rm ADV},\star }\cdot {\bf e}_r + w_0^{{n+1/2},\star }\big). \end{eqnarray} \tag{ A32 }$$

  Then,

  $$\begin{equation} \eta _\rho ^{{n+1/2},\star } = {\rm {\bf Avg}}\big(\eta _\rho ^{{\rm cart},{n+1/2},\star }\big). \end{equation} \tag{ A33 }$$

  This gives a radial cell-centered η^{n + 1/2,⋆}_ρ. To get η^{n + 1/2,⋆}_ρ at radial edges, average the two neighboring radial cell-centered values.
• D.  
  Correct ρ₀ by setting ρ^{n + 1,⋆}₀ = Avg(ρ^(2),⋆).
• E.  
  Update p₀ by calling Enforce HSE[pⁿ₀, ρ^{n + 1,⋆}₀] → [p^{n + 1,⋆}₀].
• F.  
  Compute ψ^{n + 1/2,⋆}.For planar geometry,

  $$\begin{equation} \psi _j^{{n+1/2},\star } = \case{1}{2} \big(\eta _{\rho,j-1/2}^{{n+1/2},\star } + \eta _{\rho,j+1/2}^{{n+1/2},\star }\big) g. \end{equation} \tag{ A34 }$$

  For spherical geometry, first compute

  $$\begin{eqnarray} \overline{\Gamma _1^{(1)}} &=& {\rm {\bf Avg}} \big[ \Gamma _1\big(\rho ^{(1)}, p_0^{n}, X_k^{(1)}\big) \big], \end{eqnarray} \tag{ A35 }$$

  $$\begin{eqnarray} \overline{\Gamma _1^{(2),\star }} &=& {\rm {\bf Avg}} \big[ \Gamma _1\big(\rho ^{(2),\star }, p_0^{n+1,\star }, X_k^{(2),\star }\big) \big]. \end{eqnarray} \tag{ A36 }$$

  Then, define ψ^{n + 1/2,⋆} using Equation (A18)

  $$\begin{eqnarray} \psi _j^{{n+1/2},\star } &=& \left(\frac{\overline{\Gamma _1^{(1)}}+\overline{\Gamma _1^{(2),\star }}}{2}\right)_j \left(\frac{p_0^n+p_0^{n+1,\star }}{2}\right)_j \bigg\lbrace \overline{S_j^{{n+1/2},\star }} - \frac{1}{r_j^2} \big[ \big(r^2 w_0^{{n+1/2},\star }\big)_{j+1/2} - \big(r^2 w_0^{{n+1/2},\star }\big)_{j-1/2} \big] \bigg\rbrace. \end{eqnarray} \tag{ A37 }$$
• G.  
  Update (ρh)₀. First, compute (ρh)ⁿ₀ = Avg[(ρh)⁽¹⁾]. Then, callAdvect Base Enthalpy [(ρh)ⁿ₀, w^{n + 1/2,⋆}₀, ψ^{n + 1/2,⋆}] → [(ρh)^{n + 1,⋆}₀].
• H.  
  Update the enthalpy using a discretized version of Equation (3), again omitting the reaction and heating terms since we already accounted for them in React State. This equation takes the form

  $$\begin{equation} \frac{\partial (\rho h)}{\partial t} = - \nabla \cdot ({\bf U}\rho h) + \psi + (\widetilde{{\bf U}}\cdot {\bf e}_r) \frac{\partial p_0}{\partial r}. \end{equation} \tag{ A38 }$$

  For spherical geometry, we solve the analytically equivalent form

  $$\begin{equation} \frac{\partial (\rho h)}{\partial t} = - \nabla \cdot ({\bf U}\rho h) + \psi + \nabla \cdot (\widetilde{{\bf U}}p_0) - p_0 \nabla \cdot \widetilde{{\bf U}}, \end{equation} \tag{ A39 }$$

  which experience has shown to minimize the drift from thermodynamic equilibrium. The update consists of two steps:

  • i.  
    Compute the time-centered enthalpy edge state, (ρh)^{n + 1/2,⋆,pred}, for the conservative update of (ρh)⁽¹⁾. We use Equation (18) to predict (ρh)' = (ρh)⁽¹⁾ − (ρh)ⁿ₀ to time-centered edges, using ${\bf U}= \widetilde{{\bf U}}^{{\rm ADV},\star }+w_0^{{n+1/2},\star } {\bf e}_r$, yielding $(\rho h)^{^{\prime }{n+1/2},\star,{\rm pred}}$. We convert the perturbational enthalpy to a full state enthalpy using

    $$\begin{equation} (\rho h)^{{n+1/2},\star,{\rm pred}} = (\rho h)^{^{\prime }{n+1/2},\star,{\rm pred}} + \frac{(\rho h)_0^n + (\rho h)_0^{n+1,\star }}{2}. \end{equation} \tag{ A40 }$$

    For planar geometry, we map (ρh)₀ directly to Cartesian edges. In spherical geometry, our experience has shown that a slightly different approach leads to reduced discretization errors. We first map h₀ ≡ (ρh)₀/ρ₀ and ρ₀ to Cartesian edges separately, and then multiply these terms to get (ρh)₀.
  • ii.  
    Evolve (ρh)⁽¹⁾ → (ρh)^(2),⋆.For planar geometry,

    $$\begin{eqnarray} (\rho h)^{(2),\star } &=& (\rho h)^{(1)} - \Delta t\big\lbrace \nabla \cdot \big[ \big(\widetilde{{\bf U}}^{{\rm ADV},\star }+w_0^{{n+1/2},\star } {\bf e}_r\big) (\rho h)^{{n+1/2},\star,{\rm pred}} \big] \big\rbrace + \Delta t\big(\widetilde{{\bf U}}^{{\rm ADV},\star }\cdot {\bf e}_r\big) \left(\frac{\partial p_0}{\partial r} \right)^{n} + \Delta t\psi ^{{n+1/2},\star }.\nonumber\\ \end{eqnarray} \tag{ A41 }$$

    For spherical geometry,

    $$\begin{eqnarray} (\rho h)^{(2),\star } &=& (\rho h)^{(1)} - \Delta t\big\lbrace \nabla \cdot \big[ \big(\widetilde{{\bf U}}^{{\rm ADV},\star }+w_0^{{n+1/2},\star } {\bf e}_r\big) (\rho h)^{{n+1/2},\star,{\rm pred}} \big] \big\rbrace \nonumber\\ && + \Delta t\big\lbrace \nabla \cdot \big(\widetilde{{\bf U}}^{{\rm ADV},\star }p_0^{n} \big) - p_0^{n} \nabla \cdot \widetilde{{\bf U}}^{{\rm ADV},\star }\big\rbrace + \Delta t\psi ^{{n+1/2},\star }. \end{eqnarray} \tag{ A42 }$$

    Then, for each Cartesian cell where ρ^(2),⋆ < ρ_cutoff, we recompute enthalpy using

    $$\begin{equation} (\rho h)^{(2),\star } = \rho ^{(2),\star }h\big(\rho ^{(2),\star },p_0^{n+1,\star },X_k^{(2),\star }\big). \end{equation} \tag{ A43 }$$
• I.  
  Update the temperature using the equation of state: T^(2),⋆ = T(ρ^(2),⋆, h^(2),⋆, X^(2),⋆_k) (planar geometry) or T^(2),⋆ = T(ρ^(2),⋆, p^{n + 1,⋆}₀, X^(2),⋆_k) (spherical geometry).

Step 5. React the full state through a second time interval of Δt/2.

Call React State$\big[ \rho ^{(2),\star },(\rho h)^{(2),\star }, X_k^{(2),\star }, T^{(2),\star }, (\rho H_{\rm ext})^{(2),\star }, p_0^{n+1,\star }\big] \rightarrow \big[ \rho ^{n+1,\star },(\rho h)^{n+1,\star }, X_k^{n+1,\star }, T^{n+1,\star }, (\rho \dot{\omega }_k)^{(2),\star }, (\rho H_{\rm nuc})^{(2),\star } \big].$

Step 6. Compute the time-centered expansion, S^{n + 1/2,⋆}, base state velocity, w^{n + 1/2}₀, and base state velocity forcing.

• A.  
  Compute S^{n + 1/2,⋆}. First, compute S^{n + 1,⋆} with

  $$\begin{equation} S^{n+1,\star } = -\sigma \sum _k \xi _k (\dot{\omega }_k)^{(2),\star } + \frac{1}{\rho ^{n+1,\star } p_\rho } \sum _k p_{X_k} ({\dot{\omega }}_k)^{(2),\star } + \sigma H_{\rm nuc}^{(2),\star } + \sigma H_{\rm ext}^{(2),\star }, \end{equation} \tag{ A44 }$$

  where $(\dot{\omega }_k)^{(2),\star } = (\rho \dot{\omega }_k)^{(2),\star } / \rho ^{(2),\star }$ and the thermodynamic quantities are defined using ρ^{n + 1,⋆}, X^{n + 1,⋆}_k, and T^{n + 1,⋆} as inputs to the equation of state. Then, define

  $$\begin{equation} \overline{S^{{n+1/2},\star }} = {\rm \bf Avg} (S^{{n+1/2},\star }), \qquad S^{{n+1/2}.\star } = \frac{S^n + S^{n+1,\star }}{2}. \end{equation} \tag{ A45 }$$
• B.  
  Compute w^{n + 1/2}₀. First, define

  $$\begin{equation} \overline{\Gamma _1^{{n+1/2},\star }} = \frac{\overline{\Gamma _1^n} + \overline{\Gamma _1^{n+1,\star }}}{2}, \qquad \rho _0^{{n+1/2},\star } = \frac{\rho _0^{n} + \rho _0^{n+1,\star }}{2}, \qquad p_0^{{n+1/2},\star } = \frac{p_0^{n} + p_0^{n+1,\star }}{2}, \end{equation} \tag{ A46 }$$

  with

  $$\begin{equation} \overline{\Gamma _1^{n+1,\star }} = {\rm {\bf Avg}} \big[ \Gamma _1\big(\rho ^{n+1,\star }, p_0^{n+1,\star }, X_k^{n+1,\star }\big) \big]. \end{equation} \tag{ A47 }$$

  For planar geometry, callCompute w₀ Planar$\big[\overline{S^{{n+1/2},\star }},\overline{\Gamma _1^{{n+1/2},\star }},p_0^{{n+1/2},\star },\psi ^{{n+1/2},\star }\big]\rightarrow \big[w_0^{{n+1/2}}\big]$.For spherical geometry, callCompute w₀ Spherical$\big[\overline{S^{{n+1/2},\star }},\overline{\Gamma _1^{{n+1/2},\star }},\rho _0^{{n+1/2},\star },p_0^{{n+1/2},\star },\eta _\rho ^{{n+1/2},\star }\big] \rightarrow \big[w_0^{{n+1/2}}\big]$.
• C.  
  Compute the base state velocity forcing. Rearrange Equation (A22),

  $$\begin{equation} \left(\frac{1}{\rho _0} \frac{\partial \pi _0}{\partial r} \right)^n = -\frac{w_0^{{n+1/2}} - w_0^{n-1/2}}{1/2(\Delta t^n+\Delta t^{n-1})} - w_0^n \left(\frac{\partial w_0}{\partial r}\right)^n, \end{equation} \tag{ A48 }$$

  where wⁿ₀ and (∂w₀/∂r)ⁿ are defined as

  $$\begin{eqnarray} w_0^n &=& \frac{\Delta t^{n} w_0^{{n-1/2}} + \Delta t^{n-1} w_0^{{n+1/2}}}{\Delta t^n+\Delta t^{n-1}}, \end{eqnarray} \tag{ A49 }$$

  $$\begin{eqnarray} \left(\frac{\partial w_0}{\partial r}\right)^{n} &=& \frac{1}{\Delta t^n+\Delta t^{n-1} } \left[ \Delta t^{n} \left(\frac{\partial w_0 }{ \partial r}\right)^{{n-1/2}} + \Delta t^{n-1} \left(\frac{\partial w_0 }{ \partial r}\right)^{{n+1/2}} \right]. \end{eqnarray} \tag{ A50 }$$

  If n = 0, we use Δt⁻¹ = Δt⁰.

Step 7. Construct the time-centered advective velocity on edges, $\widetilde{{\bf U}}^{\rm ADV}$.

The procedure to construct $\widetilde{{\bf U}}^{{\rm ADV},\dagger }$ is identical to the procedure for computing $\widetilde{{\bf U}}^{{\rm ADV},\dagger,\star }$ in Step 3, but uses the updated values w^{n + 1/2}₀ and πⁿ₀ rather than w^{n + 1/2,⋆}₀ and π^n,⋆₀. We note that $\widetilde{{\bf U}}^{\rm ADV}$ satisfies the discrete version of $\overline{(\widetilde{{\bf U}}^{\rm ADV}\cdot {\bf e}_r)}=0$ as well as

$$\begin{equation} \nabla \cdot \left(\beta _0^{{n+1/2},\star } \widetilde{{\bf U}}^{\rm ADV}\right) = \beta _0^{{n+1/2},\star }\left(S^{{n+1/2},\star } - \overline{S^{{n+1/2},\star }}\right), \end{equation} \tag{ A51 }$$

$$\begin{equation} \beta _0^{{n+1/2},\star } = \frac{ \beta _0^n + \beta _0^{n+1,\star } }{2}; \qquad \beta _0^{n+1,\star } = \beta _0 \left(\rho _0^{n+1,\star }, p_0^{n+1,\star }, \overline{\Gamma _1^{n+1,\star }}\right). \end{equation} \tag{ A52 }$$

Step 8. Advect the base state and full state through a time interval of Δt.

• A.  
  Update ρ₀, saving the time-centered density at radial edges by callingAdvect Base Density [ρⁿ₀, w^{n + 1/2}₀] → [ρ^(2a)₀, ρ^{n + 1/2,pred}₀].
• B.  
  Update (ρX_k). This step is identical to Step 4B except we use the updated values $w_0^{{n+1/2}}, \widetilde{{\bf U}}^{\rm ADV}$, and ρ^(2a)₀ rather than $w_0^{{n+1/2},\star }, \widetilde{{\bf U}}^{{\rm ADV},\star }$, and ρ^(2a),⋆₀. In particular,

  • i.  
    Compute the time-centered species edge states, (ρX_k)^{n + 1/2,pred}, for the conservative update of (ρX_k)⁽¹⁾. We use equations (17) and (13) to predict $\rho ^{^{\prime }(1)} = \rho ^{(1)} - \rho _0^n$ and X⁽¹⁾_k = (ρX_k)⁽¹⁾/ρ⁽¹⁾ to time-centered edges with ${\bf U}= \widetilde{{\bf U}}^{\rm ADV}+w_0^{{n+1/2}} {\bf e}_r$, yielding $\rho ^{^{\prime }{n+1/2},{\rm pred}}$ and X^{n + 1/2,pred}_k. We convert the perturbational density to a full state density using

    $$\begin{equation} \rho ^{{n+1/2},{\rm pred}} = \rho ^{^{\prime }{n+1/2},{\rm pred}} + \frac{\rho _0^n + \rho _0^{(2a)}}{2}. \end{equation} \tag{ A53 }$$

    Then, (ρX_k)^{n + 1/2,pred} = (ρ^{n + 1/2,pred}X^{n + 1/2,pred}_k).
  • ii.  
    Evolve (ρX_k)⁽¹⁾ → (ρX_k)⁽²⁾ using

    $$\begin{equation} (\rho X_k)^{(2)} = (\rho X_k)^{(1)} - \Delta t\left\lbrace \nabla \cdot \left[\left(\widetilde{{\bf U}}^{\rm ADV}+w_0^{{n+1/2}} {\bf e}_r\right) (\rho X_k)^{{n+1/2},{\rm pred}} \right] \right\rbrace, \end{equation} \tag{ A54 }$$

    $$\begin{equation} \rho ^{(2)} = \sum _k (\rho X_k)^{(2)}, \qquad X_k^{(2)} = (\rho X_k)^{(2)} / \rho ^{(2)}. \end{equation} \tag{ A55 }$$
• C.  
  Define a radial edge-centered η^{n + 1/2}_ρ.For planar geometry,

  $$\begin{eqnarray} \eta _\rho ^{{n+1/2}} &=& {\rm {\bf Avg}} \sum _k \left[\left(\widetilde{{\bf U}}^{\rm ADV}\cdot {\bf e}_r + w_0^{{n+1/2}}\right) (\rho X_k)^{{n+1/2},{\rm pred}} \right] - w_0^{{n+1/2}} \rho _0^{{n+1/2},{\rm pred}}. \end{eqnarray} \tag{ A56 }$$

  For spherical geometry, first construct η^{cart,n + 1/2}_ρ = [ρ'(U · e_r)]^{n + 1/2} on Cartesian cell centers using

  $$\begin{equation} \eta _\rho ^{{\rm cart},{n+1/2}} = \left[\left(\frac{\rho ^{(1)}+\rho ^{(2)}}{2}\right)-\left(\frac{\rho _0^n+\rho _0^{(2a)}}{2}\right)\right] \left(\widetilde{{\bf U}}^{\rm ADV}\cdot {\bf e}_r + w_0^{{n+1/2}}\right). \end{equation} \tag{ A57 }$$

  Then,

  $$\begin{equation} \eta _\rho ^{{n+1/2}} = {\rm {\bf Avg}}\left(\eta _\rho ^{{\rm cart},{n+1/2}}\right). \end{equation} \tag{ A58 }$$

  This gives a radial cell-centered η^{n + 1/2}_ρ. To get η^{n + 1/2}_ρ at radial edges, average the two neighboring cell-centered values.
• D.  
  Correct ρ₀ by setting ρ^{n + 1}₀ = Avg(ρ⁽²⁾).
• E.  
  Update p₀ by calling Enforce HSE[pⁿ₀, ρ^{n + 1}₀] → [p^{n + 1}₀].
• F.  
  Compute ψ^{n + 1/2}.For planar geometry,

  $$\begin{equation} \psi _j^{{n+1/2}} = \case{1}{2} \big(\eta _{\rho,j-1/2}^{{n+1/2}} + \eta _{\rho,j+1/2}^{{n+1/2}}\big) g. \end{equation} \tag{ A59 }$$

  For spherical geometry, first compute

  $$\begin{equation} \overline{\Gamma _1^{(2)}} = {\rm {\bf Avg}} \big[ \Gamma _1\big(\rho ^{(2)}, p_0^{n+1}, X_k^{(2)}\big) \big]. \end{equation} \tag{ A60 }$$

  Then, define ψ^{n + 1/2} using Equation (A18):

  $$\begin{eqnarray} \psi _j^{{n+1/2}} &=& \left(\frac{\overline{\Gamma _1^{(1)}}+\overline{\Gamma _1^{(2)}}}{2}\right)_j \left(\frac{p_0^n+p_0^{n+1}}{2}\right)_j \left\lbrace \overline{S_j^{{n+1/2}}} - \frac{1}{r_j^2} \left[\big(r^2 w_0^{{n+1/2}}\big)_{j+1/2} - \big(r^2 w_0^{{n+1/2}}\big)_{j-1/2} \right] \right\rbrace. \end{eqnarray} \tag{ A61 }$$
• G.  
  Update (ρh)₀ by calling Advect Base Enthalpy[(ρh)ⁿ₀, w^{n + 1/2}₀, ψ^{n + 1/2}] → [(ρh)^{n + 1}₀].
• H.  
  Update the enthalpy. This step is identical to Step 4H except we use the updated values $w_0^{{n+1/2}}, \widetilde{{\bf U}}^{\rm ADV}, \rho _0^{n+1}, (\rho h)_0^{n+1}, p_0^{n+1/2}$, and ψ^{n + 1/2} rather than$w_0^{{n+1/2},\star }, \widetilde{{\bf U}}^{{\rm ADV},\star }, \rho _0^{n+1,\star }, (\rho h)_0^{n+1,\star }, p_0^n$, and ψ^{n + 1/2,⋆}. In particular,

  • i.  
    Compute the time-centered enthalpy edge state, (ρh)^{n + 1/2,pred}, for the conservative update of (ρh)⁽¹⁾. We use Equation (18) to predict (ρh)' = (ρh)⁽¹⁾ − (ρh)ⁿ₀ to time-centered edges with ${\bf U}= \widetilde{{\bf U}}^{\rm ADV}+w_0^{{n+1/2}} {\bf e}_r$, yielding $(\rho h)^{^{\prime }{n+1/2},{\rm pred}}$. We convert the perturbational enthalpy to a full state enthalpy using

    $$\begin{equation} (\rho h)^{{n+1/2},{\rm pred}} = (\rho h)^{^{\prime }{n+1/2},{\rm pred}} + \frac{(\rho h)_0^n + (\rho h)_0^{n+1}}{2}. \end{equation} \tag{ A62 }$$
  • ii.  
    Evolve (ρh)⁽¹⁾ → (ρh)⁽²⁾.For planar geometry,

    $$\begin{eqnarray} (\rho h)^{(2)} &=& (\rho h)^{(1)} - \Delta t\big\lbrace \nabla \cdot \big[ \big(\widetilde{{\bf U}}^{\rm ADV}+w_0^{{n+1/2}} {\bf e}_r\big) (\rho h)^{{n+1/2},{\rm pred}} \big] \big\rbrace + \Delta t\big(\widetilde{{\bf U}}^{\rm ADV}\cdot {\bf e}_r\big) \left(\frac{\partial p_0}{\partial r} \right)^{n+1/2}+ \Delta t\psi ^{{n+1/2}}.\nonumber\\ \end{eqnarray} \tag{ A63 }$$

    For spherical geometry,

    $$\begin{eqnarray} (\rho h)^{(2)} &=& (\rho h)^{(1)} - \Delta t\big\lbrace \nabla \cdot \big[ \big(\widetilde{{\bf U}}^{\rm ADV}+w_0^{{n+1/2}} {\bf e}_r\big) (\rho h)^{{n+1/2},{\rm pred}} \big] \big\rbrace \nonumber\\ && + \Delta t\big[ \nabla \cdot \big(\widetilde{{\bf U}}^{\rm ADV}p_0^{{n+1/2}} \big) - p_0^{{n+1/2}} \nabla \cdot \widetilde{{\bf U}}^{\rm ADV}\big] + \Delta t\psi ^{{n+1/2}}, \end{eqnarray} \tag{ A64 }$$

    where p^{n + 1/2}₀ is defined as p^{n + 1/2}₀ = (pⁿ₀ + p^{n + 1}₀)/2.Then, for each Cartesian cell where ρ⁽²⁾ < ρ_cutoff, we recompute enthalpy using

    $$\begin{equation} (\rho h)^{(2)} = \rho ^{(2)}h\big(\rho ^{(2)},p_0^{n+1},X_k^{(2)}\big). \end{equation} \tag{ A65 }$$
• I.  
  Update the temperature using the equation of state: T⁽²⁾ = T(ρ⁽²⁾, h⁽²⁾, X⁽²⁾_k) (planar geometry) or T⁽²⁾ = T(ρ⁽²⁾, p^{n + 1}₀, X⁽²⁾_k) (spherical geometry).

Step 9. React the full state through a second time interval of Δt/2.

Call React State$[\rho ^{(2)},(\rho h)^{(2)}, X_k^{(2)},T^{(2)}, (\rho H_{\rm ext})^{(2)}, p_0^{n+1}] \rightarrow [\rho ^{n+1}, (\rho h)^{n+1}, X_k^{n+1}, T^{n+1}, (\rho \dot{\omega }_k)^{(2)}, (\rho H_{\rm nuc})^{(2)} ].$

Step 10. Define the new time expansion, S^{n + 1}, and $\overline{\Gamma _1^{n+1}}$.

• A.  
  Define

  $$\begin{equation} S^{n+1} = -\sigma \sum _k \xi _k (\dot{\omega }_k)^{(2)} + \sigma H_{\rm nuc}^{(2)} + \frac{1}{\rho ^{n+1} p_\rho } \sum _k p_{X_k} ({\dot{\omega }}_k)^{(2)} + \sigma H_{\rm ext}^{(2)}, \end{equation} \tag{ A66 }$$

  where $(\dot{\omega }_k)^{(2)} = (\rho \dot{\omega }_k)^{(2)} / \rho ^{(2)}$ and the thermodynamic quantities are defined using ρ^{n + 1}, X^{n + 1}_k, and T^{n + 1} as inputs to the equation of state. Then, compute

  $$\begin{equation} \overline{S^{n+1}} = {\rm \bf Avg} (S^{n+1}). \end{equation} \tag{ A67 }$$
• B.  
  Define

  $$\begin{equation} \overline{\Gamma _1^{n+1}} = {\rm {\bf Avg}}\left[\Gamma _1\left(\rho ^{n+1}, p_0^{n+1}, X_k^{n+1}\right) \right]. \end{equation} \tag{ A68 }$$

Step 11. Update the velocity.

First, we compute the time-centered edge velocities, $\widetilde{{\bf U}}^{{n+1/2},{\rm pred}}$. Then, we define

$$\begin{equation} \rho ^{n+1/2}= \frac{\rho ^n + \rho ^{n+1}}{2}, \qquad \rho _0^{n+1/2}= \frac{\rho _0^n + \rho _0^{n+1}}{2}. \end{equation} \tag{ A69 }$$

We update the velocity field $\widetilde{{\bf U}}^n$ to $\widetilde{{\bf U}}^{n+1,\dagger }$ by discretizing Equation (10) as

$$\begin{eqnarray} \widetilde{{\bf U}}^{n+1,\dagger } &=& \widetilde{{\bf U}}^n - \Delta t\big[\big(\widetilde{{\bf U}}^{\rm ADV}+ w_0^{{n+1/2}} {\bf e}_r\big) \cdot \nabla \widetilde{{\bf U}}^{{n+1/2},{\rm pred}} \big] - \Delta t\big(\widetilde{{\bf U}}^{\rm ADV}\cdot {\bf e}_r\big) \left(\frac{\partial w_0}{\partial r} \right)^{n+1/2}{\bf e}_r \nonumber \\ && + \Delta t\left[ - \frac{1}{\rho ^{n+1/2}} \mathbf {G} \pi ^{n-1/2}+ \left(\frac{1}{\rho _0}\frac{\partial \pi _0}{\partial r}\right)^n {\bf e}_r - \frac{\big(\rho ^{n+1/2}-\rho _0^{n+1/2}\big)}{\rho ^{n+1/2}} g^{{n+1/2}} {\bf e}_r \right], \end{eqnarray} \tag{ A70 }$$

where G approximates a cell-centered gradient from nodal data. Again, the † superscript refers to the fact that the updated velocity does not satisfy the divergence constraint.

Finally, we use an approximate nodal projection to define $\widetilde{{\bf U}}^{n+1}$ from $\widetilde{{\bf U}}^{n+1,\dagger },$ such that $\widetilde{{\bf U}}^{n+1}$ approximately satisfies

$$\begin{equation} \nabla \cdot \big(\beta _0^{{n+1/2}} \widetilde{{\bf U}}^{n+1} \big) = \beta _0^{{n+1/2}} \big(S^{n+1} - \overline{S^{n+1}} \big), \end{equation} \tag{ A71 }$$

where β^{n + 1/2}₀ is defined as

$$\begin{equation} \beta _0^{{n+1/2}} = \frac{\beta _0^n + \beta _0^{n+1}}{2}; \qquad \beta _0^{n+1} = \beta \big(\rho _0^{n+1}, p_0^{n+1}, \overline{\Gamma _1^{n+1}}, g^{n+1}\big). \end{equation} \tag{ A72 }$$

As part of the projection we also define the new-time perturbational pressure, π^{n + 1/2}. This projection necessarily differs from the MAC projection used in Step 3 and Step 7 because the velocities in those steps are defined on edges and $\widetilde{{\bf U}}^{n+1}$ is defined at cell centers, requiring different divergence and gradient operators. Details of the approximate projection are given in Paper III.

Step 12. Compute a new Δt.

Compute Δt for the next time step with the procedure described in Section 3.4 of Paper III using w₀ as computed in Step 6 and $\widetilde{{\bf U}}^{n+1}$ as computed in Step 11.

This completes one step of the algorithm.

 Figure 14 is a flow chart summarizing the 12 step algorithm, including the notation used as we advance the solution by Δt. Figure 15 is a flow chart of the advection steps (Steps 4 and 8), which includes the notation we use as we advect the solution through a time interval of Δt.

As discussed in Paper IV, in order to prevent the velocity from becoming too large in low density regions far from the center of the star, we impose a cutoff at a moderately small density, ρ_cutoff, and hold the density at this constant value outside of the star. The cutoff affects the evolution in the following ways:

• 1.  
  After advancing enthalpy in the advection step (Steps 4Hii and 8Hii in Section A.4), we recompute the enthalpy using the equation of state if ρ ⩽ ρ_cutoff.
• 2.  
  When computing gravity (Section A.2) we only add ρ₀ to m_encl if ρ₀ > ρ_cutoff in order to prevent an unphysical amount of mass from contributing to the calculation.
• 3.  
  When computing p₀ in Enforce HSE (Section A.3.3), we hold p₀ constant once ρ₀ ⩽ ρ_cutoff.
• 4.  
  In React State (Section A.3.1), we set $\dot{\omega }_k=0$ and ρH_nuc = 0 if ρ ⩽ ρ_cutoff.
• 5.  
  When computing ψ (Steps 4F and 8F in Section A.4), we set ψ = 0 if ρ₀ ⩽ ρ_cutoff.
• 6.  
  When we compute the velocity forcing in Steps 3, 7, 11, and 12, we set the buoyancy term (the term proportional to ρ − ρ₀) to zero if ρ < 5ρ_cutoff.

Additionally, we use an anelastic cutoff density, ρ_anelastic, in the computation of β₀ (Steps 3, 7, and 11 in Section A.4). When ρ_0,j ⩽ ρ_anelastic, we set β_0,j = (ρ_0,j/ρ_{0,j − 1})β_{0,j − 1}.