CASTRO: A NEW COMPRESSIBLE ASTROPHYSICAL SOLVER. I. HYDRODYNAMICS AND SELF-GRAVITY

At the beginning of each time step, we assume that, in the case of self-gravity, g is defined consistently with the current mass distribution in U. The algorithm at a single level of refinement is composed of the following steps:

Step 1. Advance the nuclear reaction network through a time interval of Δt/2.

Define U⁽¹⁾ = Uⁿ with the exception of

$$\begin{eqnarray} (\rho E)^{(1)} &=& (\rho E)^n + \frac{{\Delta t}}{2}(\rho {H_{\rm nuc}})^n, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} (\rho X_k)^{(1)} &=& (\rho X_k)^n + \frac{{\Delta t}}{2}(\rho {\dot{\omega }}_k)^n, \end{eqnarray} \tag{ 14 }$$

where (ρH_nuc)ⁿ and $(\rho {\dot{\omega }}_k)^n$ are computed using calls to the user-defined reaction network. Note that ρ is unchanged during this step.

Step 2. Advect the solution through Δt.

Advance the solution using time-centered fluxes and an explicit representation of the source term, neglecting the contribution from reactions which are taken into account in steps 1 and 4 (the asterisk superscript notation indicates that we will later correct this state to effectively time center the source terms):

$$\begin{equation} {\bf U}^{(2,\ast \ast)} = {\bf U}^{(1)} - {\Delta t}\nabla \cdot {\bf F}^{n+{1/2}}+ {\Delta t}{\bf S}^{(1)}, \end{equation} \tag{ 15 }$$

where

$$\begin{equation} {\bf S}_{{\bf U}}^{(1)} = \left(\begin{array}{@{}c@{}} S_{\rho } \\ {\bf S}_{\rho {\bf u}} \\ S_{\rho E} \\ S_{\rho X_k} \\ S_{\rho C^{\rm adv}_k} \\ S_{\rho C^{\rm aux}_k} \end{array}\right)^{(1)} = \left(\begin{array}{@{}c@{}} S_{{\rm ext},\rho } \\ (\rho {\bf g})^{(1)} + {\bf S}_{{\rm ext},\rho {\bf u}} \\ (\rho {\bf u}\cdot {\bf g})^{(1)} + S_{{\rm ext},\rho E} \\ S_{{\rm ext},\rho X_k} \\ S_{{\rm ext},\rho C^{\rm adv}_k} \\ S_{{\rm ext},\rho C^{\rm aux}_k} \end{array}\right)^{(1)}. \end{equation} \tag{ 16 }$$

The construction of the fluxes is described in detail in Section 5.2. We note that in the single-level algorithm we can use the gravitational forcing computed in step 3 of the previous time step, since the density has not changed.

After the advective update, we ensure that the solution is physically meaningful by forcing the density to exceed a non-negative, user-defined minimum value. We also ensure that the mass fractions are all non-negative and sum to one.

We also have an option for a user-defined sponge in order to prevent the velocities in the upper atmosphere from becoming too large, and subsequently, the time step from becoming too small. We multiply the velocity by 1/(1 + Δt κ f_damp(ρ)), where κ is the sponge strength, and f_damp is a smooth function of density that varies from 0 to 1. Full details of the sponge are given in Zingale et al. (2009). Finally, we adjust (ρE)^(2,**) to be consistent with u^(2,**).

Step 3. Correct the solution with time-centered source terms and compute gravity at tⁿ⁺¹.

We correct the solution by effectively time-centering the source terms. First, we correct U with updated external sources:

$$\begin{equation} {\bf U}^{(2),\ast } = {\bf U}^{(2,\ast \ast)} + \frac{{\Delta t}}{2}\big(S_{{\rm ext},{\bf U}}^{(2,\ast \ast)} - S_{{\rm ext},{\bf U}}^{(1)}\big). \end{equation} \tag{ 17 }$$

Next, we evaluate gravity using ρ^(2,*). If using full gravity we solve

$$\begin{equation} {\bf g}^{(2,\ast)} = -\nabla \phi ^{(2, \ast)}, \qquad \nabla ^2 \phi ^{(2, \ast)} = 4\pi G\rho ^{(2, \ast)}, \end{equation} \tag{ 18 }$$

where we supply an initial guess for ϕ^(2,*) from the previous solve. In the single-level algorithm described here, g^(2,*) is saved to be used as g⁽¹⁾ in step 2 of the next time step. This suffices in the single-level algorithm because ρ does not change between the end of step 3 of one time step and the start of step 2 of the next time step.

We then correct the solution with the updated gravity:

$$\begin{eqnarray} (\rho {\bf u})^{(2)} &=& (\rho {\bf u})^{(2,\ast)} + \frac{{\Delta t}}{2}[(\rho {\bf g})^{(2,\ast)} - (\rho {\bf g})^{(1)}], \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} (\rho E)^{(2)} &=& (\rho E)^{(2,\ast)} + \frac{{\Delta t}}{2} [(\rho {\bf u}\, {\bm \cdot}\, {\bf g})^{(2,\ast)} - (\rho {\bf u}\, {\bm \cdot}\, {\bf g})^{(1)}]. \end{eqnarray} \tag{ 20 }$$

For all other conserved variables other than ρu and ρE, U⁽²⁾ = U^(2,*). We note here that the time discretization of the gravitational forcing terms differs from that in the FLASH (Fryxell et al. 2000) and ENZO (O'Shea et al. 2005) codes, where the gravitational forcing at t^n+1/2 is computed by extrapolation from values at tⁿ and tⁿ⁻¹ (see also Bryan et al. 1995). Our discretization of the gravitational terms is consistent with our predictor-corrector approach in the handling of other source terms.

Step 4. Advance the nuclear reaction network through a time interval of Δt/2.

Define Uⁿ⁺¹ = U⁽²⁾ with the exception of

$$\begin{eqnarray} (\rho E)^{n+1} &=& (\rho E)^{(2)} + \frac{{\Delta t}}{2}(\rho {H_{\rm nuc}})^{(2)}, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} (\rho X_k)^{n+1} &=& (\rho X_k)^{(2)} + \frac{{\Delta t}}{2}(\rho {\dot{\omega }}_k)^{(2)}. \end{eqnarray} \tag{ 22 }$$

We also include an option to modify any component of the new-time state as needed to account for special user requirements.

Step 5. Compute the new time step.

The time step is computed using the standard CFL condition for explicit methods, with additional constraints (such as one based on rate of burning) possible as needed. The user sets a CFL factor, σ^CFL, between 0 and 1. The sound speed, c, is computed by the EOS, and for a calculation in n_dim dimensions,

$$\begin{equation} {\Delta t}= \sigma ^\mathrm{CFL} \ \ \min _{i=1\ldots n_\mathrm{dim}} \ \ \left\lbrace {\Delta t}_i \right\rbrace, \end{equation} \tag{ 23 }$$

where

$$\begin{equation} {\Delta t}_i = \min _{\bf x} \left\lbrace \frac{\Delta x_i}{|{\bf {u}}_i| + c \; } \right\rbrace. \end{equation} \tag{ 24 }$$

min_x is the minimum taken over all computational grid cells in the domain.

This concludes the single-level algorithm description. We note that whenever the kinetic energy dominates the total energy, making the calculation of e from E numerically unreliable, we use a method similar to the dual-energy approach described in Bryan et al. (1995) to compute the internal energy with sufficient precision. In practice, this involves evolving ρe in time and using this solution when appropriate.

