MAESTROeX: A Massively Parallel Low Mach Number Astrophysical Solver

The spatial discretization and AMR methodology remains unchanged from Paper V. We now summarize some of the key points here before describing the new temporal integrator in the next section. We recommend the reader review Section 3 of Paper V for further details.

We shall refer to local atmospheric flows in two and three dimensions as problems in "planar" geometry, and full-star flows in three dimensions as problems in "spherical" geometry. The solution in both cases consists of the Cartesian grid solution and the one-dimensional base state solution. Figure 1 illustrates the relationship between the base state and the Cartesian grid state for both planar and spherical geometries in the presence of spatially adaptive grids. One of the key numerical modules is the "lateral average," which computes the average over a layer of constant radius of a Cartesian grid variable and stores the result in a one-dimensional base state array. In planar geometries, this is a straightforward arithmetic average of values in cells at a particular height since the base state cell centers are in alignment with the Cartesian grid cell centers. However for spherical problems, the procedure is much more complicated. In Section 4 of Paper V, we describe how there is a finite, easily computable set of radii that any three-dimensional Cartesian cell center can map to. Specifically, for every three-dimensional Cartesian cell, there exists an integer m such that the distance from the cell center to the center of the star is given by

$$\begin{eqnarray}&&{\hat{r}}_{m}={\rm{\Delta }}x\sqrt{0.75+2m}.\end{eqnarray} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

Figure 2 is a two-dimensional illustration (two dimensions is chosen in the figure for ease of exposition; this mapping is only used for three-dimensional spherical stars) of the relationship between the Cartesian grid state and the one-dimensional base state array. We compute the lateral average by first summing the values in all the cells associated with each radius, dividing by the number of contributing cells to obtain the arithmetic average, and then using quadratic interpolation to map this data onto a one-dimensional base state. Previously, for spherical problems MAESTRO only allowed for a base state with constant Δr (typically equal to Δx/5).

Zoom In Zoom Out Reset image size

The companion "fill" module maps a base state array onto the full Cartesian grid state. For planar problems, direct injection can be used due to the perfect alignment of the base state and Cartesian grid state. For spherical problems, quadratic interpolation of the base state is used to assign values to each Cartesian cell center.

In this paper we explore a new option to retain an irregularly spaced base state to eliminate mapping errors from the "fill" module. For the lateral average, as before we first sum the values in all the cells associated with each radius and divide by the number of contributing cells to obtain the arithmetic average. However we do not interpolate this onto a uniformly spaced base state and retain the use of this irregularly spaced base state. The advantage with this approach is that the fill module does not require any interpolation. A potential benefit to eliminating the mapping error is to consider a spherical star in hydrostatic equilibrium at rest. In the absence of reactions, the star should remain at rest. The buoyancy forcing term in the momentum equation contains $\rho \mbox{--}{\rho }_{0}$. With the original scheme, interpolation errors in computing ρ₀ by averaging would cause artificial acceleration in the velocity field due to the interpolation error from the Cartesian grid to and from the radial base state. By retaining the radial base state as an irregularly spaced array, the truncation error in the mapping scheme is completely eliminated. The only source of error is due to machine precision effects resulting from averaging a large number of values, which is not significant over the course of any simulation. We note that Δr decreases as the base state moves further from the center of the star, which results in far more total cells in the irregularly spaced array than the previous uniformly spaced array.

We now describe the new temporal integration scheme, noting that it can be used for the original base state mapping (with uniform base state grid spacing) as well as the new irregularly spaced base state mapping. Previously we adopted an approach where we split the velocity into a base state component, ${w}_{0}(r,t)$, and a local velocity $\widetilde{{\boldsymbol{U}}}({\boldsymbol{x}},t)$, so that

$$\begin{eqnarray}&&{\boldsymbol{U}}=\widetilde{{\boldsymbol{U}}}({\boldsymbol{x}},t)+{w}_{0}(r,t){{\boldsymbol{e}}}_{r},\end{eqnarray} \tag{ 10 }$$

where ${{\boldsymbol{e}}}_{r}$ is the normal vector in the outward radial direction. We used w₀ to provide an estimate of the base state density evolution over a time step. This resulted in some unnecessary complications to the temporal integration scheme including base state advection modules for density, enthalpy, and velocity, as well as more cumbersome split velocity dynamics evolution equations. Our new temporal integration scheme uses full velocities for scalar and velocity advection, and only uses the above splitting to satisfy the velocity divergence constraint due to boundary considerations at the edge of the star. This results in a much simpler numerical scheme than the one from Paper V since we use the velocity directly rather than more complex terms involving the perturbational velocity. Additionally, the new scheme uses a simpler predictor-corrector approach to the base state density and pressure that no longer requires evolution equations and numerical discretizations to update the base state, greatly simplifying the algorithm while retaining the same overall second-order accuracy in time.

At the beginning of each time step we have the cell-centered Cartesian grid state, ${({\boldsymbol{U}},\rho {X}_{k},\rho h)}^{n}$, and nodal Cartesian grid state, ${\pi }^{n-1/2}$, and base state ${({\rho }_{0},{p}_{0})}^{n}$. At any time, the associated density, composition, and enthalpy can be trivially computed using, e.g.,

