MAESTRO: An Adaptive Low Mach Number Hydrodynamics Algorithm for Stellar Flows

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

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{\mathchoice{\kern 0.0pt\raise 2.15277pt\hbox{\sevenrm 3}\kern-1.49994pt/\kern-1.49994pt\lower 1.07639pt\hbox{\sevenrm 4}}{\kern 0.0pt\raise 2.15277pt\hbox{\sevenrm 3}\kern-1.49994pt/\kern-1.49994pt\lower 1.07639pt\hbox{\sevenrm 4}}{\kern 0.0pt\raise 2.15277pt\hbox{\fiverm 3}\kern-1.99997pt/\kern-1.49994pt\lower 1.07639pt\hbox{\fiverm 4}}{3\!/4}+2m}$ from the center of the star, where $m\geq 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{\mathchoice{\kern 0.0pt\raise 2.15277pt\hbox{\sevenrm 3}\kern-1.49994pt/\kern-1.49994pt\lower 1.07639pt\hbox{\sevenrm 4}}{\kern 0.0pt\raise 2.15277pt\hbox{\sevenrm 3}\kern-1.49994pt/\kern-1.49994pt\lower 1.07639pt\hbox{\sevenrm 4}}{\kern 0.0pt\raise 2.15277pt\hbox{\fiverm 3}\kern-1.99997pt/\kern-1.49994pt\lower 1.07639pt\hbox{\fiverm 4}}{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. 1.

   Create an itemized list, $\hat{\phi}_{m}$, where each element is associated with a radius $\hat{r}_{m}=\Delta x\sqrt{\mathchoice{\kern 0.0pt\raise 2.15277pt\hbox{\sevenrm 3}\kern-1.49994pt/\kern-1.49994pt\lower 1.07639pt\hbox{\sevenrm 4}}{\kern 0.0pt\raise 2.15277pt\hbox{\sevenrm 3}\kern-1.49994pt/\kern-1.49994pt\lower 1.07639pt\hbox{\sevenrm 4}}{\kern 0.0pt\raise 2.15277pt\hbox{\fiverm 3}\kern-1.99997pt/\kern-1.49994pt\lower 1.07639pt\hbox{\fiverm 4}}{3\!/4}+2m}$ from the center of the star.

2. 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. 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:

   $$\overline{\phi(r)}=\hat{\phi}_{k_{-}}+\frac{\hat{\phi}_{k}-\hat{\phi}_{k_{-}}}{\hat{r}_{k}-\hat{r}_{k_{-}}}(r-\hat{r}_{k_{-}})+\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}),$$

   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_{+}}$. 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.

In §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.

planar:

1. 1.

   To map a cell-centered radial array onto Cartesian cell centers, we use direct-injection since the radial cell centers are in alignment with the Cartesian cell centers.

2. 2.

   To map an edge-centered radial array onto Cartesian cell centers, we average the two nearest radial edge-centered values.

3. 3.

   To map a cell-centered radial array onto Cartesian edges with normal in the radial direction, we use 4th order spatial interpolation. For example, in two dimensions,

   $$\phi_{i,j+\mathchoice{\kern 0.0pt\raise 1.50694pt\hbox{\sevenrm 1}\kern-1.22911pt/\kern-1.22911pt\lower 0.75346pt\hbox{\sevenrm 2}}{\kern 0.0pt\raise 1.50694pt\hbox{\sevenrm 1}\kern-1.22911pt/\kern-1.22911pt\lower 0.75346pt\hbox{\sevenrm 2}}{\kern 0.0pt\raise 1.50694pt\hbox{\fiverm 1}\kern-1.63885pt/\kern-1.22911pt\lower 0.75346pt\hbox{\fiverm 2}}{1\!/2}}=\frac{7}{12}(\phi_{j}+\phi_{j+1})-\frac{1}{12}(\phi_{j-1}+\phi_{j+2}).$$

   (21)

   We constrain $\phi_{i,j+\mathchoice{\kern 0.0pt\raise 1.50694pt\hbox{\sevenrm 1}\kern-1.22911pt/\kern-1.22911pt\lower 0.75346pt\hbox{\sevenrm 2}}{\kern 0.0pt\raise 1.50694pt\hbox{\sevenrm 1}\kern-1.22911pt/\kern-1.22911pt\lower 0.75346pt\hbox{\sevenrm 2}}{\kern 0.0pt\raise 1.50694pt\hbox{\fiverm 1}\kern-1.63885pt/\kern-1.22911pt\lower 0.75346pt\hbox{\fiverm 2}}{1\!/2}}$ to lie between the interpolated values, and lower the order of interpolation near domain boundaries. For the Cartesian edges transverse to the base state direction, we use direct-injection since the radial cell centers are in alignment with these Cartesian edges.