We use an unsplit Godunov method with characteristic tracing and full corner coupling in three dimensions (Miller & Colella 2002) to compute time-centered edge states. We have replaced the PPM limiters in Miller & Colella (2002) with an updated PPM algorithm that is designed to preserve accuracy at smooth extrema and is insensitive to asymmetries caused by roundoff error (Colella & Sekora 2008; McCorquodale & Colella 2010). CASTRO also has options to use the unsplit piecewise-linear algorithm described in Colella (1990), Saltzman (1994), or to retain the PPM limiters in Miller & Colella (2002), which were originally developed in Colella & Woodward (1984) using a split integrator.

There are four major steps in the construction of the face-centered fluxes, F^n+1/2, that are used in step 2 in Section 5.1 to update the solution. We also include details on the solution of the Riemann problem. In summary,

Step 2.1. Rewrite the state, U⁽¹⁾, in terms of primitive variables, Q⁽¹⁾.

Step 2.2. Construct a piecewise parabolic approximation of Q⁽¹⁾ within each cell.

Step 2.3. Predict average values of Q⁽¹⁾ on edges over the time step using characteristic extrapolation.

Step 2.4. Compute fluxes, F^n+1/2, using an approximate Riemann problem solver.

We expand each of these steps in more detail below.

Step 2.1. Compute primitive variables and source terms.

We define Q⁽¹⁾ = (ρ, u, p, ρe, X_k, C^adv_k, C^aux_k)⁽¹⁾. The pressure, p, is computed through a call to the EOS using ρ, e, and X_k. Note that we also include ρe in Q; this quantity is used in the approximate Riemann solver to avoid an EOS call to evaluate the energy flux, analogous to the effective dynamics for γ = p/(ρe) + 1 in the Colella & Glaz (1985) approximate Riemann solver.

For the overall integration algorithm, we want to include the effect of source terms except for reactions in the characteristic tracing (step 2.2). (Reactions are treated separately using a symmetric operator-split approach in steps 1 and 4 of the algorithm.) The time evolution equations written in terms of the primitive variables, Q, and omitting contributions from reactions, are

$$\begin{equation*} \frac{\partial {\bf Q}}{\partial t} = \left(\begin{array}{@{}c@{}} -{\bf u}\cdot \nabla \rho - \rho \nabla \cdot {\bf u}\\ -{\bf u}\cdot \nabla {\bf u}- \frac{1}{\rho }\nabla p \\ -{\bf u}\cdot \nabla p - \rho c^2\nabla \cdot {\bf u}\nonumber \\ - {\bf u}\cdot \nabla (\rho e) - (\rho e+p)\nabla \cdot {\bf u}\\ -{\bf u}\cdot \nabla X_k \\ -{\bf u}\cdot \nabla C^{\rm adv}_k \\ -{\bf u}\cdot \nabla C^{\rm aux}_k \end{array}\right) + {\bf S}_{{\bf Q}}, \end{equation*}$$

where

$$\begin{equation} {\bf S}_{{\bf Q}} = \left(\begin{array}{@{}c@{}} S_\rho \\ {\bf S}_{{\bf u}} \\ S_p \\ S_{\rho e} \\ S_{X_k} \\ S_{C^{\rm adv}_k} \\ S_{C^{\rm aux}_k} \end{array}\right) = \left(\begin{array}{@{}c@{}} S_{{\rm ext},\rho } \\ {\bf g}+ \frac{1}{\rho }{\bf S}_{{\rm ext},\rho {\bf u}} \\ \frac{p_e}{\rho } S_{{\rm ext},\rho E} + p_\rho S_{{\rm ext}\rho } + \frac{p_{X_k}}{\rho } S_{{\rm ext},\rho X_k}\\ \ S_{{\rm ext},\rho E} \\ \frac{1}{\rho }S_{{\rm ext},\rho X_k} \\ \frac{1}{\rho }S_{{\rm ext},\rho C^{\rm adv}_k} \\ \frac{1}{\rho }S_{{\rm ext},\rho C^{\rm aux}_k} \end{array}\right). \end{equation} \tag{ 25 }$$

Here, c is the sound speed, defined as $c = \sqrt{\Gamma _1 p /\rho }$, with Γ₁ = dlog p/dlog ρ|_s, with s the entropy. The remaining thermodynamic derivatives are $p_e = \partial p / \partial e |_{\rho, X_k}$, $p_\rho = \partial p / \partial \rho |_{e, X_k}$, and $p_{X_k} = \partial p / \partial X_k |_{\rho, e, X_{j,(j\ne k)}}$. Often, the EOS is a function of ρ, T, and X_k, and returns derivatives with these quantities held constant. In terms of the latter derivatives, our required thermodynamic derivatives are

$$\begin{eqnarray} p_e &\,{=}\,& \left(\left. \frac{\partial e}{\partial T} \right|_{\rho, X_k} \right)^{-1} \left. \frac{\partial p}{\partial T} \right|_{\rho,X_k}, \end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray} p_\rho &\,{=}\,& \left. \frac{\partial p}{\partial \rho } \right|_{T,X_k} - \left(\left. \frac{\partial e}{\partial T} \right|_{\rho, X_k} \right)^{-1} \left. \frac{\partial p}{\partial T} \right|_{\rho,X_k} \left. \frac{\partial e}{\partial \rho } \right|_{T,X_k}\!, \quad \end{eqnarray} \tag{ 26b }$$

$$\begin{eqnarray} p_{X_k} &\,{=}\,& \left. \frac{\partial p}{\partial X_k} \right|_{\rho,T,X_{j,(j\ne k)}} \,{-}\, \left(\!\left. \frac{\partial e}{\partial T} \right|_{\rho, X_k} \!\right)^{-1}\!\! \left. \frac{\partial p}{\partial T} \right|_{\rho,X_k} \!\! \left. \frac{\partial e}{\partial X_k} \right|_{\rho,T,X_{j,(j \ne k)}}\!.\nonumber\\ \end{eqnarray} \tag{ 26c }$$

Step 2.2. Reconstruct parabolic profiles within each cell.

In this step, we construct a limited piecewise parabolic profile of each q in Q (we use q to denote an arbitrary primitive variable from Q). These constructions are performed in each coordinate direction separately. The default option in CASTRO is to use a new limiting procedure that avoids reducing the order of the reconstruction at smooth local extrema. The details of this construction are given in Colella & Sekora (2008) and McCorquodale & Colella (2010). In summary,

• 1.  
  Step 2.2a. For each cell, we compute the spatial interpolation of qⁿ to the high and low faces of cell q_i using a limited cubic interpolation formula. These interpolants are denoted by q_i,+ and q_i,−.
• 2.  
  Step 2.2b. Construct quadratic profiles using q_i,−, q_i, and q_i,+.

  $$\begin{equation} q_i^{\rm quad}(x) = q_{i,-} + \xi (x)\lbrace q_{i,+} - q_{i,-} + q_{6,i}[1-\xi (x)]\rbrace, \end{equation} \tag{ 27 }$$

  $$\begin{equation} q_6 = 6q_{i} - 3(q_{i,-}+q_{i,+}), \end{equation} \tag{ 28 }$$

  $$\begin{equation} \xi (x) = \frac{x - ih}{h}, \quad 0 \le \xi (x) \le 1 \;, \end{equation} \tag{ 29 }$$

  where h is the mesh spacing in the direction of interpolation. Also, as in Miller & Colella (2002), we compute a flattening coefficient, χ ∈ [0, 1], used in the edge state prediction to further limit slopes near strong shocks. The computation of χ is identical to the approach used in FLASH (Fryxell et al. 2000), except that a flattening coefficient of 1 indicates that no additional limiting takes place, whereas a flattening coefficient of 0 means we effectively drop order to a first-order Godunov scheme, which is opposite of the convention used in FLASH.

Step 2.3. Characteristic extrapolation.