$$\begin{eqnarray}&&{\rho }^{n}=\displaystyle \sum _{k}{\left(\rho {X}_{k}\right)}^{n},\quad {X}_{k}^{n}={\left(\rho {X}_{k}\right)}^{n}{/\rho }^{n},\quad {h}^{n}={\left(\rho h\right)}^{n}{/\rho }^{n}.\end{eqnarray} \tag{ 11 }$$

Temperature is computed using the EOS,³ e.g., $T\,=T(\rho ,{p}_{0},{X}_{k})$, where p₀ has been mapped to the Cartesian grid using the fill module, and (${\overline{{\rm{\Gamma }}}}_{1},{\beta }_{0})$ are similarly computed from $(\rho ,{p}_{0},{X}_{k});$ see Appendix A of Paper I and Appendix C of Paper III for details on how β₀ is computed.

The overall flow of the algorithm begins with a second-order Strang splitting approach to integrate the advection-reaction system for the thermodynamic variables $(\rho {X}_{k},\rho h)$, followed by a second-order projection methodology to integrate the velocities subject to a divergence constraint. Within the thermodynamic variable update we use a predictor-corrector approach to achieve second-order accuracy in time. To summarize:

• 1.  
  In Step 1 we react the thermodynamic variables over the first Δt/2 interval.
• 2.  
  In Steps 2–4 we advect the thermodynamic variables over Δt. Specifically, we compute an estimate for the expansion term, S, compute face-centered, time-centered velocities that satisfy the divergence constraint, and then advect the thermodynamic variables.
• 3.  
  In Step 5 we react the thermodynamic variables over the second Δt/2 interval.⁴
• 4.  
  In Steps 6–8 we redo the advection in Steps 2–4 but are able to use the trapezoidal rule to time-center certain quantities such as S, ρ₀, etc.
• 5.  
  In Step 9 we redo the reactions from Step 5 beginning with the improved results from Steps 6–8.
• 6.  
  In Steps 10–11 we advect the velocity, and then project these velocities so they satisfy the divergence constraint while updating π.

There are a few key numerical modules we use in each time step.

• 1.  
  Average $[\phi ]\to [\overline{\phi }]$ computes the lateral average of a quantity over a layer at constant radius r, as described above in Section 3.1.
• 2.  
  Enforce HSE $[{\rho }_{0}]\to [{p}_{0}]$ computes the base state pressure, p₀, from a base state density, ρ₀ by integrating the hydrostatic equilibrium condition in one dimension. This follows equation (A10) in Paper V, noting that for the irregularly spaced base state case, Δr is not constant, where Δr_j+1/2 = r_j+1−r_j for cell face with index j+1/2. The base state pressure remains equal to a constant value at the location of a prescribed cutoff density outward for the entire simulation.
• 3.  
  React State $[{(\rho {X}_{k})}^{\mathrm{in}},\,{(\rho h)}^{\mathrm{in}},\,{p}_{0}]\,\to \,[{(\rho {X}_{k})}^{\mathrm{out}},\,{(\rho h)}^{\mathrm{out}},\,(\rho \dot{\omega }),(\rho {H}_{\mathrm{nuc}})]$ uses the multistep variable-coefficient Adams–Moulton and Backward Differentiation Formula methods in the VODE (Brown et al. 1989) package to integrate the species and enthalpy due to reactions over Δt/2 by solving

  $$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dX}}_{k}}{{dt}} & = & {\dot{\omega }}_{k}(\rho ,{X}_{k},T);\\ \displaystyle \frac{{dT}}{{dt}} & = & \displaystyle \frac{1}{{c}_{p}}\left(-\displaystyle \sum _{k}{\xi }_{k}{\dot{\omega }}_{k}+{H}_{\mathrm{nuc}})\right).\end{array}\end{eqnarray} \tag{ 12 }$$

  The inputs are the species, enthalpy, and base state pressure, and the outputs are the species, enthalpy, reaction rates, and nuclear energy generation rate. See Paper III for details.

Each time step is constrained by the standard advective CFL condition,

$$\begin{eqnarray}&&{\rm{\Delta }}t={\sigma }^{\mathrm{CFL}}\mathop{\min }\limits_{i}({\rm{\Delta }}x/{U}_{i}),\end{eqnarray} \tag{ 13 }$$

where for our simulations we typically use ${\sigma }^{\mathrm{CFL}}\sim 0.7$ and the minimum is taken over all spatial directions over all cells. For simulations that are initialized with zero velocity, we use the buoyancy force in the momentum equation to give a proper timescale to start the simulation; see Section 3.4 in Paper III for details.

In stratified low Mach number models, due to the extreme variation in density, the velocity can become very large in low density regions at the edge of the star. These large velocities can severely affect the time step, so throughout Papers II–V, we have employed two techniques to help control these dynamics without significantly affecting the dynamics in the convective region. First, we use a cutoff density technique, where we hold the density constant outside a specified radius (typically near where the density is ∼4 orders of magnitude smaller than the largest densities in the simulation). Second, we employ a sponge technique where we artificially damp the velocities near and beyond the cutoff region. For more details, refer to Paper V and the previous references cited within.

Beginning with ${({\boldsymbol{U}},\rho {X}_{k},\rho h)}^{n}$, ${\pi }^{n-1/2}$, and ${({\rho }_{0},{p}_{0})}^{n}$, the temporal integration scheme contains the following steps:

• Step 1: react the thermodynamic variables over the first Δt/2 interval. Call React State [(ρX_k)ⁿ, (ρh)ⁿ, ${p}_{0}^{n}$] → [${\left(\rho {X}_{k}\right)}^{\left(1\right)},{\left(\rho h\right)}^{\left(1\right)}$, ${\left(\rho {\dot{\omega }}_{k}\right)}^{\left(1\right)},{\left(\rho {H}_{\mathrm{nuc}}\right)}^{\left(1\right)}$].
• Step 2: compute the time-centered expansion term, ${S}^{n+1/2,\star }$.We compute an estimate for the time-centered expansion term in the velocity divergence constraint. Following Bell et al. (2004), we extrapolate to the half-time using S at the previous and current time steps,

  $$\begin{eqnarray}&&{S}^{n+1/2,\star }={S}^{n}+\displaystyle \frac{{\rm{\Delta }}{t}^{n}}{2}\left(\displaystyle \frac{{S}^{n}-{S}^{n-1}}{{\rm{\Delta }}{t}^{n-1}}\right).\end{eqnarray} \tag{ 14 }$$

  Note that in the first time step we average S⁰ and S¹ from the initialization step.
• Step 3: construct a face-centered, time-centered advective velocity, ${{\boldsymbol{U}}}^{\mathrm{ADV},\star }$.The construction of face-centered time-centered states used to discretize the advection terms for velocity, species, and enthalpy, are performed using a standard multidimensional corner transport upwind approach (Colella 1990; Saltzman 1994) with the piecewise-parabolic method one-dimensional tracing (Colella & Woodward 1984). The full details of this Godunov advection approach for all steps in this algorithm are described in Appendix A of Zingale et al. (2015). Here we use Equation (2) to compute face-centered, time-centered velocities, ${{\boldsymbol{U}}}^{\mathrm{ADV},\dagger ,\star }$. The † superscript refers to the fact that the predicted velocity field does not satisfy the divergence constraint,

  $$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left({\beta }_{0}^{n}{{\boldsymbol{U}}}^{\mathrm{ADV},\star }\right)={\beta }_{0}^{n}\left[{S}^{n+1/2* * ,\star }-\displaystyle \frac{1}{{\overline{{\rm{\Gamma }}}}_{1}^{n}{p}_{0}^{n}}{\left(\displaystyle \frac{\partial {p}_{0}}{\partial t}\right)}^{n-1/2}\right].\end{eqnarray} \tag{ 15 }$$

  We project ${{\boldsymbol{U}}}^{\mathrm{ADV},\dagger ,\star }$ onto the space of velocities that satisfy the constraint to obtain ${{\boldsymbol{U}}}^{\mathrm{ADV},\star }$. Each projection step in the algorithm involves the solution of a variable-coefficient Poisson solve using multigrid. Note that we still employ velocity-splitting as described by Equation (10) for this step in order to enforce the appropriate behavior of the system near the edge of the star as determined by the cutoff density. The details of this "MAC" projection are provided in Appendix.