4. 4.

   To map an edge-centered radial array onto Cartesian edges, we use direct-injection on Cartesian edges normal to the base state direction since the radial edges are in alignment with these Cartesian edges. For the remaining Cartesian edges, we average the two nearest radial edge-centered values.

spherical:

1. 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 4a). This is a departure from Paper IV, in which we used piecewise constant interpolation.

2. 2.

   To map an edge-centered radial array onto Cartesian cell centers, we use linear interpolation from the nearest two points (see Figure 4b).

3. 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 4c).

4. 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 4c).

At each time step the state data is defined on a nested hierarchy of grids, ranging from the base level ($\ell=1$), which covers the entire computational domain, to the finest level ($\ell={\ell}_{\mathrm{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 ($\Delta x=\Delta y=\Delta 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 $\ell+1$ is completely contained in the union of grids at level $\ell$. Additionally, in the interior, we require that each grid at level $\ell+1$ be a distance of at least two level $\ell$ cells from the boundary between level $\ell$ and level $\ell-1$ grids; this allows us to always fill “ghost cells” at level $\ell+1$ from the level $\ell$ 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,¹¹1The “Step” notation is used in describing the full details of the algorithm in Appendix A.4 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 Step 3 and 7 enforce that $\widetilde{{\bf{U}}}^{\mathrm{ADV},\star}$ and $\widetilde{{\bf{U}}}^{\mathrm{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 $\phi$ 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:

• •

  Before any explicit operation at level $\ell>1,$ data in ghost cells at that level are filled by interpolating from level $\ell-1$, or imposing physical boundary conditions, as appropriate.

• •

  Edge-based fluxes at level $\ell<\ell_{max}$ that underlie edges at level $\ell+1$ are defined to be the average of the fluxes on level $\ell+1$ at that edge. This enforces conservation.

• •

  After any update to the solution, data at finer levels is 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 vs. 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 $\Delta r\neq\Delta x$. Both of these methods pose problems, as they generate oscillations in perturbational quantities (such as $\rho^{\prime},(\rho h)^{\prime}$ 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.

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 ²²2The boldface notation refers to numerical modules we have described in Appendix A.3. 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_{0}$ using Enforce HSE. To compute $\psi$ and $w_{0}$ 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 $\Delta r$ for all levels (see Figure 5). As in the single-level algorithm, we choose $\Delta r=\Delta x/5$, but here, $\Delta 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 §6.6, as an example). In §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(\rho,h,X_{k})$ (for planar problems) or $T=T(\rho,p_{0},X_{k})$ (for spherical problems). Then, we recompute $\overline{\Gamma_{1}}$ and $\beta_{0}$ 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^{3}$ 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 §4, we first create a unit cube with $384^{3}$ resolution and no refinement. We create a radial array with $\Delta r=\Delta x/5$, and initialize the radial array cell centers with the Gaussian profile $\phi_{\rm exact}(r)=e^{-10r^{2}}$. We map $\phi_{\rm exact}$ to Cartesian cell centers with the fill operation, then compute the lateral average of this Cartesian field, ${\bf Avg}(\phi)$. We repeat this process by choosing the grid with two levels of refinement used in the full-star simulation test in §6.6, shown in Figure 12.

Figure 6 (left) shows the relative error between $\phi_{\rm exact}$ and ${\bf Avg}(\phi$) 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\in[0,0.036]$ and is $\mathcal{O}(10^{-13})$ for $r\in[0.036,0.7]$. For the test with two levels of refinement, the relative error is $\mathcal{O}(10^{-8})$ for $r\in[0,0.7]$, which is still a vast improvement when compared to the Paper IV mapping applied to a single-level simulation.

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\times 10^{9}~\mathrm{g~cm^{-3}}$), temperature ($6\times 10^{8}$ K), and a uniform composition ($X(^{12}\mathrm{C})=0.3$, $X(^{16}\mathrm{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:

We begin with a guess of $\rho_{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+\mathchoice{\kern 0.0pt\raise 1.50694pt\hbox{\sevenrm 1}\kern-1.22911pt/\kern-1.22911pt\lower 0.75346pt\hbox{\sevenrm 2}}{\kern 0.0pt\raise 1.50694pt\hbox{\sevenrm 1}\kern-1.22911pt/\kern-1.22911pt\lower 0.75346pt\hbox{\sevenrm 2}}{\kern 0.0pt\raise 1.50694pt\hbox{\fiverm 1}\kern-1.63885pt/\kern-1.22911pt\lower 0.75346pt\hbox{\fiverm 2}}{1\!/2}}$, is updated each iteration based on the current value of the density. Once the temperature falls below $10^{7}$ 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\times 10^{8}$ cm of the model onto a one-dimensional array with 1280 elements with $\Delta r=3.90625\times 10^{5}$ cm. For the full-star three-dimensional simulation, we map the model onto a $5\times 10^{8}$ cm³ domain with 256³ Cartesian grid cells with $\Delta x=5\Delta r=19.53125\times 10^{5}$ cm. For the one and three-dimensional MAESTRO calculations, we use cutoff densities (see Appendix A.5) of $\rho_{\rm cutoff}=10^{5}~\mathrm{g~cm^{-3}}$ and $\rho_{\rm anelastic}=10^{6}~\mathrm{g~cm^{-3}}$, corresponding to radii of approximately $1.8\times 10^{8}$ cm and $1.9\times 10^{8}$ 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_{\rm sp}$ equal to the radius where $\rho_{0}=10^{7}~\mathrm{g~cm^{-3}}$, $r_{\rm md}$ equal to the radius where $\rho_{0}=3\times 10^{6}~\mathrm{g~cm^{-3}}$, and $\kappa=10$ s^-1. 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 seconds and look at the solution at 2.0 seconds (after the compressible solution has had time to re-equilibrate). We use $H_{\rm ext}=H_{0}e^{-(r/10^{7}~\mathrm{cm})^{2}}$, with $H_{0}=10^{16}~\mathrm{erg~g^{-1}~s^{-1}}$ 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 $\rho_{0},p_{0}$, 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.

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\times 10^{7}$ cm by $2.88\times 10^{8}$ cm and generate a high resolution initial model with $\Delta r=7.03125\times 10^{4}$ cm (which is equal to $\Delta x$ for our highest resolution, “exact” solution) using the method described in Appendix C with $r_{\mathrm{base}}=0$. For each of the remaining, lower resolution simulations for this test, we generate an initial model using $\Delta r$ equal to the effective $\Delta x$ of each simulation by linearly interpolating values from the high resolution model. Next, we add a temperature perturbation of the form:

where $\sigma=2.5\times 10^{6}~$cm and $d_{ij}$ is the physical distance between the cell center corresponding to cell $(i,j)$ and the location $(3.6\times 10^{7}~$cm$,3.2\times 10^{7}~$cm). Then, we call the equation of state to compute a consistent $\rho,h=\rho,h(p_{0},T,X_{k})$ everywhere. We use the reaction network described in § 4.2 in Paper III. We use cutoff densities of $\rho_{\rm cutoff}=\rho_{\rm anelastic}=3\times 10^{6}~\mathrm{g~cm^{-3}}$, and a sponge with $r_{\rm sp}$ equal to the radius where $\rho_{0}=10\rho_{\rm cutoff}$, $r_{\rm md}$ equal to the radius where $\rho_{0}=\rho_{\rm cutoff}$, and $\kappa=10$ s^-1. 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\times 4096$ cells and $\Delta t=3.125\times 10^{-3}$ s to be the exact solution. We perform three single-level simulations using resolutions of $64\times 256,128\times 512$, and $256\times 1024$ grid cells using fixed time steps of $\Delta 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\times 512$ and $256\times 1024$ grid cells with fixed time steps of $\Delta 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\in[1.8\times 10^{7},5.4\times 10^{7}]$ cm, so effectively we have refined 1/8^th 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_{0}$ from $\rho_{0}$, we use $r=1.8\times 10^{7}$ 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_{1}$ error norm for each simulation, we average the data from the exact solution onto a grid with corresponding resolution. We measure the $L_{1}$ error norm in the physical space corresponding to the refined region using

where $n_{\rm cell}$ is the number of cells we sum over. This form of the $L_{1}$ error norm gives us the average error per cell. We compute the convergence rate, $p$, between a coarser and finer simulation using

Tables 1 and 2 show the $L_{1}$ 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_{1}$ 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 §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\times 10^{7}$ cm by $7.2\times 10^{7}$ cm $2.88\times 10^{8}$ cm and for each simulation, we generate an initial model with $\Delta r=5.625\times 10^{5}$ cm (which is equal to the effective $\Delta x$ for both of the simulations in this test) using the method described in Appendix C with $r_{\mathrm{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\times 10^{7}~$cm$,3.6\times 10^{7}~$cm$,3.2\times 10^{7}~$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\times 128\times 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 seconds 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 seconds 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 second 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\times 10^{7}$ cm.

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\times 10^{8}$ cm by $4\times 10^{8}$ cm and for each simulation, we generate an initial model with $\Delta r$ equal to the effective $\Delta x$ for that simulation using the method described in Appendix C with $r_{\mathrm{base}}=5\times 10^{7}$ 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

with $r_{\rm layer}=1.25\times 10^{8}$ cm, $H_{0}=2.5\times 10^{16}$ erg g^-1 s^-1, $L_{x}=2.5\times 10^{8}$ 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 $\Psi=(0,\pi/3,\pi/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 $\rho_{\rm cutoff}=\rho_{\rm anelastic}=3\times 10^{6}~\mathrm{g~cm^{-3}}$ and a sponge with $r_{\rm sp}$ equal to the radius where $\rho_{0}=10^{8}~\mathrm{g~cm^{-3}}$, $r_{\rm md}$ equal to the radius where $\rho_{0}=\rho_{\rm cutoff}$, and $\kappa=100$ s^-1. 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\times 128$ and $320\times 512$ cells, and a simulation with two levels of refinement and an effective resolution of $320\times 512$ cells. For the multi-level simulation, we fix the refined grids, ensuring that $r\in[9.375\times 10^{7},1.5626\times 10^{8}]$ 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 $\rho_{0},p_{0}$, and $\overline{T}$ after 30 seconds 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 4 processors; the single-level fine grid run required approximately 1.9 seconds per time step and the multi-level run required approximately 0.6 seconds per time step, for a factor of 3 speedup.

Figure 11 shows the temperature profile after 3 and 4 seconds for each of the three simulations. The vertical distance shown is from $z=5\times 10^{7}$ cm to $2.2\times 10^{8}$ 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.

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 §2.4 in Paper IV. We map this one-dimensional model into the center of a computational domain of $5\times 10^{8}$ 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_{12}$ value of 2.93 describe in that paper. We initialize the simulation with a velocity perturbation described exactly as in §2.4 in Paper IV. We use the same cutoff densities, sponge parameters, and boundary conditions as in §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 $\rho>10^{8}~\mathrm{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.

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 non-linear, and the reaction rates scale with temperature as $\sim T^{23}$ (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 upwards 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 non-linear, 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 seconds 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 seconds 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.

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.