We begin by extrapolating Q⁽¹⁾ to edges at t^n+1/2. The edge states are dual valued, i.e., at each face, there is a left state and a right state estimate, denoted q_L,i+1/2 and q_R,i+1/2 (we write the equations in one dimension for simplicity). The spatial extrapolation is one dimensional, i.e., transverse derivatives are omitted and accounted for later.

• 1.  
  Step 2.3a. Integrate the quadratic profiles. We are essentially computing the average value swept out by the quadratic profile across the face assuming the profile is moving at a speed λ_k, where λ_k is a standard wave speed associated with gas dynamics.Define the following integrals, where σ_k = |λ_k|Δt/h:

  $$\begin{eqnarray} \mathcal {I}_{i,+}(\sigma _k) &=& \frac{1}{\sigma _k h}\int _{(i+{1/2})h-\sigma _k h}^{(i+{1/2})h}q_i^{\rm quad}(x)dx, \end{eqnarray} \tag{ 30a }$$

  $$\begin{eqnarray} \mathcal {I}_{i,-}(\sigma _k) &=& \frac{1}{\sigma _k h}\int _{(i-{1/2})h}^{(i-{1/2})h+\sigma _k h}q_i^{\rm quad}(x)dx. \end{eqnarray} \tag{ 30b }$$

  Substituting Equation (27) gives

  $$\begin{eqnarray} \mathcal {I}_{i,+}(\sigma _k) &=& q_{i,+} - \frac{\sigma _k}{2}\left[q_{i,+}-q_{i,-}-\left(1-\frac{2}{3}\sigma _k\right)q_{6,i}\right],\nonumber\\ \end{eqnarray} \tag{ 31a }$$

  $$\begin{eqnarray} \mathcal {I}_{i,-}(\sigma _k) &=& q_{i,-} + \frac{\sigma _k}{2}\left[q_{i,+}-q_{i,-}+\left(1-\frac{2}{3}\sigma _k\right)q_{6,i}\right].\nonumber\\ \end{eqnarray} \tag{ 31b }$$
• 2.  
  Step 2.3b. Obtain a left and right edge state at t^n+1/2 by applying a characteristic tracing operator (with flattening) to the integrated quadratic profiles. Note that we also include the explicit source term contribution:

  $$\begin{eqnarray} q_{L,i+{1/2}} &=& q_i - \chi _i\sum _{k:\lambda _k \ge 0}{\bf l}_k\cdot [q_i-\mathcal {I}_{i,+}(\sigma _k)]{\bf r}_k + \frac{{\Delta t}}{2}S_{q,i}^n,\nonumber\\ \end{eqnarray} \tag{ 32a }$$

  $$\begin{eqnarray} q_{R,i-{1/2}} &=& q_i - \chi _i\sum _{k:\lambda _k \le 0}{\bf l}_k\cdot [q_i-\mathcal {I}_{i,-}(\sigma _k)]{\bf r}_k + \frac{{\Delta t}}{2}S_{q,i}^n.\nonumber\\ \end{eqnarray} \tag{ 32b }$$

  In non-Cartesian coordinates, volume source terms are added to the traced states. Here, r_k and l_k are the standard right column and left row eigenvectors associated with the equations of gas dynamics (see Toro 1997).An unsplit approximation that includes full corner coupling is constructed by constructing increasingly accurate approximations to the transverse derivatives. The details follow exactly as given in Section 4.2.1 in Miller & Colella (2002), except for the solution of the Riemann problem, which is described in step 2.4.

Step 2.4. Compute fluxes

The fluxes are computed using an approximate Riemann solver. The solver used here is essentially the same as that used in P. Colella et al. (1997, unpublished), which is based on ideas discussed in Bell et al. (1989). This solver is computationally faster and considerably simpler than the approximate Riemann solver introduced by Colella & Glaz (1985). The Colella and Glaz solver was based on an effective dynamics for γ and was designed for real gases that are well approximated by this type of model. The approximate Riemann solver used in CASTRO is suitable for a more general convex EOS.

As with other approximate Riemann solvers, an important design principle is to avoid additional evaluations of the EOS when constructing the numerical flux. For that reason, we include ρe in Q and compute (ρe)_L,R. The information carried in ρe is overspecified but it allows us to compute an energy flux without an inverse call to the EOS.

The numerical flux computation is based on approximating the solution to the Riemann problem and evaluating the flux along the x/t = 0 ray. The procedure is basically a two-step process in which we first approximate the solution in phase space and then interpret the phase space solution in real space.

• 1.  
  Step 2.4a. To compute the phase space solution, we first solve for p* and u*, the pressure between the two acoustic waves and the velocity of the contact discontinuity, respectively. These quantities are computed using a linearized approximation to the Rankine–Hugoniot relations. We first define Γ_1,L/R by using the cell-centered values on either side of the interface. Next, we compute Lagrangian sound speeds, $W_L = \sqrt{ \Gamma _{1,L} p_L \rho _L}$ and $W_R = \sqrt{ \Gamma _{1,R} p_R \rho _R}$, and the corresponding Eulerian sound speeds $c_{L,R} = \sqrt{ \Gamma _{1,L,R} p_{L,R} / \rho _{L,R}}$. Then,

  $$\begin{equation} p^* = \frac{W_L p_R + W_R p_L + W_L W_R(u_L-u_R)}{W_L + W_R}, \end{equation} \tag{ 33a }$$

  $$\begin{equation} u^* = \frac{W_L u_L + W_R u_R+ (p_L - p_R)}{W_L + W_R}. \end{equation} \tag{ 33b }$$

  From u* and p*, we can compute

  $$\begin{eqnarray} \rho ^*_{L,R} &=& \rho _{L,R} + \frac{p^* - p_{L,R}}{c^2_{L,R}}, \end{eqnarray} \tag{ 34a }$$

  $$\begin{eqnarray} (c^*_{L,R})^2 &=& \Gamma _{1,L,R} p^*_{L,R} / \rho ^*_{L,R}, \end{eqnarray} \tag{ 34b }$$

  $$\begin{eqnarray} (\rho e)^*_{L,R} &=& (\rho e)_{L,R} + (p^* - p_{L,R})\frac{(e+p/ \rho)_{L,R}}{c^2_{L,R}},\quad \ \ \ \end{eqnarray} \tag{ 34c }$$

  $$\begin{eqnarray} v^*_{L,R} &=& v_{L,R}, \end{eqnarray} \tag{ 34d }$$

  where v generically represents advected quantities (which includes transverse velocity components). Here, the notation *_L,R refers to values on the left and right side of the contact discontinuity.
• 2.  
  Step 2.4b. The next step in the approximate Riemann solver is to interpret this phase space solution. If u*>0 then the contact discontinuity is moving to the right and numerical flux depends on the speed and structure of the acoustic wave connecting Q_L and Q*_L associated with the λ = u − c eigenvalue. Similarly, if u* < 0 then the contact is moving to the left and the numerical flux depends on the speed and structure of the acoustic wave connecting Q_R and Q*_R associated with the λ = u + c eigenvalue. Here, we discuss in detail the case in which u*>0; the other case is treated analogously.For u*>0, we define λ_L = u_L − c_L and λ*_L = u*_L − c*_L. If p*_L>p_L then the wave is a shock wave and we define a shock speed σ = 1/2(λ_L + λ*_L). For that case if σ>0 then the shock is moving to the right and we define the Godunov state Q_G = Q_L; otherwise Q_G = Q*_L. The rarefaction case is somewhat more complex. If both λ_L and λ*_L are negative, then the rarefaction fan is moving to the left and Q_G = Q*_L. Similarly, if both λ_L and λ*_L are positive, then the rarefaction fan is moving to the right and Q_G = Q_L. However, in the case in which λ_L < 0 < λ*_L, the rarefaction spans the x/t = 0 ray and we need to interpolate the solution. For this case, we define

  $$\begin{equation} {\bf Q}_G = \alpha {\bf Q}^*_L + (1-\alpha) {\bf Q}_L, \end{equation} \tag{ 35 }$$

  where α = λ_L/(λ_L − λ*_L). This choice of α corresponds to linearly interpolating Q through the rarefaction to approximate the state that propagates with zero speed.As noted above the case in which u* < 0 is treated analogously. When u* = 0 we compute Q_G by averaging Q*_L and Q*_R. For the Riemann problem approximation, we allow for user-specified floors for ρ, p, and c to prevent the creation of non-physical values.