• Step 4: advect the thermodynamic variables over a time interval of ${\rm{\Delta }}t.$

  • A.  
    Update $(\rho {X}_{k})$ using a discretized version of

    $$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {X}_{k})}{\partial t}=-{\rm{\nabla }}\cdot (\rho {X}_{k}{\boldsymbol{U}}),\end{eqnarray} \tag{ 16 }$$

    where the reaction terms have been omitted since they were already accounted for in React State. The update consists of two steps:

    • i.  
      Compute the face-centered, time-centered species, ${(\rho {X}_{k})}^{n+1/2,\mathrm{pred},\star }$, for the conservative update of ${(\rho {X}_{k})}^{(1)}$ using a Godunov approach (Zingale et al. 2015). As described in Paper V, for robust numerical slope limiting we predict ${\rho }^{{\prime} n}={\rho }^{n}-{\rho }_{0}^{n}$ and X_kⁿ to faces and here we spatially interpolate ${\rho }_{0}^{n}$ to faces to assemble the fluxes.
    • ii.  
      Evolve (ρX_k)⁽¹⁾ → (ρX_k)^(2),⋆ using

      $$\begin{eqnarray}\begin{array}{rcl}{\left(\rho {X}_{k}\right)}^{(2),\star } & = & {\left(\rho {X}_{k}\right)}^{\left(1\right)}\\ & & -{\rm{\Delta }}t\left\{{\rm{\nabla }}\cdot \left[{{\boldsymbol{U}}}^{\mathrm{ADV},\star }{\left(\rho {X}_{k}\right)}^{n+1/2,\mathrm{pred},\star }\right]\right\}.\end{array}\end{eqnarray} \tag{ 17 }$$
  • B.  
    Update ρ₀ by calling Average $[{\rho }^{(2),\star }]\to [{\rho }_{0}^{n+1,\star }]$.
  • C.  
    Update p₀ by calling Enforce HSE $[{\rho }_{0}^{n+1,\star }]\,\to \,[{p}_{0}^{n+1,\star }]$.
  • D.  
    Update the enthalpy using a discretized version of equation

    $$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho h)}{\partial t}=-{\rm{\nabla }}\cdot (\rho h{\boldsymbol{U}})+\displaystyle \frac{{{Dp}}_{0}}{{Dt}}+\rho {H}_{\mathrm{nuc}},\end{eqnarray} \tag{ 18 }$$

    again omitting the reaction terms since we already accounted for them in React State. This equation takes the form:

    $$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial (\rho h)}{\partial t} & = & -{\rm{\nabla }}\cdot (\rho h{\boldsymbol{U}})+\displaystyle \frac{\partial {p}_{0}}{\partial t}\\ & & +({\boldsymbol{U}}\cdot {{\boldsymbol{e}}}_{r})\displaystyle \frac{\partial {p}_{0}}{\partial r}.\end{array}\end{eqnarray} \tag{ 19 }$$

    For spherical geometry, we solve the analytically equivalent form,

    $$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial (\rho h)}{\partial t} & = & -{\rm{\nabla }}\cdot (\rho h{\boldsymbol{U}})+\displaystyle \frac{\partial {p}_{0}}{\partial t}\\ & & +{\rm{\nabla }}\cdot ({\boldsymbol{U}}{p}_{0})-{p}_{0}{\rm{\nabla }}\cdot {\boldsymbol{U}}.\end{array}\end{eqnarray} \tag{ 20 }$$

    The update consists of two steps:

    • i.  
      Compute the face-centered, time-centered enthalpy, ${(\rho h)}^{n+1/2,\mathrm{pred},\star },$ for the conservative update of ${(\rho h)}^{(1)}$ using a Godunov approach (Zingale et al. 2015). As described in Paper V, for robust numerical slope limiting we predict ${(\rho h)}^{{\prime} n}\,={(\rho h)}^{n}-(\rho h{)}_{0}^{n}$ to faces, where $(\rho h{)}_{0}^{n}$ is obtained by calling Average $[{(\rho h)}^{n}]\to [(\rho h{)}_{0}^{n}]$, and here we spatially interpolate $(\rho h{)}_{0}^{n}$ to faces to assemble the fluxes.
    • ii.  
      Evolve (ρh)⁽¹⁾ → (ρh)^(2),⋆ using

      $$\begin{eqnarray}\begin{array}{rcl}{\left(\rho h\right)}^{\left(2\right),\star } & = & {\left(\rho h\right)}^{\left(1\right)}-{\rm{\Delta }}t\left\{{\rm{\nabla }}\cdot \left[{{\boldsymbol{U}}}^{\mathrm{ADV},\star }{\left(\rho h\right)}^{n+1/2,\mathrm{pred},\star }\right]\right\}\\ & & +{\rm{\Delta }}t\displaystyle \frac{{{Dp}}_{0}}{{Dt}}\end{array}\end{eqnarray} \tag{ 21 }$$

      where here

      $$\begin{eqnarray}\begin{array}{r}\displaystyle \frac{{Dp}_{0}}{Dt}=\left\{\begin{array}{cc}{\textstyle \tfrac{{p}_{0}^{n+1,\ast }-{p}_{0}^{n}}{{\rm{\Delta }}t}}+\left({{\boldsymbol{U}}}^{{\rm{A}}{\rm{D}}{\rm{V}},\star }\cdot {{\boldsymbol{e}}}_{r}\right){\left({\textstyle \tfrac{{\rm{\partial }}{p}_{0}}{{\rm{\partial }}r}}\right)}^{n} & ({\rm{p}}{\rm{l}}{\rm{a}}{\rm{n}}{\rm{a}}{\rm{r}})\\ {\textstyle \tfrac{{p}_{0}^{n+1,\ast }-{p}_{0}^{n}}{{\rm{\Delta }}t}}+\left[{\rm{\nabla }}\cdot \left({{\boldsymbol{U}}}^{{\rm{A}}{\rm{D}}{\rm{V}},\star }{p}_{0}^{n}\right)-{p}_{0}^{n}{\rm{\nabla }}\cdot {{\boldsymbol{U}}}^{{\rm{A}}{\rm{D}}{\rm{V}},\star }\right] & ({\rm{s}}{\rm{p}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{l}})\end{array}\right.,\end{array}\end{eqnarray} \tag{ 22 }$$

      and ${p}_{0}^{n+1/2}=({p}_{0}^{n}+{p}_{0}^{n+1,* })/2$.
• Step 5: react the thermodynamic variables over the second Δt/2 interval. Call React State [(ρ X_k)^(2),⋆, (ρh)^(2),⋆, ${p}_{0}^{n+1,\star }]$ → [(ρX_k)^n+1,⋆, (ρh)^n+1,⋆, ${\left(\rho {\dot{\omega }}_{k}\right)}^{n+1,\star },{\left(\rho {H}_{\mathrm{nuc}}\right)}^{n+1,\star }].$
• Step 6: compute the time-centered expansion term, ${S}^{n+1/2,\star }$. First, compute ${S}^{n+1,\star }$ with

  $$\begin{eqnarray}&&{S}^{n+1,\star }={\left(-\sigma \displaystyle \sum _{k}{\xi }_{k}{\dot{\omega }}_{k}+\displaystyle \frac{1}{\rho {p}_{\rho }}\displaystyle \sum _{k}{p}_{{X}_{k}}{\dot{\omega }}_{k}+\sigma {H}_{\mathrm{nuc}}\right)}^{n+1,\star }.\end{eqnarray} \tag{ 23 }$$

  Then, define

  $$\begin{eqnarray}&&{S}^{n+1/2}=\displaystyle \frac{{S}^{n}+{S}^{n+1,\star }}{2},\end{eqnarray} \tag{ 24 }$$
• Step 7: construct a face-centered, time-centered advective velocity, ${{\boldsymbol{U}}}^{\mathrm{ADV}}$. The procedure to construct ${{\boldsymbol{U}}}^{\mathrm{ADV},\dagger }$ is identical to the Godunov procedure for computing ${{\boldsymbol{U}}}^{\mathrm{ADV},\dagger ,\star }$ in Step 3, but uses the updated value ${S}^{n+1/2}$ rather than ${S}^{n+1/2,\star }$. The † superscript refers to the fact that the predicted velocity field does not satisfy the divergence constraint,

  $$\begin{eqnarray}&&\begin{array}{l}\quad {\rm{\nabla }}\cdot \left({\beta }_{0}^{n+1/2}{{\boldsymbol{U}}}^{\mathrm{ADV}}\right)={\beta }_{0}^{n+1/2}\\ \quad \times \,\left[{S}^{n+1/2}-\displaystyle \frac{1}{{\overline{{\rm{\Gamma }}}}_{1}^{n+1/2}{p}_{0}^{n+1/2}}{\left(\displaystyle \frac{\partial {p}_{0}}{\partial t}\right)}^{n+1/2}\right],\end{array}\end{eqnarray} \tag{ 25 }$$

  with

  $$\begin{eqnarray}\begin{array}{rcl}{\beta }_{0}^{n+1/2} & = & \displaystyle \frac{{\beta }_{0}^{n}+{\beta }_{0}^{n+1,\star }}{2},\\ {\overline{{\rm{\Gamma }}}}_{1}^{n+1/2} & = & \displaystyle \frac{{\overline{{\rm{\Gamma }}}}_{1}^{n}+{\overline{{\rm{\Gamma }}}}_{1}^{n+1,\star }}{2}.\end{array}\end{eqnarray} \tag{ 26 }$$

  we project ${{\boldsymbol{U}}}^{\mathrm{ADV},\dagger }$ onto the space of velocities that satisfy the constraint to obtain ${{\boldsymbol{U}}}^{\mathrm{ADV}}$ using a MAC projection (see Appendix).
• Step 8: advect the thermodynamic variables over a time interval of ${\rm{\Delta }}t.$

  • A.  
    Update $(\rho {X}_{k})$. This step is identical to Step 4A except we use the updated values ${{\boldsymbol{U}}}^{\mathrm{ADV}}$ and ${\rho }_{0}^{n+1,\star }$ rather than ${{\boldsymbol{U}}}^{\mathrm{ADV},\star }$ and ${\rho }_{0}^{n}$. In particular:

    • i.  
      Compute the face-centered, time-centered species, ${(\rho {X}_{k})}^{n+1/2,\mathrm{pred}}$, for the conservative update of ${(\rho {X}_{k})}^{(1)}$ using a Godunov approach (Zingale et al. 2015). Again, we predict ${\rho }^{{\prime} n}={\rho }^{n}-{\rho }_{0}^{n}$ and X_kⁿ to faces but here we spatially interpolate $({\rho }_{0}^{n}+{\rho }_{0}^{n+1,* })/2$ to faces to assemble the fluxes.
    • ii.  
      Evolve (ρX_k)⁽¹⁾ → (ρX_k)⁽²⁾ using

      $$\begin{eqnarray}\begin{array}{rcl}{\left(\rho {X}_{k}\right)}^{\left(2\right)} & = & {\left(\rho {X}_{k}\right)}^{\left(1\right)}\\ & & -{\rm{\Delta }}t\left\{{\rm{\nabla }}\cdot \left[{{\boldsymbol{U}}}^{\mathrm{ADV}}{\left(\rho {X}_{k}\right)}^{n+1/2,\mathrm{pred}}\right]\right\},\end{array}\end{eqnarray} \tag{ 27 }$$
  • B.  
    Update ρ₀ by calling Average $[{\rho }^{(2)}]\to [{\rho }_{0}^{n+1}]$.
  • C.  
    Update p₀ by calling Enforce HSE $[{\rho }_{0}^{n+1}]\to [{p}_{0}^{n+1}]$.
  • D.  
    Update the enthalpy. This step is identical to Step 4D except we use the updated values ${{\boldsymbol{U}}}^{\mathrm{ADV}}$, ${\rho }_{0}^{n+1}$, $(\rho h{)}_{0}^{n+1}$, and ${p}_{0}^{n+1}$ rather than ${{\boldsymbol{U}}}^{\mathrm{ADV},\star },{\rho }_{0}^{n}$, $(\rho h{)}_{0}^{n}$, and p₀ⁿ. In particular:

    • i.  
      Compute the face-centered, time-centered enthalpy, ${(\rho h)}^{n+1/2,\mathrm{pred}},$ for the conservative update of ${(\rho h)}^{(1)}$ using a Godunov approach (Zingale et al. 2015). Again, we predict ${(\rho h)}^{{\prime} n}={(\rho h)}^{n}-(\rho h{)}_{0}^{n}$ to faces but here we spatially interpolate $[(\rho h{)}_{0}^{n})+(\rho h{)}_{0}^{n+1,* }]/2$ to faces to assemble the fluxes.
    • ii.  
      Evolve (ρh)⁽¹⁾ → (ρh)⁽²⁾.

      $$\begin{eqnarray}\begin{array}{rcl}{\left(\rho h\right)}^{\left(2\right)} & = & {\left(\rho h\right)}^{\left(1\right)}-{\rm{\Delta }}t\left\{{\rm{\nabla }}\cdot \left[{{\boldsymbol{U}}}^{\mathrm{ADV}}{\left(\rho h\right)}^{n+1/2,\mathrm{pred}}\right]\right\}\\ & & +{\rm{\Delta }}t\displaystyle \frac{{{Dp}}_{0}}{{Dt}}\end{array}\end{eqnarray} \tag{ 28 }$$

      where here

      $$\begin{eqnarray}\begin{array}{r}\displaystyle \frac{{Dp}_{0}}{Dt}=\left\{\begin{array}{lc}{\textstyle \tfrac{{p}_{0}^{n+1}-{p}_{0}^{n}}{{\rm{\Delta }}t}}+\left({{\boldsymbol{U}}}^{{\rm{A}}{\rm{D}}{\rm{V}}}\cdot {{\boldsymbol{e}}}_{r}\right){\left({\textstyle \tfrac{{\rm{\partial }}{p}_{0}}{{\rm{\partial }}r}}\right)}^{n+1/2} & ({\rm{p}}{\rm{l}}{\rm{a}}{\rm{n}}{\rm{a}}{\rm{r}})\\ {\textstyle \tfrac{{p}_{0}^{n+1}-{p}_{0}^{n}}{{\rm{\Delta }}t}}+\left[{\rm{\nabla }}\cdot \left({{\boldsymbol{U}}}^{{\rm{A}}{\rm{D}}{\rm{V}}}{p}_{0}^{n+1/2}\right)-{p}_{0}^{n+1/2}{\rm{\nabla }}\cdot {{\boldsymbol{U}}}^{{\rm{A}}{\rm{D}}{\rm{V}}}\right] & ({\rm{s}}{\rm{p}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{l}})\end{array}\right.,\end{array}\end{eqnarray} \tag{ 29 }$$

      and ${p}_{0}^{n+1/2}=({p}_{0}^{n}+{p}_{0}^{n+1})/2$.
• Step 9: react the thermodynamic variables over the second Δt/2 interval. Call React State [(ρX_k)⁽²⁾, (ρh)⁽²⁾, ${p}_{0}^{n+1}$] → [(ρX_k)ⁿ⁺¹, (ρh)ⁿ⁺¹, ${\left(\rho {\dot{\omega }}_{k}\right)}^{n+1},{\left(\rho {H}_{\mathrm{nuc}}\right)}^{n+1}$].
• Step 10: define the new-time expansion term, Sⁿ⁺¹.

  • A.  
    Define

    $$\begin{eqnarray}&&{S}^{n+1}={\left(-\sigma \displaystyle \sum _{k}{\xi }_{k}{\dot{\omega }}_{k}+\sigma {H}_{\mathrm{nuc}}+\displaystyle \frac{1}{\rho {p}_{\rho }}\displaystyle \sum _{k}{p}_{{X}_{k}}{\dot{\omega }}_{k}\right)}^{n+1}.\end{eqnarray} \tag{ 30 }$$
• Step 11: update the velocity. First, we compute the face-centered, time-centered velocities, ${{\boldsymbol{U}}}^{n+1/2,\mathrm{pred}}$ using a Godunov approach (Zingale et al. 2015). Then, we update the velocity field ${{\boldsymbol{U}}}^{n}$ to ${{\boldsymbol{U}}}^{n+1,\dagger }$ by discretizing Equation (2) as

  $$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{U}}}^{n+1,\dagger }={{\boldsymbol{U}}}^{n}-{\rm{\Delta }}t\left[{{\boldsymbol{U}}}^{\mathrm{ADV}}\cdot {\rm{\nabla }}{{\boldsymbol{U}}}^{n+1/2,\mathrm{pred}}\right]\\ \quad -\,{\rm{\Delta }}t\left[\displaystyle \frac{{\beta }_{0}^{n+1/2}}{{\rho }^{n+1/2}}{\rm{\nabla }}\left(\displaystyle \frac{{\pi }^{n-1/2}}{{\beta }_{0}^{n-1/2}}\right)\right.\\ \quad \left.+\,\displaystyle \frac{\left({\rho }^{n+1/2}-{\rho }_{0}^{n+1/2}\right)}{{\rho }^{n+1/2}}{g}^{n+1/2}{{\boldsymbol{e}}}_{r}\right],\end{array}\end{eqnarray} \tag{ 31 }$$

  where

  $$\begin{eqnarray}\begin{array}{rcl}{\rho }^{n+1/2} & = & \displaystyle \frac{{\rho }^{n}+{\rho }^{n+1}}{2},\\ {\rho }_{0}^{n+1/2} & = & \displaystyle \frac{{\rho }_{0}^{n}+{\rho }_{0}^{n+1}}{2}.\end{array}\end{eqnarray} \tag{ 32 }$$

  Again, the † superscript refers to the fact that the updated velocity does not satisfy the divergence constraint,

  $$\begin{eqnarray}&&\begin{array}{l}{\rm{\nabla }}\cdot \left({\beta }_{0}^{n+1}{{\boldsymbol{U}}}^{n+1}\right)\\ =\,{\beta }_{0}^{n+1}\left[{S}^{n+1}-\displaystyle \frac{1}{{\overline{{\rm{\Gamma }}}}_{1}^{n+1}{p}_{0}^{n+1}}{\left(\displaystyle \frac{\partial {p}_{0}}{\partial t}\right)}^{n+1/2}\right].\end{array}\end{eqnarray} \tag{ 33 }$$

  We use an approximate projection to project ${{\boldsymbol{U}}}^{n+1,\dagger }$ onto the space of velocities that satisfy the constraint to obtain ${{\boldsymbol{U}}}^{n+1}$ using a "nodal" projection. This projection necessarily differs from the MAC projection used in Step 3 and Step 7 because the velocities in those steps are defined on faces and ${{\boldsymbol{U}}}^{n+1}$ is defined at cell centers, requiring different divergence and gradient operators. Furthermore, as part of the nodal projection, we also define a nodal new-time perturbational pressure, ${\pi }^{n+1/2}$. Refer to Appendix for more details.

This completes one step of the algorithm.

To initialize the simulation we use the same procedure described in Paper III. At the beginning of each simulation, we define $({\boldsymbol{U}},\rho {X}_{k},\rho h)$. We set initial values for ${\boldsymbol{U}},\rho {X}_{k}$, and $\rho h$ and perform a sequence of projections (to ensure the velocity field satisfies the divergence constraint) followed by a small number of steps of the temporal advancement scheme to iteratively find initial values for ${\pi }^{n-1/2}$ and S⁰ and S¹ for use in the first time step.

Our approach to AMR is algorithmically the same as the treatment described in Section 5 of Paper V; we refer the reader there for details. MAESTROeX supports refinement ratios of 2 between levels. We note that for spherical problems, AMR is only available for the case of a uniformly spaced base state.

We perform weak scaling tests for simulations of convection preceding ignition in a spherical, full-star Chandrasekhar mass white dwarf. The simulation setup remains the same as reported in Section 3 of Nonaka et al. (2011) and originally used in Zingale et al. (2011), and thus we emphasize that these scaling tests are performed using meaningful, scientific calculations. Here, we perform simulations using ${256}^{3},{512}^{3},{768}^{3},{1024}^{3},{1280}^{3}$, and 1536³ grid cells on a spatially uniform grid (no AMR). We divide each simulation into 64³ grids, so these simulations contain between 64 grids (256³) and 13,824 grids (1536³). These simulations were performed using the NERSC cori system on the Intel Xeon Phi (KNL) partition. Each node contains 68 cores, each capable of supporting up to four hardware threads (i.e., a maximum of 272 hardware threads per node). For these tests, we assign four MPI tasks to each node, and 16 OpenMP threads per MPI process. Each MPI task is assigned to a single grid, so our tests use between 64 and 13,824 MPI processes (i.e., between 1024 and 221,184 total OpenMP threads). For 64³ grids we discovered that using more than 16 OpenMP threads did not decrease the wallclock time due to a lack of work available per grid; in principle, one could use larger grids, fewer MPI processes, and more threads per MPI process to obtain a flatter weak scaling curve; however, the overall wallclock time would increase except for extremely large numbers of MPI processes (beyond the range we tested here). Thus, the more accurate measure of weak scaling is to consider the number of MPI processes, since the scaling plot would look virtually identical for larger thread counts. Note that the largest simulation used roughly 36% of the entire computational system.