The fluxes can then be evaluated from the final Q_G. A small quadratic artificial viscosity that is proportional to the divergence of the velocity field is added to the flux in order to add additional dissipation at strong compressions. We also scale all the species fluxes so that they sum to the density flux, as in the sCMA algorithm described by Plewa & Müller (1999).

6.1.1. Overview

The grid hierarchy is composed of different levels of refinement ranging from coarsest (ℓ = 0) to finest (ℓ = ℓ_finest). The maximum number of levels of refinement allowed, ℓ_max, is specified at the start of a calculation. At any given time in the calculation, there may not be that many levels in the hierarchy, i.e., ℓ_finest can change dynamically as the calculation proceeds as long as ℓ_finest ⩽ ℓ_max. Each level is represented by the union of non-overlapping rectangular grids of a given resolution. Each grid is composed of an even number of cells in each coordinate direction; cells are the same size in each coordinate direction but grids may have different numbers of cells in each direction. Figure 1 shows a cartoon of AMR grids in two dimensions with two levels of refinement.

Zoom In Zoom Out Reset image size

In this implementation, the refinement ratio between levels ℓ and ℓ + 1, which we call r_ℓ, is always two or four, with the same factor of refinement in each coordinate direction. The grids are properly nested, in the sense that the union of grids at level ℓ + 1 is contained in the union of grids at level ℓ. Furthermore, the containment is strict in the sense that, except at physical boundaries, the level ℓ grids are large enough to guarantee that there is a border at least n_proper level ℓ cells wide surrounding each level ℓ + 1 grid (grids at all levels are allowed to extend to the physical boundaries so the proper nesting is not strict there). The parameter n_proper is two for factor 2 refinement, and one for factor 4 refinement, since four ghost cells are needed for the PPM algorithm.

6.1.2. Error Estimation and Regridding

6.1.3. Enlarging the Domain

6.2.1. Overview

6.2.2. Monopole Gravity

6.2.3. Full Gravity Solve

Overview. Solving the Poisson equation for self-gravity on a multilevel grid hierarchy introduces additional complications. We start by defining some necessary notation. We define L^ℓ as an approximation to ∇² at level ℓ, with the assumption that Dirichlet boundary conditions are supplied on the boundary of the union of level ℓ grids (we allow more general boundary conditions at physical boundaries), and define a level solve as the process of solving

$$\begin{eqnarray*} &&L^{\ell } \phi ^{\ell } = 4\pi G\rho ^{\ell }\end{eqnarray*}$$

at level ℓ.

We define L^comp_ℓ,m as the composite grid approximation to ∇² on levels ℓ through m, and define a composite solve as the process of solving

$$\begin{eqnarray*} &&L_{\ell,m}^{\rm comp} \phi ^{\rm comp} = 4\pi G\rho ^{\rm comp} \end{eqnarray*}$$

on levels ℓ through m. The solution to the composite solve satisfies

$$\begin{eqnarray*} &&L^{m} \phi ^{\rm comp} = 4\pi G\rho ^{m}\end{eqnarray*}$$

at level m, but satisfies

$$\begin{eqnarray*} &&L^{{\ell ^\prime }} \phi ^{{\ell ^\prime }} = 4\pi G\rho ^{{\ell ^\prime }}\end{eqnarray*}$$

for ℓ ⩽ ℓ' < m only on the regions of each level not covered by finer grids or adjacent to the boundary of the finer grid region. In regions of a level ℓ' grid covered by level ℓ' + 1 grids, the solution is defined as the volume average of the solution at ℓ' + 1; in level ℓ' cells immediately adjacent to the boundary of the union of level ℓ' + 1 grids, a modified interface operator is used that reflects the geometry of the interface (see, e.g., Almgren et al. 1998 for details of the multilevel cell-centered interface stencil).

In an algorithm without subcycling one can perform a composite solve at every time step, as described in Ricker (2008), to solve for ϕ on all levels. Because the CASTRO algorithm uses subcycling in time, however, we must use level solves at times when the solution is not defined at all levels, and then synchronize the solutions at different levels as appropriate. Even without changes in ρ due to volume averaging and refluxing, replacing a composite solve by separate level solves generates a mismatch in the normal gradient of ϕ at the boundary between each level. We correct these mismatches with a multilevel correction solve, which is a two-level composite solve for a correction to ϕ. In addition to correcting the solutions once the mismatch is detected, we add a correction term to later level solve solutions in order to minimize the magnitude of the correction that will be needed.

Multilevel algorithm. At the start of a calculation, we perform a composite solve from level 0 through ℓ_finest to compute ϕ at all levels. In addition, after every regridding step that creates new grids at level ℓ + 1 and higher, a composite solve from level ℓ through ℓ_finest is used to compute ϕ at those levels.

Following an approach similar to that described in Miniati & Colella (2007), at the start and end of each level ℓ time step we perform a level solve to compute ϕ^ℓ. The difference between ϕ^comp_ℓ and ϕ^ℓ at the start of the time step is stored in ϕ^ℓ,corr. This difference is added to ϕ^ℓ at the beginning and end of this level ℓ time step. Thus, ϕ^ℓ + ϕ^ℓ,corr is identical to ϕ^comp_ℓ at the start of the time step; at the end of the time step it is an approximation to what the solution to the composite solve would be. In the event that the density does not change over the course of the time step, the effect of this lagged correction is to make ϕ^ℓ + ϕ^ℓ,corr at the end of the time step identical to ϕ^comp_ℓ at that time, thus there is no mismatch between levels to correct. In general, when the density is not constant, the effect of the lagged correction is to make the correction solve that follows the end of the time step much quicker. We now describe the two-level correction step. In the discussion below, we will refer to the two levels involved in a correction solve as the "coarse" and "fine" levels.

At the end of r_ℓ level ℓ + 1 time steps, when the level ℓ + 1 solution has reached the same point in time as the level ℓ solution, and after the volume averaging and refluxing steps above have been performed, we define two quantities on the coarse grid. The first is the cell-centered quantity, (δρ)^c, which carries the change in density at the coarse level due only to refluxing. The second is the face-centered flux register,

$$\begin{equation} \delta F_\phi ^\ell = - A^c \frac{\partial \phi ^c}{\partial n} + \sum A^f \frac{\partial \phi ^f}{\partial n}, \end{equation} \tag{ 36 }$$

which accounts for the mismatch in the normal gradient of ϕ at coarse–fine interfaces. Here, A^c and A^f represent area weighting factors on the coarse and fine levels, respectively. We define the composite residual, R^comp, to be zero in all fine cells and in all coarse cells away from the union of fine grids, and

$$\begin{equation} R^{\rm comp} = 4 \pi G (\delta \rho)^c - (\nabla \cdot \delta F_\phi) |^c, \end{equation} \tag{ 37 }$$

on all cells adjacent to the union of fine grids, where (∇·)|^c refers to the discrete divergence at the coarse level, where the only non-zero contribution comes from δF_ϕ on the coarse–fine interface. We then solve

$$\begin{equation} L_{\ell,\ell +1}^{\rm comp} \; \delta \phi = R^{\rm comp} \end{equation} \tag{ 38 }$$

and define the update to gravity at both levels,

$$\begin{equation} \delta {\bf g}= -\nabla (\delta \phi). \end{equation} \tag{ 39 }$$