In the left panel of Figure 3 we compare the wallclock time per time step as a function of total core count (in this case, the total number of OpenMP threads) for the original FBoxLib-based MAESTRO implementation to the AMReX MAESTROeX implementation. These tests were performed using the original temporal integration strategy in Nonaka et al. (2010), noting that the new temporal integration with and without the irregular base state gives essentially the same results. We also include a plot of MAESTROeX without base state evolution. Comparing the original and new implementations, we see similar scaling results except for the largest simulation, where MAESTROeX performs better. We see that the increase in wallclock time from the smallest to largest simulations is roughly 42%. We also note that without base state evolution, the code runs 14% faster for small problems, and scales much better with wallclock time from the smallest to largest simulation increasing by only 13%. This is quite remarkable since there are three linear solves per time step (two cell-centered Poisson solves used in the MAC projection, and a nodal Poisson solve used to compute the updated cell-centered velocities). Contrary to our prior assumptions, the linear solves are not the primary scaling bottleneck in this code. In the right panel of Figure 3, we isolate the wallclock time required for these linear solves and see that (i) the linear solves only use 20%–23% of the total computational time, and (ii) the increase in the solver wallclock time from the smallest to largest simulations is only 28%. Further profiling reveals that the primary scaling bottleneck is the average operator. The averaging operator requires collecting the sum of Cartesian data onto one-dimensional arrays holding every possible mapping radius. This amounts to at least 24,384 double precision values (for the 256³ simulation) up to 883,584 values (for the 1536³ simulation). The averaging operator requires a global sum reduction over all processors, and the communication of this data is the primary scaling bottleneck. For the simulation with base state evolution, this averaging operator is only called once per time step (as opposed to 14 times per time step when base state evolution is included). The difference in total wallclock times with and without base state evolution is almost entirely due to the averaging. Note that as expected, advection, reactions, and calls to the EOS scale almost perfectly, since there is only a single parallel communication call to fill ghost cells.