This update is used to correct the gravitational source terms. We define the new-time state after volume averaging but before refluxing as $(\overline{\rho },\overline{\bf u}, \overline{\rho E},\ldots),$ and the contributions to the solution on the coarse grid from refluxing as ((δρ)^c, δ(ρu)^c, δ(ρE)^c, ...). Then we can define the sync sources for momentum on the coarse and fine levels, S^sync,c_ρu, and S^sync,f_ρu, respectively as follows:

$$\begin{eqnarray*} S^{{\rm sync},c}_{\rho {\bf u}} &=& (\overline{\rho }^c + (\delta \rho)^c) ({\bf g}^{c,n+1} + \delta {\bf g}^{c}) - \overline{\rho }^c \; {\bf g}^{c,n+1} \\ &=& [(\delta \rho)^c {\bf g}^{c,n+1} + (\overline{\rho }^{c} + (\delta \rho)^c) \; \delta {\bf g}^c) ], \\ S^{{\rm sync},f}_{\rho {\bf u}} &=& \overline{\rho }^{f} \; \delta {\bf g}^f. \end{eqnarray*}$$

These momentum sources lead to the following energy sources:

$$\begin{eqnarray*} S^{{\rm sync},c}_{\rho E} &=& S^{{\rm sync},c}_{\rho {\bf u}} \cdot \left(\overline{{\bf u}}^{c} + {1/2}\; \Delta t_c \; S^{{\rm sync},c}_{\rho {\bf u}} / \overline{\rho }^{c} \right), \\ S^{{\rm sync},f}_{\rho E} &=& S^{{\rm sync},f}_{\rho {\bf u}} \cdot \left(\overline{{\bf u}}^{f} + {1/2}\; \Delta t_f \; S^{{\rm sync},f}_{\rho {\bf u}} / \overline{\rho }^{f} \right). \end{eqnarray*}$$

The state at the coarse and fine levels is then updated using

$$\begin{eqnarray*} (\rho {\bf u})^{c,n+1} &=& (\rho {\bf u})^{c} + \delta (\rho {\bf u})^c + {1/2}\Delta t_c S^{{\rm sync},c}_{\rho {\bf u}},\nonumber\\ (\rho {\bf u})^{f,n+1} &=& (\overline{\rho {\bf u}})^{f} + {1/2}\Delta t_f S^{{\rm sync},f}_{\rho {\bf u}}, \\ (\rho E)^{c,n+1} &=& (\rho E)^{c} + \delta (\rho E)^c + {1/2}\Delta t_c S^{{\rm sync},c}_{\rho E},\nonumber\\ (\rho E)^{f,n+1} &=& (\overline{\rho E})^{f} + {1/2}\Delta t_f S^{{\rm sync},f}_{\rho E}. \end{eqnarray*}$$

(The factor of 1/2 follows from the time-centering of the sources.)

To complete the correction step,

• 1.  
  We add δϕ directly to ϕ^ℓ and ϕ^ℓ+1 and interpolate δϕ to any finer levels and add it to the current ϕ at those levels. We note that at this point ϕ at levels ℓ and ℓ + 1 is identical to the solution that would have been computed using a two-level composite solve with the current values of density. Thus the new, corrected, ϕ at each level plays the role of ϕ^comp in the next time step.
• 2.  
  If level ℓ>0, we transmit the effect of this change in ϕ to the coarser levels by updating the flux register between level ℓ and level ℓ − 1. In particular, we set

  $$\begin{equation} \delta {F_\phi }^{\ell -1} = \delta {F_\phi }^{\ell -1} + \sum A^c \frac{\partial (\delta \phi)^{c-f}}{\partial n}. \end{equation} \tag{ 40 }$$

Performance issues. The multilevel algorithm is not as computationally expensive as it might appear. Because multigrid is an iterative solver, the cost of each solve is proportional to the number of V-cycles, which is a function of the desired reduction in residual. We can reduce the number of V-cycles needed in two ways. First, we can supply a good initial guess for the solution; second, we can lower the desired reduction in residual.

In the case of level solves, we always use ϕ from a previous level solve, when available, as a guess in the current level solve. Thus, even in a single-level calculation, ϕ from the beginning of the time step is used as a guess in the level solve at the end of the time step. If no regridding occurs, then ϕ at the end of one time step can be used as a guess for the level solve at the start of the next time step. The extent to which ρ changes in a time step dictates the extent to which a new computation of gravity is needed, but this also dictates the cost of the update.

Similarly, there is no point in solving for δϕ to greater accuracy than we solve for ϕ. When we do the correction solve for δϕ, we require only that the residual be reduced to the magnitude of the final residual from the level solve, not that we reduce the correction residual by the same factor. Thus, if the right-hand side for the correction solve is already small, the cost of the correction solve will be significantly less than that of the initial level solve.

CASTRO is implemented within the BoxLib framework, a hybrid C++/Fortran90 software system that provides support for the development of parallel structured-grid AMR applications. The basic parallelization strategy uses a hierarchical programming approach for multicore architectures based on both MPI and OpenMP. In the pure-MPI instantiation, at least one grid at each level is distributed to each core, and each core communicates with every other core using only MPI. In the hybrid approach, where on each socket there are n cores which all access the same memory, we can instead have one larger grid per socket, with the work associated with that grid distributed among the n cores using OpenMP.

In BoxLib, memory management, flow control, parallel communications, and I/O are expressed in the C++ portions of the program. The numerically intensive portions of the computation, including the multigrid solvers, are handled in Fortran90. The fundamental parallel abstraction in both the C++ and the Fortran90 is the MultiFab, which holds the data on the union of grids at a level. A MultiFab is composed of FABs; each FAB is an array of data on a single grid. During each MultiFab operation the FABs composing that MultiFab are distributed among the cores. MultiFabs at each level of refinement are distributed independently. The software supports two data distribution schemes, as well as a dynamic switching scheme that decides which approach to use based on the number of grids at a level and the number of processors. The first scheme is based on a heuristic knapsack algorithm as described in Crutchfield (1991) and in Rendleman et al. (2000). The second is based on the use of a Morton-ordering space-filling curve.

Each processor contains metadata that is needed to fully specify the geometry and processor assignments of the MultiFabs. At a minimum, this requires the storage of an array of boxes specifying the index space region for each AMR level of refinement. One of the advantages of computing with fewer, larger grids in the hybrid OpenMP–MPI approach is that the size of the metadata is substantially reduced.

Data for checkpoints and analysis are written in a self-describing format that consists of a directory for each time step written. Checkpoint directories contain all necessary data to restart the calculation from that time step. Plotfile directories contain data for postprocessing, visualization, and analytics, which can be read using amrvis, a customized visualization package developed at LBNL for visualizing data on AMR grids, or VisIt (VisIt User's Manual 2005). Within each checkpoint or plotfile directory is an ASCII header file and subdirectories for each AMR level. The header describes the AMR hierarchy, including number of levels, the grid boxes at each level, the problem size, refinement ratio between levels, step time, etc. Within each level directory are the MultiFab files for each AMR level. Checkpoint and plotfile directories are written at user-specified intervals.

For output, each processor writes its own data to the appropriate MultiFab files. The output streams are coordinated to allow only one processor to write to a file at one time and to try to maintain maximum performance by keeping the number of open data streams, which is set at run time, equal to the number of files being written. Data files typically contain data from multiple processors, so each processor writes data from its associated grid(s) to one file, then another processor can write data from its associated grid(s) to that file. A designated I/O Processor writes the header files and coordinates which processors are allowed to write to which files and when. The only communication between processors is for signaling when processors can start writing and for the exchange of header information. We also use the C++ setbuf function for good single file performance. While I/O performance even during a single run can be erratic, recent timings on the Franklin machine (XT4) at NERSC indicate that CASTRO's I/O performance, when run with a single level composed of multiple uniformly sized grids, matches some of the top results for the N5 IOR benchmark (roughly 13 GB s⁻¹; Franklin Performance Monitoring 2010). For more realistic simulations with multiple grids at multiple levels, CASTRO is able to write data at approximately 5 GB s⁻¹ sustained, over half of the average I/O benchmark reported speed.

Restarting a calculation can present some difficult issues for reading data efficiently. In the worst case, all processors would need data from all files. If multiple processors try to read from the same file at the same time, performance problems can result, with extreme cases causing file system thrashing. Since the number of files is generally not equal to the number of processors and each processor may need data from multiple files, input during restart is coordinated to efficiently read the data. Each data file is only opened by one processor at a time. The IOProcessor creates a database for mapping files to processors, coordinates the read queues, and interleaves reading its own data. Each processor reads all data it needs from the file it currently has open. The code tries to maintain the number of input streams to be equal to the number of files at all times.

Checkpoint and plotfiles are portable to machines with a different byte ordering and precision from the machine that wrote the files. Byte order and precision translations are done automatically, if required, when the data are read.

In Figure 2, we show the scaling behavior of the CASTRO code, using only MPI-based parallelism, on the jaguarpf machine at the Oak Ridge Leadership Computing Facility (OLCF). A weak scaling study was performed, so that for each run there was exactly one 64³ grid per processor. We ran the code with gravity turned off, with the monopole approximation to gravity, and with the Poisson solve for gravity. The monopole approximation to gravity adds very little to the run time of the code; with and without the monopole approximation the code scales excellently from 8 to 64,000 processors. For the 64,000 processor case without gravity, the time for a single core to advance one cell for one time step is 24.8 μs.

Zoom In Zoom Out Reset image size

Good scaling of linear solves is known to be much more difficult to achieve; we report relatively good scaling up to only 13,824 processors in the pure-MPI approach. An early strong scaling study contrasting the pure-MPI and the hybrid-MPI–OpenMP approaches for a 768³ domain shows that one can achieve at least a factor of 3 improvement in linear solver time by using the hybrid approach at large numbers of processors. Improving the performance of the linear solves on the new multicore architectures is an area of active research; more extensive development and testing is underway.

We also ran a scaling study with a single level of local refinement using the monopole gravity approximation. In this MPI-only study, there is one 64³ grid at each level for each processor. Because of subcycling in time, a coarse time step consists of a single step on the coarse grid and two steps on the fine grid. Thus, we would expect that the time to advance the multilevel solution by one coarse time step would be a factor of 3 greater than the time to advance the single-level coarse solution by one coarse time step, plus any additional overhead associated with AMR. From the data in the figure, we conclude that AMR introduces a modest overhead, ranging from approximately 5% for the 8 processor case to 19% for the 64,000 processor case. By contrast, advancing a single-level calculation at the finer resolution by the same total time, i.e., two fine time steps, would require a factor of 16 more resources than advancing the coarse single-level solution.

To test the behavior of the hydrodynamics solver, we run several different one-dimensional shock-tube problems. The setup for these problems consists of a left and right state, with the interface in the center of the domain. All calculations use a gamma-law EOS with γ = 1.4. We show results from each problem run using one-dimensional Cartesian coordinates, but we have verified that the results are identical when each problem is run in two-dimensional or three-dimensional Cartesian coordinates and the interface is normal to a coordinate axis. The length of the domain is always taken as 1.0, with the interface in the center. We use a base grid of 32 cells, with two additional levels of factor 2 refinement, for an effective resolution of 128 cells. The refinement criteria are based on gradients of density and velocity. In the case of the double rarefaction, we also present results from runs with two levels of factor 4 refinement (effective resolution of 512 cells) and three levels of factor 4 refinement (effective resolution of 2048 cells). In each case, analytic solutions are found using the exact Riemann solver from Toro (1997). All calculations are run with the new PPM limiters and a CFL number of 0.9. For each problem, we show density, pressure, velocity, and internal energy.

8.1.1. Sod's Problem

Sod's problem (Sod 1978) is a simple shock-tube problem that exhibits a shock, contact discontinuity, and a rarefaction wave. The non-dimensionalized initial conditions are

$$\begin{equation} \begin{array}{l} \rho _L = 1 \\ u_L = 0 \\ p_L = 1 \end{array} \qquad \begin{array}{l} \rho _R = 0.125, \\ u_R = 0, \\ p_R = 0.1. \end{array} \end{equation} \tag{ 41 }$$

This results in a rightward moving shock and contact discontinuity, and a leftward moving rarefaction wave. Figure 3 shows the resulting pressure, density, velocity, and internal energy at t = 0.2 s. We see excellent agreement with the exact solution.

Zoom In Zoom Out Reset image size

8.1.2. Double Rarefaction

The double rarefaction problem tests the behavior of the hydrodynamics algorithm in regions where a vacuum is created. We run the problem as described in Toro (1997). The non-dimensionalized initial conditions are

$$\begin{equation} \begin{array}{l} \rho _L = 1 \\ u_L = -2 \\ p_L = 0.4 \end{array} \qquad \begin{array}{l} \rho _R = 1, \\ u_R = 2, \\ p_R = 0.4. \end{array} \end{equation} \tag{ 42 }$$

This results in two rarefaction waves propagating in opposite directions away from the center. As a result, matter is evacuated from the center, leaving behind a vacuum. Figure 4 shows the CASTRO solutions at t = 0.15 s. The agreement with the exact solution is excellent at the 128-cell resolution for density, pressure, and velocity; the internal energy is more sensitive, but clearly converges to the analytic solution except at the center line. This is a very common pathology for this problem, since the internal energy, e, is derived from Equation (7) using values of p and ρ which are both approaching zero in the center of the domain (Toro 1997).

Zoom In Zoom Out Reset image size

8.1.3. Strong Shock

The final shock-tube problem we try is a strong shock. We initialize the problem as described in Toro (1997). The initial conditions are

$$\begin{equation} \begin{array}{l} \rho _L = 1 \\ u_L = 0 \\ p_L = 1000 \end{array} \qquad \begin{array}{l} \rho _R = 1, \\ u_R = 0, \\ p_R = 0.01. \end{array} \end{equation} \tag{ 43 }$$

The initial pressure jump of 6 orders of magnitude results in a strong rightward moving shock. This large dynamic range can cause trouble for some hydrodynamics solvers. The shock is followed very closely by a contact discontinuity. A leftward moving rarefaction is also present. Figure 5 shows the CASTRO results at t = 0.012 s. We see good agreement between the CASTRO results and the exact solution.

Zoom In Zoom Out Reset image size

Another standard hydrodynamics test is the Sedov–Taylor blast wave. The problem setup is very simple: a large amount of energy is deposited into the center of a uniform domain. This drives a blast wave (spherical or cylindrical, depending on the domain geometry). An analytic solution is provided by Sedov (1959). We use a publicly available code described by Kamm & Timmes (2007) to generate the exact solutions.

The Sedov explosion can test the geometrical factors in the hydrodynamics scheme. A cylindrical blast wave (e.g., a point explosion in a two-dimensional plane) can be modeled in two-dimensional Cartesian coordinates. A spherical blast wave can be modeled in one-dimensional spherical, two-dimensional axisymmetric (cylindrical r–z), or three-dimensional Cartesian coordinates.

In the Sedov problem, the explosion energy, $\mathcal {E}_{\rm exp}$ (in units of energy, not energy/mass or energy/volume), is deposited into a single point, in a medium of uniform ambient density, ρ_ambient, and pressure, p_ambient. Initializing the problem can be difficult because the small volume is typically only one cell in extent, which can lead to grid imprinting in the solution. A standard approach (see, for example, Fryxell et al. 2000, Omang et al. 2006, and the references therein) is to convert the explosion energy into a pressure contained within a certain volume, V_init, of radius r_init as

$$\begin{equation} p_\mathrm{init} = \frac{(\gamma - 1) \mathcal {E}_{\rm exp}}{V_\mathrm{init}}. \end{equation} \tag{ 44 }$$

This pressure is then initialized to p_init in all of the cells where r < r_init. We use the gamma-law EOS with γ = 1.4.

To further minimize any grid effects, we do subsampling in each cell: each cell is divided it into N_sub subcells in each coordinate direction, each subcell is initialized independently, and then the subcells are averaged together (using volume weighting for spherical or cylindrical coordinates) to determine the initial state of the full cell.

For these runs, we use ρ_ambient = 1 g cm⁻³, p_ambient = 10⁻⁵ dyn cm⁻², $\mathcal {E}_{\rm exp}= 1\,\mathrm{erg}$, r_init = 0.01 cm, and N_sub = 10. A base grid with Δx = 0.03125 cm is used with three levels of factor 2 refinement. For most geometries, we model the explosion in a domain ranging from 0 to 1 cm in each coordinate direction. In this case, the base grid would have 32 cells in each coordinate direction and the finest mesh would correspond to 256 cells in each coordinate direction. For the two-dimensional axisymmetric case, we model only one quadrant, and the domain ranges from 0 to 0.5 cm. All calculations were run with a CFL number of 0.5, and the initial time step was shrunk by a factor of 100 to allow the point explosion to develop. We refine on regions where ρ>3 g cm⁻³, ∇ρ>0.01 g cm⁻³ cm⁻¹, p>3 dyn cm⁻², or ∇p>0.01 dyn cm⁻² cm⁻¹.

Figure 6 shows the CASTRO solution to a spherical Sedov explosion at time t = 0.01 s, run in one-dimensional spherical, two-dimensional cylindrical, and three-dimensional Cartesian coordinates. For the two-dimensional and three-dimensional solutions, we compute the radial profile by mapping each cell into its corresponding radial bin and averaging. The radial bin width was picked to match the width of a cell at the finest level of refinement in the CASTRO solution. The density, velocity, and pressure plots match the exact solution well. As with the double rarefaction problem, the internal energy is again the most difficult quantity to match due to the vacuum region created at the origin. Figure 7 shows the same set of calculations run with four levels of factor 2 refinement. Here, the agreement is even better. Figure 8 shows the CASTRO solution at time t = 0.1 s to a cylindrical Sedov explosion, run in two-dimensional Cartesian coordinates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The R–T instability results when a dense fluid is placed over a less-dense fluid in a gravitational field (Taylor 1950; Layzer 1955; Sharp 1984). The interface is unstable and a small perturbation will result in the growth a buoyant uprising bubbles and dense, falling spikes of fluid. This instability provides a mechanism for mixing in many astrophysical systems. Despite its seemingly simplistic nature, only the linear growth regime is understood analytically (see, for example, Chandrasekhar 1961). In the nonlinear regime, R–T instability calculations are often used as a means of code validation (Dimonte et al. 2004).

For our purposes, the R–T instability provides a good basis to compare different choices of the advection algorithm. We model a single-mode R–T instability—a perturbation consisting of a single wavelength that disturbs the initial interface. Short-wavelength perturbations have a faster growth rate than long-wavelength perturbations, so grid effects can easily drive the instability on smaller scales than our initial perturbation. No viscous terms are explicitly modeled.

We choose the density of the dense fluid to be ρ₂ = 2 g cm⁻³ and the light fluid is ρ₁ = 1 g cm⁻³. The gravitational acceleration is taken to be g = −1 cm s⁻² in the vertical direction. The gamma-law EOS is used with γ = 1.4. Our domain has a width of L_x = 0.5 cm and a height of L_y = 1 cm. The initial interface separating the high and low-density fluid is centered vertically at L_y/2, with the density in the top half taken to be ρ₂ and the density in the lower half ρ₁. Since g and ρ₁, ρ₂ are constant, we can analytically integrate the equation of hydrostatic equilibrium to get the pressure in both the high- and low-density regions of the domain:

$$\begin{equation} p(y) = \left\lbrace \begin{array}{@{}ll@{}} p_\mathrm{base} + \rho _1 g y & y < L_y / 2 \\ p_\mathrm{base} + \rho _1 g L_y/2 + \rho _2 g (y - L_y/2) \quad & y > L_y/2, \end{array} \right. \end{equation} \tag{ 45 }$$

where y is the vertical coordinate and p_base is the pressure at the base of the domain. We take p_base = 5 dyn cm⁻².

To initiate the instability, the interface is perturbed by slightly shifting the density, keeping the interface centered vertically in the domain. We define the perturbed interface height, ψ, to be a function of position in the x-direction as

$$\begin{equation} \psi (x) = \frac{A}{2} \left[ \cos \left(\frac{2\pi x}{L_x} \right) + \cos \left(\frac{2\pi (L_x - x)}{L_x} \right) \right] + \frac{L_y}{2}, \end{equation} \tag{ 46 }$$

with the amplitude A = 0.01 cm. We note that the cosine part of the perturbation is done symmetrically, to prevent roundoff error from introducing an asymmetry in the flow. The density is then perturbed as

$$\begin{equation} \rho (x,y) = \rho _1 + \frac{\rho _2 - \rho _1}{2} \left[ 1 + \tanh \left(\frac{y - \psi (x)}{h} \right) \right]. \end{equation} \tag{ 47 }$$

The tanh  profile provides a slight smearing of the initial interface, over a smoothing length h. We take h = 0.005 cm.

In Figure 9, we present simulation results for the R–T problem at t = 2.5 s for several different variants of the hydrodynamics. All calculations were run with 256 × 512 grid cells. In the bottom right image, we show the results obtained using the unsplit PPM with the new limiter used in CASTRO. The left and middle images on the bottom row are results using the unsplit piecewise-linear method and unsplit PPM with limiters as in Miller & Colella (2002), respectively. The results with all three methods are reasonably good; however, the piecewise linear and original PPM limiter both exhibit mild anomalies at the tip of both the bubble and the spike.

Zoom In Zoom Out Reset image size

In the upper row, we present results for the R–T problem using operator-split analogs of the unsplit methods. The details of the algorithms such as limiters, Riemann solver, etc., are the same as in the unsplit methods; the only difference is the use of operator splitting. We note that all three of the operator-split methods produce spurious secondary instabilities. This behavior is a direct result of the operator-split approach. Physically, for these low Mach number flows, the density field is advected by a nearly incompressible flow field, and remains essentially unchanged along Lagrangian trajectories. However, in regions where there is significant variation in the local strain rate, an operator-split integration approach alternately compresses and expands the fluid between subsequent sweeps. This alternating compression/expansion provides the seed for the anomalies observed with operator-split methods.

We note that both the CPU time and the memory usage are roughly a factor of 2 larger for the unsplit algorithm than for the split algorithm in this two-dimensional implementation. For a pure-hydrodynamics problem with gamma-law EOS this factor is nontrivial; for a simulation that uses the full self-gravity solver, a realistic reaction network, a costly EOS, or significant additional physics, the additional cost of the hydrodynamic solver may be negligible.

In three dimensions, one might expect the ratio of CPU time for the unsplit algorithm relative to the split algorithm to be even larger than in two dimensions because of the additional Riemann solves required to construct the transverse terms. However, this effect is counterbalanced by the need to advance ghost cells in the split algorithm to provide boundary conditions for subsequent sweeps. Consequently, we observe an increase in CPU time that is slightly less than the factor of 2 observed in two dimensions. The three-dimensional implementation of the unsplit algorithm in CASTRO uses a strip-mining approach that only stores extra data on a few planes at a time, so we see an increase of less than 10% in the memory required for the unsplit integrator compared to the split integrator in three dimensions.

A challenging problem for a hydrodynamics code is to keep a star in hydrostatic equilibrium. Because of the different treatment of the pressure, density, and gravitational acceleration by the hydrodynamics algorithm, small motions can be driven by the inexact cancellation of ∇p and ρg. This is further exaggerated by modeling a spherical star on a three-dimensional Cartesian grid. Here, we test the ability of CASTRO to maintain hydrostatic equilibrium for a spherical, self-gravitating star.

Our initial model is a nearly Chandrasekhar mass, carbon–oxygen white dwarf, which is generated by specifying a core density (2.6 × 10⁹ g cm⁻³), 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{eqnarray} p_{0,j+1} &=& p_{0,j} + \frac{1}{2} \Delta r(\rho _{0,j} + \rho _{0,j+1}) g_{j+{1/2}}, \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} s_{0,j+1} &=& s_{0,j}, \end{eqnarray} \tag{ 49 }$$