Zoom In Zoom Out Reset image size

To explore the accuracy of the new temporal algorithm, we now analyze in detail three-dimensional, full-star calculations of convection preceding ignition in a white dwarf. Again, we refer the reader to Nonaka et al. (2011) and Zingale et al. (2011) for setup details. We implement both uniformly and irregularly spaced base state with the new temporal algorithm, while only uniform base state spacing is used in the original algorithm. As in Section 3 of Nonaka et al. (2011), we choose the peak temperature and peak Mach number as the two diagnostics to compare the simulations. Figure 4 shows the evolution of both peak temperature and peak Mach number until time of ignition on a single-level grid with resolution of 256³. The simulation using the new temporal scheme with uniformly spaced base state gives the same qualitative results as the original scheme, and predicts a similar time of ignition (t = 7810 s compared to t = 7850 s for original algorithm). The simulation using the new temporal scheme with irregularly spaced base state displays a slightly different peak temperature behavior during the initial transition period $t\lt 150\,{\rm{s}}$, which results in the difference between the curves post transition. We strongly suspect that this is a result of using a different initial model file (the resolution near the center of the star is much coarser with the irregular spacing than the uniform spacing). Fortunately, the simulation with irregular base state spacing still follows the same trend as with uniform spacing, and the star is shown to ignite at an earlier time t = 6840 s.