with Δr = 1.653125 × 10⁵ cm. We begin with a guess of ρ_0,j+1 and T_0,j+1 and use the EOS 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 until the density falls to 10⁻⁴ g cm⁻³, after which we hold the density constant. This uniquely determines the initial model. We note that this is the same procedure we follow to initialize a convecting white dwarf for the multilevel low Mach number code, MAESTRO, described in Nonaka et al. (2010).

We map the model onto a (5 × 10⁸ cm)³ domain with 192³, 384³, and 768³ grid cells, and center the star in the domain. We let the simulation run to 1 s, and compare the maximum magnitude of velocity versus time and the magnitude of velocity versus radius at t = 1 s, a time greater than two sound-crossing times. We only consider regions of the star at r < 1.8 × 10⁸ cm, which corresponds to a density of ρ ≈ 5.4 × 10⁵ g cm⁻³. Note that the density reaches the floor of 10⁻⁴ g cm⁻³ at r = 1.9 × 10⁸ cm. We turn on the sponge at the radius where ρ = 100 g cm⁻³ and the sponge reaches its full strength at the radius where ρ = 10⁻⁴ g cm⁻³ with a sponge strength of κ = 1000 s⁻¹. We use a CFL of 0.9 and no refinement. We use the Helmholtz EOS (Timmes & Swesty 2000; Fryxell et al. 2000) and no reactions are modeled.