Zoom In Zoom Out Reset image size

Figure 5 shows the peak temperature evolution over the first 1000 s on two grids of differing resolutions, 256³ and 512³. Limited allocations prevented us from running this simulation further. As previously suspected, the simulation using irregularly spaced base state agrees much closer with the results from using uniform spacing as the resolution increases. This is most likely due to the increased resolution of the initial model, which more closely matches the uniformly spaced counterpart. This is especially important when computing the base state pressure from base state density, which is particularly sensitive to coarse resolution near the center of the star.

Zoom In Zoom Out Reset image size

In terms of efficiency, all three simulations on the 256³ single-level grid were run on Cori haswell with 64 processors and 8 threads per core and their run times were compared. As a result of simplifying the algorithm by eliminating the evolution equations for the base state density and pressure, the simulation using the new temporal algorithm took only 6.75 s per time step with uniformly spaced base state, which is 13% faster than the 7.77 s per time step when using the original scheme. However, we do observe a 25% increase in run time of 9.72 s when using irregularly spaced base state with the new algorithm. This can be explained by the irregularly spaced base state array being much larger in size than its uniformly spaced counterpart, and thus require additional communication and computation time. One possible strategy to significantly reduce the run time is to consider truncating the base state beyond the cutoff density radius.

We now test the performance of MAESTROeX for adaptive, three-dimensional simulations to track localized regions of interest over time. Figure 6 illustrates the initial grid structures with two levels of refinement for both planar and spherical geometries where the grid is refined according to the temperature and density profiles, respectively. For each of the problems we tested, the single-level simulation was run using the original temporal scheme and the adaptive simulations using the original and new temporal algorithms. We want to show that the adaptive simulation can give similar results to the single-level simulation and in a more computationally efficient manner.

Zoom In Zoom Out Reset image size

In the planar case, we use the same problem setup for a hot bubble rising in a white dwarf environment as described in Section 6 of Paper V. Here we use a domain size of 3.6 × 10⁷ cm by 3.6 × 10⁷ cm by 2.88 × 10⁸ cm, and allow the grid structure to change with time. The single-level simulation at a resolution of 128² × 1024 was run on Cori haswell with 48 processors and took approximately 33.5 s per time step (averaged over 10 time steps) using either the original or new temporal algorithm. The adaptive simulation has a resolution of 32² × 256 at the coarsest level, resulting in the same effective resolution at the finest level as the single-level simulation. We tag cells that satisfy $T\mbox{--}\bar{T}\gt 3\times {10}^{7}$ K as well as all cells at that height. The adaptive run took only 3.7 s per time step, and this 89% decrease in run time is mostly due to the fact that initially only 6.25% of the cells (1,048,576 out of 128² × 1024 cells) are refined at the finest level. Figure 7 shows a series of planar slices of the temperature profile at time intervals of 1.25 s, and verifies that the adaptive simulation is able to capture the same dynamics as the single-level simulation at much lower computational cost.

Zoom In Zoom Out Reset image size

We continue to use the full-star white dwarf problem described in Section 4.2 to test adaptive simulations on spherical geometry. The adaptive grid is refined twice by tagging the density at ρ > 10⁵ g cm⁻³ on the first level and ρ > 10⁸ g cm⁻³ on the second level. These tagging values have been shown to work well previously in Paper V, but we have found that the code may encounter numerical difficulties when the tagging values are too close to each other in subsequent levels of refinement. The simulation on a single-level grid of 512³ resolution took 12.7 s per time step (again averaged over 10 time steps). The adaptive grid has a resolution of 128³ at the coarsest level and an effective resolution of 512³ at the finest level. On this grid, 27.8% of the cells (4,666,880 out of 256³ cells) are refined at the first level and 5.3% (7,176,192 out of 512³ cells) at the second. The adaptive simulation took 5.61 s, resulting in more than a factor of two in speedup. Both simulations are computed to t = 2000 s and we choose to use the peak temperature as the diagnostic to compare the results. Figure 8 shows the evolution of the peak temperature for all three runs and shows that the adaptive simulation gives the same qualitative result as the single-level simulation. We do not expect the curves to match up exactly because the governing equations are highly nonlinear, and slight differences in the solution caused by solver tolerance and discretization error can change the details of the results. Each simulation was run on Cori haswell with 512 processors and 4 threads per core.

Zoom In Zoom Out Reset image size