Figure 10 shows a plot of the maximum magnitude of velocity versus time. At each of the three resolutions, we show the results using a monopole gravity approximation and Poisson solve for gravity. We note that in each simulation, the maximum velocity is not strictly increasing, leading us to believe that over longer periods of time the velocities will remain small. We note that sound speed at the center of the star is approximately 9.4 × 10⁸ cm s⁻¹, so at the highest resolution, the peak velocity is less than 1% of the sound speed. The monopole and Poisson cases match up very well, except for the finest resolution. The reason why we see larger peak velocities in the finest resolution Poisson solver simulation is due to the large velocities at the edge of the star.

Zoom In Zoom Out Reset image size

Figure 11 shows a plot of the magnitude of velocity versus radius at t = 1 s. Again, at each of the three resolutions, we show the results using a monopole gravity approximation and Poisson solve for gravity. Here, we see clear second-order convergence in the max norm, and the monopole and Poisson simulations agree best at the highest resolution. We also see how in the finest resolution runs, the velocities at the edge of the star can become large, but this is likely outside the region of interest for a typical simulation.

Zoom In Zoom Out Reset image size

As a second test of the gravity solver in CASTRO, we implement the homologous dust collapse test problem, a "pressure-less" configuration that collapses under its own self-gravity. An analytic solution that describes the radius of the sphere as a function of time is found in Colgate & White (1966). Our implementation of this problem follows that described in FLASH 3.2 User's Guide (2009) and Monchmeyer & Muller (1989). The problem is initialized with a sphere with a large, uniform density, ρ₀, of radius r₀. The pressure everywhere should be negligible, i.e., the sound-crossing time should be much longer than the free-fall collapse time (see, for example, FLASH 3.2 User's Guide 2009). Colgate & White (1966) use p = 0. We choose a value that does not appear to affect the dynamics. As the sphere collapses, the density inside should remain spatially constant, but increase in value with time.

Following FLASH 3.2 User's Guide (2009), we take ρ₀ = 10⁹ g cm⁻³ and r₀ = 6.5 × 10⁸ cm. The pressure is not specified, so we take it to be 10¹⁵ dyn cm⁻². Outside of the sphere, we set the density to ρ_ambient = 10⁻⁵ g cm⁻³. Finally, since the sharp cutoff at the edge of the sphere is unphysical, we smooth the initial profile by setting

$$\begin{equation} \rho = \rho _0 - \frac{\rho _0 - \rho _\mathrm{ambient}}{2} \left[ 1 + \tanh \left(\frac{r - r_0}{h} \right) \right] \end{equation} \tag{ 50 }$$

with the smoothing length, h = 4 × 10⁶ ≪ r₀. We use the gamma-law EOS with γ = 1.66.

Figure 12 shows the radius versus time for the one-dimensional, two-dimensional, and three-dimensional simulations as compared to the exact solution. In all three cases, we see excellent agreement with the exact solution.

Zoom In Zoom Out Reset image size

A final test is a code comparison of the evolution of three reacting bubbles in a plane-parallel stellar atmosphere. This problem is almost identical to the setup described in Section 4.2 of Almgren et al. (2008) with two minor differences. First, we eliminate the stably stratified layer at the base of the atmosphere by setting the lower y extrema of the domain to 5.00625 × 10⁷ cm—this way, the bottom-most row of cells in the domain is initialized with the specified base density (2.6 × 10⁹ g cm⁻³) and temperature. Second, we set the base temperature of the atmosphere to 6 × 10⁸ K (instead of 7 × 10⁸ K) to minimize the amount of reactions occurring near the lower domain boundary. Three temperature perturbations are seeded in pressure equilibrium with a range of heights and widths as specified by Equation (87) and Table 1 of Almgren et al. (2008). We use a uniform computation grid of 384 × 576 cells and a domain width of 2.16 × 10⁸ cm.

We compare the evolution to the FLASH code (Fryxell et al. 2000), version 2.5, using the standard dimensionally split PPM hydrodynamics module that comes with FLASH. The lower boundary condition in both cases provides hydrostatic support by integrating the equation of hydrostatic equilibrium together with the equations of state into the ghost cells, assuming a constant temperature, as described in Zingale et al. (2002). The left and right boundary is periodic. We use the same single step (¹²C + ¹²C → ²⁴Mg) reaction module described in Almgren et al. (2008). Both codes use the general stellar EOS described in Fryxell et al. (2000) and Timmes & Swesty (2000) with the Coulomb corrections enabled.

Figures 13 and 14 show contours of the temperature and X(²⁴Mg) after 2.5 s of evolution for both FLASH and CASTRO. We see excellent agreement between the two codes in terms of bubble heights and contour levels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a final example, in Figure 15 we show a two-dimensional snapshot of temperature from a three-dimensional calculation of a Type Ia supernova (H. Ma & A. J. Aspden 2010, private communication; Ma et al. 2010). This simulation uses a realistic stellar EOS and a turbulent flame model, and is typical of more realistic CASTRO applications. The domain is 5.12 ×10⁸ cm on a side, and is covered with 512 64³ grids. There are two levels of factor 2 refinement, with approximately 1.8% of the domain covered by level 2 grids with an effective resolution of 2.5 ×10⁵ cm. Figure 16 is a close-up of the center of the domain so that the level 2 grids are more visible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size