Fornax: A Flexible Code for Multiphysics Astrophysical Simulations

As is standard practice, we denote the contravariant components of a vector with raised indices, covariant components with lowered indices, covariant differentiation with a semicolon, and partial differentiation with a comma, and make use of Einstein notation for summation over repeated indices. Here and throughout this work, we adopt a coordinate basis. In this notation, the equations of Newtonian hydrodynamics can be written

$$\begin{eqnarray}&&{\rho }_{,t}+{(\rho {v}^{i})}_{;i}=0,\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{(\rho {v}_{j})}_{,t}+{(\rho {v}^{i}{v}_{j}+P{{\delta }^{i}}_{j})}_{;i}={S}_{j},\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&{\left(\rho e\right)}_{,t}+{\left[\rho {v}^{i}\left(e+\displaystyle \frac{P}{\rho }\right)\right]}_{;i}={S}_{E},\end{eqnarray} \tag{ 1c }$$

$$\begin{eqnarray}&&{(\rho X)}_{,t}+{(\rho {{Xv}}^{i})}_{;i}={S}_{X},\end{eqnarray} \tag{ 1d }$$

where e is the specific total energy of the gas, X is an arbitrary scalar quantity that may represent, e.g., composition, and S_j, S_E, and S_X are source terms that account for additional physics. The contravariant components of the velocity are ${v}^{i}={{dx}}^{i}/{dt}$, i.e., they are coordinate velocities. In a coordinate basis, the covariant derivatives can be expanded to yield

$$\begin{eqnarray}&&{\rho }_{,t}+\displaystyle \frac{1}{\sqrt{g}}{(\sqrt{g}\rho {v}^{i})}_{,i}=0,\end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&{(\rho {v}_{j})}_{,t}+\displaystyle \frac{1}{\sqrt{g}}{\left[\sqrt{g}{{T}^{i}}_{j}\right]}_{,i}={{\rm{\Gamma }}}_{{jk}}^{l}{{T}^{k}}_{l}+{S}_{j},\end{eqnarray} \tag{ 2b }$$

$$\begin{eqnarray}&&{\left(\rho e\right)}_{,t}+\displaystyle \frac{1}{\sqrt{g}}{\left[\sqrt{g}\rho {v}^{i}\left(e+\displaystyle \frac{P}{\rho }\right)\right]}_{,i}={S}_{E},\end{eqnarray} \tag{ 2c }$$

$$\begin{eqnarray}&&{(\rho X)}_{,t}+\displaystyle \frac{1}{\sqrt{g}}{\left(\sqrt{g}\rho {{Xv}}^{i}\right)}_{,i}={S}_{X},\end{eqnarray} \tag{ 2d }$$

where g is the determinant of the metric, ${{\rm{\Gamma }}}_{\ {jk}}^{l}$ are the Christoffel symbols, defined in terms of derivatives of the metric, and ${{T}^{k}}_{l}\equiv \rho {v}^{i}{v}_{j}+P{{\delta }^{i}}_{j}$ is the fluid stress tensor.

Note that we have chosen to express the momentum equation as a conservation law for the covariant components of the momentum. There is good reason to do so. Written this way, the geometric source terms, ${{\rm{\Gamma }}}_{{jk}}^{l}{{T}^{k}}_{l}$, vanish identically for components associated with ignorable coordinates in the metric. A good example is in spherical (r, θ, ϕ) coordinates, where, since ϕ does not explicitly enter into the metric, the geometric source terms vanish for the ρv_ϕ equation. Physically, ρv_ϕ is the angular momentum, so we are left with an explicit expression of angular momentum conservation that the numerics will satisfy to machine precision, rather than to the level of truncation error. In general, the covariant expression of the momentum equation respects the geometry of the problem without special consideration or coordinate-specific modifications of the code.

Fornax evolves the zeroth and first moments of the frequency-dependent comoving-frame radiation transport equation. Keeping all terms to ${ \mathcal O }(v/c)$ and dropping terms proportional to the fluid acceleration, the monochromatic radiation moment equations can be written

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\nu ,t}+{({F}_{\nu }^{i}+{v}^{i}{E}_{\nu })}_{;i}\\ +\,{{v}^{i}}_{;j}\left[{P}_{\nu i}^{j}-{\partial }_{\nu }\left(\nu {P}_{\nu i}^{j}\right)\right]={R}_{\nu E},\end{array}\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&\begin{array}{l}{F}_{\nu j,t}+{({c}^{2}{P}_{\nu j}^{i}+{v}^{i}{F}_{\nu j})}_{;i}\\ +\,{{v}^{i}}_{;j}{F}_{\nu i}-{{v}^{i}}_{;k}\,{\partial }_{\nu }\left(\nu {Q}_{\nu {ji}}^{k}\right)={R}_{\nu j}.\end{array}\end{eqnarray} \tag{ 3b }$$

In Equations (3), E_ν and F_νj denote the monochromatic energy density and flux of the radiation field at frequency ν in the comoving frame, ${P}_{\nu i}^{j}$ is the radiation pressure tensor (second moment), ${Q}_{\nu {ji}}^{k}$ is the heat-flux tensor (third moment), and R_νE and R_νj are source terms that account for interactions between the radiation and matter. These interaction terms are written

$$\begin{eqnarray}&&{R}_{\nu E}={j}_{\nu }-c{\kappa }_{\nu }{E}_{\nu },\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{R}_{\nu j}=-c({\kappa }_{\nu }+{\sigma }_{\nu }){F}_{\nu j},\end{eqnarray} \tag{ 4b }$$

where j_ν is the emissivity, ${\kappa }_{\nu }$ is the absorption coefficient, and σ_ν is the scattering coefficient. Correspondingly, there are energy and momentum source terms in the fluid equations:

$$\begin{eqnarray}&&{S}_{j}=-\displaystyle \frac{1}{{c}^{2}}{\int }_{0}^{\infty }{R}_{\nu j}\,d\nu ,\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{S}_{E}=-{\int }_{0}^{\infty }\left({R}_{\nu E}+\displaystyle \frac{{v}^{i}}{{c}^{2}}{R}_{\nu i}\right)d\nu .\end{eqnarray} \tag{ 5b }$$

As in the fluid sector, we rewrite these covariant derivatives as partial derivatives, introducing geometric source terms in the radiation momentum equation.

Fornax can treat either photon or neutrino radiation fields. For neutrinos, Equations (3) are solved separately for each species, and Equations (5) are summed over species. Additionally, the electron fraction is evolved according to Equation 2(d), with X = Y_e and

$$\begin{eqnarray}&&{S}_{X}=\displaystyle \sum _{s}{\int }_{0}^{\infty }{\xi }_{s\nu }({j}_{s\nu }-c{\kappa }_{s\nu }{E}_{s\nu })\,d\nu ,\end{eqnarray} \tag{ 6 }$$

where s refers to the neutrino species and

$$\begin{eqnarray}{\xi }_{s\nu }=\left\{\begin{array}{cc}-{\left({N}_{{\rm{A}}}\nu \right)}^{-1}, & s={\nu }_{e}\\ {\left({N}_{{\rm{A}}}\nu \right)}^{-1}, & s={\bar{\nu }}_{e}\\ 0, & s={\nu }_{x}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

where N_A is Avogadro's number.

The hydrodynamics in Fornax is evolved using a directionally unsplit, high-order Godunov-type scheme. Having already described the underlying spatial discretization of the variables as well as the time-stepping scheme, the essential remaining element is a method for computing fluxes on cell faces. Since Fornax employs Runge–Kutta time stepping, there is no need for the characteristic tracing step or transverse flux gradient corrections required in single-step unsplit schemes. This approach has the great advantage of relative simplicity, especially in the multiphysics context. Fornax follows the standard approach used by similar codes, first reconstructing the state profile at cell faces and then solving the resulting Riemann problem to obtain fluxes.

There are many potential approaches to reconstructing face data based on cell values. One particularly popular approach is to reconstruct the primitive variables (ρ, P, and vⁱ) using the method described by Woodward & Colella (1984), the so-called piecewise-parabolic method (PPM). This method is based on the idea of reconstructing profiles of volume-averaged data, with curvilinear coordinates incorporated by reconstructing the profiles in volume coordinates rather than in the original coordinates themselves. This approach can be problematic in the vicinity of coordinate singularities. Blondin & Lufkin (1993) suggested a generalized approach to PPM in cylindrical and spherical coordinates. Unfortunately, it can be shown that their approach cannot be used generically and, in fact, fails even for standard spherical coordinates. Here, we present an alternative approach that works in arbitrary geometries and coordinates, which contains PPM as a special case.

We begin by writing down an expression for a general nth-order polynomial:

$$\begin{eqnarray}&&p(x)=\displaystyle \sum _{j=0}^{n}{c}_{j}{x}^{j},\end{eqnarray} \tag{ 13 }$$

where c_j is the coefficient of the $j\mathrm{th}$-order monomial. As in PPM, we wish to construct a polynomial p(x) consistent with p_i, the volume-average of the quantity p in cell i. In our formulation, the volume-average of p(x) must then satisfy a set of constraint equations of the form

$$\begin{eqnarray}\begin{array}{rcl}{p}_{i} & = & \langle p{\rangle }_{i}\\ & = & \displaystyle \frac{1}{V}{\displaystyle \int }_{{x}_{i-1/2}}^{{x}_{i+1/2}}p(x)\sqrt{g}\,{dx}\\ & = & \displaystyle \sum _{j=0}^{n}{c}_{j}\langle {x}^{j}{\rangle }_{i},\end{array}\end{eqnarray} \tag{ 14 }$$

where ${x}_{i\pm 1/2}\equiv {x}_{i}\pm {\rm{\Delta }}x/2$ and the quantities $\langle {x}^{j}{\rangle }_{i}$ can be precomputed for each cell either directly or via Romberg integration, then stored in one-dimensional arrays. For the case n = 0, i.e., for piecewise-constant reconstruction, this yields the trivial result that p(x) = p_i. For n > 0, we must further constrain the reconstruction. To proceed, we first distinguish between cases where n is even and n is odd. For odd n, we form a symmetric stencil around each edge containing data in n + 1 cells and require that $\langle p{\rangle }_{i}={p}_{i}$ in each cell i. This necessarily yields continuous left and right states at cell interfaces, but that continuity may be broken by the subsequent application of a slope limiter. For even n, we similarly constrain the polynomial within a symmetric stencil of n + 1 cells, this time centered on a given cell, yielding potentially discontinuous data at cell faces even before the application of a slope limiter.

In Fornax, since we most often employ third-order piecewise-parabolic reconstruction (n = 2), we now describe it in more detail. The reconstruction of p(x) in cell i depends on the data values p_i−1, p_i, and p_i+1. Note that this requires the same number of ghost cells as linear reconstruction, is of similar cost, and being of higher order, produces superior results for many problems, as we show below and in Section 8.

To begin, if p_i is a local extremum of the data, then the reconstruction in cell i is taken to be piecewise-constant. Otherwise, three constraints in the form of Equation (14) are used to form the linear system

$$\begin{eqnarray}\left[\begin{array}{lll}1 & \langle x{\rangle }_{i-1} & \langle {x}^{2}{\rangle }_{i-1}\\ 1 & \langle x{\rangle }_{i} & \langle {x}^{2}{\rangle }_{i}\\ 1 & \langle x{\rangle }_{i+1} & \langle {x}^{2}{\rangle }_{i+1}\end{array}\right]\left[\begin{array}{l}{c}_{0}\\ {c}_{1}\\ {c}_{2}\end{array}\right]=\left[\begin{array}{l}{p}_{i-1}\\ {p}_{i}\\ {p}_{i+1}\end{array}\right],\end{eqnarray} \tag{ 15 }$$

which is then solved for the unique interpolation coefficients c₀, c₁, and c₂. The resulting polynomial is then used to determine ${p}_{R,i-1/2}\equiv p({x}_{i-1/2})$, the right state at the lower face of cell i, and ${p}_{L,i+1/2}\equiv p({x}_{i+1/2})$, the left state at the upper face of cell i. If ${p}_{R,i-1/2}$ is sufficiently close to either p_i−1 or p_i such that p(x) has an extremum for $x\in ({x}_{i-1},{x}_{i})$, then ${p}_{R,i-1/2}$ is reset to ensure that p(x) is monotone on this interval and has zero slope at x_i−1. The value of ${p}_{L,i+1/2}$ is reset in an analogous manner as needed. Next, p(x) is reconstructed once more according to the solution of the linear system

$$\begin{eqnarray}\left[\begin{array}{lll}1 & {x}_{i-1/2} & {x}_{i-1/2}^{2}\\ 1 & \langle x{\rangle }_{i} & \langle {x}^{2}{\rangle }_{i}\\ 1 & {x}_{i+1/2} & {x}_{i+1/2}^{2}\end{array}\right]\left[\begin{array}{l}{c}_{0}\\ {c}_{1}\\ {c}_{2}\end{array}\right]=\left[\begin{array}{l}{p}_{R,i-1/2}\\ {p}_{i}\\ {p}_{L,i+1/2}\end{array}\right].\end{eqnarray} \tag{ 16 }$$

If neither state has been reset by prior monotonicity constraints, then the solutions of Equations (15) and (16) are identical. Finally, if p(x) has an extremum in cell i, then either ${p}_{R,i-1/2}$ or ${p}_{L,i+1/2}$—whichever state is nearest to the extremum—is reset to ensure that p(x) is monotone in cell i with zero slope at that face.

For smooth data away from extrema, this process yields a third-order-accurate piecewise-parabolic reconstruction at the lower and upper faces of each cell. Once the reconstruction step is completed in all cells and for all variables, the resulting left and right states at each face, which are derived from distinct interpolation parabolas, are then passed to a Riemann solver to obtain the numerical flux at that face as described in the next section.

Fornax currently implements three choices for computing fluxes as a function of the reconstructed left and right states: the local Lax–Friedrichs, HLLE, and HLLC solvers. The HLLC solver, incorporating the most complete information on the modal structure of the equations, is the least diffusive option and is the default choice. Unfortunately, so-called three-wave solvers are susceptible to the carbuncle or odd–even instability. Many schemes have been proposed to inhibit the development of this numerical instability. In Fornax, we tag interfaces determined to be within shocks (by the pressure jump across the interface) and switch to the HLLE solver in orthogonal directions for cells adjacent to the tagged interface. In all our tests, this very simple approach has been successful at preventing the instability, incurs essentially no cost, and is typically used so sparingly that the additional diffusion introduced into the solution is likely negligible.

To ${ \mathcal O }({\rm{\Delta }}{x}^{2})$, we compute face-averaged fluxes between cells using approximate Riemann solvers applied to reconstructed states at the area centroid of each face. For the hydrodynamic fluxes, we use the HLLC solver of Batten et al. (1997) and Toro et al. (1994). For the fluxes of the radiation moments (specifically, the Fⁱ and ${c}^{2}{P}_{j}^{i}$ terms), we use the HLL solver as in Harten et al. (1983). As has been pointed out by multiple authors, the HLL solver applied to our radiation subsystem yields fluxes that fail to preserve the asymptotic diffusion limit. To recover this limit, we have found it sufficient to introduce a correction to the energy fluxes of the form

$$\begin{eqnarray}&&{{ \mathcal F }}_{E,\mathrm{corrected}}^{\mathrm{HLL}}=\displaystyle \frac{{S}_{R}{F}_{L}-{S}_{L}{F}_{R}+\epsilon {S}_{R}{S}_{L}({E}_{R}-{E}_{L})}{{S}_{R}-{S}_{L}},\end{eqnarray} \tag{ 17 }$$

with

$$\begin{eqnarray}&&\epsilon \equiv \min \left(1,\displaystyle \frac{1}{{\tau }_{\mathrm{cell}}}\right).\end{eqnarray} \tag{ 18 }$$

In Equation (18), τ_cell is an approximation to the optical depth of a cell in the direction normal to the face, computed using a simple arithmetic average of the total opacity (absorption plus scattering) of the cells on either side of the face. In Equation (17), S_L and S_R are the wavespeed estimates for the fastest left- and right-going waves, F_L and F_R are the normal components of the reconstructed radiation fluxes on the left and right sides of the face, and E_L and E_R are the reconstructed radiation energy densities on the left and right sides of the face. Note that when τ_cell ≤ 1, the corrected numerical fluxes are identical to the standard HLL fluxes. Note further that from here on, we drop indices referring to radiation species (photon or neutrino) and frequency group for clarity, though each of these radiation terms are computed per species and per group.

Since we evolve the comoving-frame radiation moments, the intercell fluxes include advective terms that depend on the fluid velocity component normal to the faces. For consistency and accuracy, we make use of the contact (i.e., particle) velocity (v*) computed during the construction of the hydrodynamic fluxes in the HLLC Riemann solver. At a given cell face, these advective fluxes are then set to

$$\begin{eqnarray}({vE})=\left\{\begin{array}{cc}{v}^{* }{E}_{R}, & {v}^{* }\lt 0\\ {v}^{* }{E}_{L}, & {v}^{* }\geqslant 0\end{array}\right.,\end{eqnarray} \tag{ 19a }$$

$$\begin{eqnarray}(v{\boldsymbol{F}})=\left\{\begin{array}{cc}{v}^{* }{{\boldsymbol{F}}}_{R}, & {v}^{* }\lt 0\\ {v}^{* }{{\boldsymbol{F}}}_{L}, & {v}^{* }\geqslant 0\end{array}\right.,\end{eqnarray} \tag{ 19b }$$

where the subscripts L and R indicate quantities reconstructed to the left or right side of the cell face, respectively. In other words, we use the reconstructed quantities in the upwind direction, defined with respect to v*, the particle velocity at a given cell face computed by the Riemann solver for the hydrodynamics.

Several terms in the radiation moment equations depend upon the gradient of the velocity field. Ultimately, volume averages of each of these velocity-gradient-dependent terms are needed. These are approximated, as in the calculation of the geometric source terms, by combining separately volume-averaged components as

$$\begin{eqnarray}&&\langle {{v}^{i}}_{;j}\rangle \approx \langle {{v}^{i}}_{,j}\rangle +\langle {{\rm{\Gamma }}}_{{jk}}^{i}\rangle \langle {v}^{k}\rangle .\end{eqnarray} \tag{ 20 }$$

As before with the advective fluxes, we again make use of the particle velocities computed during the Riemann solve for the hydrodynamic fluxes. In addition to the normal components, we also need the transverse components; these are available in the construction of the hydrodynamic fluxes as well. Thus, given the velocity components on every face of a cell, we assume a linear profile (in our generic coordinates xⁱ) for the intracell velocities and compute the volume-averaged partial derivatives directly via finite differences across the cell. Finally, to construct the volume-averaged velocities $\langle {v}^{k}\rangle$, we compute the volume-average of the linear profile in each direction k.

For the terms involving frequency-space derivatives, we use an approach similar to Vaytet et al. (2011), considering our multigroup treatment as a generalized finite-volume discretization of frequency space. The evolution equations for these terms can be written as

$$\begin{eqnarray}&&{E}_{g,t}-{{v}^{i}}_{;j}\left[{(\nu {P}_{\nu i}^{j})}_{g+1/2}-{(\nu {P}_{\nu i}^{j})}_{g-1/2}\right]=0,\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&{F}_{{gj},t}-{{v}^{i}}_{;k}\left[{(\nu {Q}_{\nu {ji}}^{k})}_{g+1/2}-{(\nu {Q}_{\nu {ji}}^{k})}_{g-1/2}\right]=0.\end{eqnarray} \tag{ 21b }$$

As in Vaytet et al. (2011), we adopt a simple upwind formulation of these terms, defining the intergroup flux at a given group boundary ${\nu }_{g-1/2}$ as

$$\begin{eqnarray}{(\nu {P}_{\nu i}^{j})}_{g-1/2}=\left\{\begin{array}{cc}{\nu }_{g-1/2}\,\langle {P}_{\nu i}^{j}{\rangle }_{L}, & {{v}^{i}}_{;j}\lt 0\\ {\nu }_{g-1/2}\,\langle {P}_{\nu i}^{j}{\rangle }_{R}, & {{v}^{i}}_{;j}\geqslant 0\end{array}\right.,\end{eqnarray} \tag{ 22a }$$

$$\begin{eqnarray}{(\nu {Q}_{\nu {ji}}^{k})}_{g-1/2}=\left\{\begin{array}{cc}{\nu }_{g-1/2}\,\langle {Q}_{\nu {ji}}^{k}{\rangle }_{L}, & {{v}^{i}}_{;k}\lt 0\\ {\nu }_{g-1/2}\,\langle {Q}_{\nu {ji}}^{k}{\rangle }_{R}, & {{v}^{i}}_{;k}\geqslant 0\end{array}\right.,\end{eqnarray} \tag{ 22b }$$

where for a given frequency group g, $\langle {P}_{\nu i}^{j}{\rangle }_{g}={P}_{{gi}}^{j}/{\rm{\Delta }}{\nu }_{g}$ is the group-averaged spectral density, Δν_g is the group width, and the subscripts L and R indicate the reconstruction of the spectral densities to the left or right side of the group boundary, respectively.

The interaction terms that couple the radiation and matter can introduce significant stiffness into the system, requiring an implicit treatment to provide stability to their numerical evolution on hydrodynamic timescales. In Fornax, we use operator splitting to advance the system:

$$\begin{eqnarray}&&{(\rho {v}_{j})}_{,t}={S}_{j},\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}&&{(\rho e)}_{,t}={S}_{E},\end{eqnarray} \tag{ 23b }$$

$$\begin{eqnarray}&&{(\rho X)}_{,t}={S}_{X},\end{eqnarray} \tag{ 23c }$$

$$\begin{eqnarray}&&{E}_{{sg},t}={R}_{{sgE}},\end{eqnarray} \tag{ 23d }$$

$$\begin{eqnarray}&&{F}_{{sgj},t}={R}_{{sgj}},\end{eqnarray} \tag{ 23e }$$

after the explicit transport update described above is applied. It is essential to treat these terms after the explicit update to capture stiff equilibria properly, as obtained in optically thick environments.

Equations 23(b)–(d) treat the terms describing the emission and absorption of radiation and their coupling to the material energy and composition. The essential task is to solve the system

$$\begin{eqnarray}&&\displaystyle \frac{{u}^{n+1}-{u}^{-}}{{\rm{\Delta }}t}=-\displaystyle \sum _{s}\displaystyle \sum _{g}({j}_{{sg}}^{n+1}-c{\kappa }_{{sg}}^{n+1}{E}_{{sg}}^{n+1}),\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{X}^{n+1}-{X}^{-}}{{\rm{\Delta }}t}=\displaystyle \sum _{s}\displaystyle \sum _{g}{\xi }_{{sg}}({j}_{{sg}}^{n+1}-c{\kappa }_{{sg}}^{n+1}{E}_{{sg}}^{n+1}),\end{eqnarray} \tag{ 24b }$$

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{{sg}}^{n+1}-{E}_{{sg}}^{-}}{{\rm{\Delta }}t}={j}_{{sg}}^{n+1}-c{\kappa }_{{sg}}^{n+1}{E}_{{sg}}^{n+1},\end{eqnarray} \tag{ 24c }$$

where u denotes the material internal energy density and the "–" superscript denotes the state just after the explicit update described above. Although Equations (24) are spatially decoupled, i.e., their solution depends entirely on local data, if solved directly, they would represent a large system of stiff nonlinear equations, especially in the multispecies, multigroup context. Fortunately, the solution of this system can be made simpler by separating the updates into "inner" and "outer" parts. In the inner update, the frequency groups are decoupled from one another and can be updated using a direct implicit scheme. Then, in the outer update, all that remains is to find the roots of just a few equations, independent of the number of groups or species. This can be accomplished using an iterative scheme. In the photon radiation context, there is just a single equation, and in the CCSN context, there are two. Without loss of generality, we continue our description for the core-collapse case, where X represents the electron fraction, Y_e, and the opacities and emissivities depend on both T and Y_e.

In the inner update, each individual frequency group is updated using a fully implicit backward Euler scheme of the form

$$\begin{eqnarray}&&{E}_{{sg}}^{k}=\displaystyle \frac{{E}_{{sg}}^{-}+{\rm{\Delta }}t\,{j}_{{sg}}^{k}}{1+c\,{\rm{\Delta }}t\,{\kappa }_{{sg}}^{k}},\end{eqnarray} \tag{ 25 }$$

where the opacities and emissivities are computed using the values of T^k and ${Y}_{e}^{k}$ at outer iteration k. Equation (25) has the important property wherein for $c\,{\rm{\Delta }}t\,{\kappa }_{{sg}}^{k}\gg 1$, i.e., for large optical depths, the updated energy, ${E}_{{sg}}^{k}$, approaches the thermal equilibrium solution, ${j}_{{sg}}^{k}/(c\,{\kappa }_{{sg}}^{k})$.

In the outer update, the hydrodynamic quantities u and ρY_e are updated implicitly according to the sum over groups (and of species) of the emission and absorption terms. The integrated changes in energy and composition are given by

$$\begin{eqnarray}&&{\rm{\Delta }}{E}^{k}=\displaystyle \sum _{s}\displaystyle \sum _{g}({E}_{{sg}}^{k}-{E}_{{sg}}^{-}),\end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{X}^{k}=\displaystyle \sum _{s}\displaystyle \sum _{g}{\xi }_{{sg}}({E}_{{sg}}^{k}-{E}_{{sg}}^{-}).\end{eqnarray} \tag{ 26b }$$

Finally, we compute the residuals

$$\begin{eqnarray}&&{r}_{E}^{k}={u}^{k}-{u}^{-}+{\rm{\Delta }}{E}^{k},\end{eqnarray} \tag{ 27a }$$

$$\begin{eqnarray}&&{r}_{X}^{k}=\rho ({Y}_{e}^{k}-{Y}_{e}^{-})-{\rm{\Delta }}{X}^{k},\end{eqnarray} \tag{ 27b }$$

where, like the opacities and emissivities, the internal energy density u^k depends on both T^k and ${Y}_{e}^{k}$. The solution to our implicit system is given by the values of T and Y_e for which ${r}_{E}^{k}={r}_{X}^{k}=0$.⁸ To find these roots, we use a Newton–Raphson iteration, using finite differences to form the Jacobian and employing a backtracking line search algorithm to improve robustness.⁹ At each iteration, a linear system is solved for $\delta {T}^{k+1}\equiv {T}^{k+1}-{T}^{k}$ and $\delta {Ye}\equiv {Y}_{e}^{k+1}-{Y}_{e}^{k}$, the changes in temperature and electron fraction, respectively, between iterations k and k + 1. It is in solving this linear system that our reduction to two equations has obvious virtues. Convergence is checked by requiring that the relative change of the temperature and composition between iterations is below a user-specified tolerance, typically 10⁻⁶. In practice, this procedure yields converged solutions within a few iterations under a wide range of conditions.

Finally, Equations 23(a), (b), and (e) treat the momentum coupling represented by the absorption and scattering of radiative flux. For each frequency group, we first update the jth component of the flux via the implicit backward Euler scheme,

$$\begin{eqnarray}&&{F}_{{sgj}}^{n+1}=\displaystyle \frac{{F}_{{sgj}}^{-}}{1+c\,{\rm{\Delta }}t({\kappa }_{{sg}}^{n+1}+{\sigma }_{{sg}}^{n+1})},\end{eqnarray} \tag{ 28 }$$

which is crucial for maintaining the correct asymptotic behavior in the diffusion regime. To see this, consider the relevant case of a near-equilibrium, static atmosphere with large cell optical depths, as is approximately the case inside the protoneutron star in the core-collapse problem. Under these conditions, $\parallel F\parallel /({cE})\ll 1$, so that the pressure tensor under the M1 closure approaches the Eddington closure, ${{P}^{i}}_{j}=(E/3){{\delta }^{i}}_{j}$. The flux after the explicit transport update described above is then given by

$$\begin{eqnarray}&&{F}_{{sgj}}^{-}={F}_{{sgj}}^{n}-{\rm{\Delta }}t\displaystyle \frac{{c}^{2}}{3}{E}_{{sg},j}^{n}\approx -{\rm{\Delta }}t\displaystyle \frac{{c}^{2}}{3}{E}_{{sg},j}^{n},\end{eqnarray} \tag{ 29 }$$

where the last approximate equality holds in a near-equilibrium, optically thick atmosphere in which $F\sim {cE}/\tau \ll {c}^{2}{E}_{,i}$. In this case, our implicit update yields the updated flux

$$\begin{eqnarray}\begin{array}{rcl}{F}_{{sgj}}^{n+1} & = & \displaystyle \frac{{F}_{{sgj}}^{-}}{1+c{\rm{\Delta }}t({\kappa }_{{sg}}^{n+1}+{\sigma }_{{sg}}^{n+1})},\\ & & \approx -\left[\displaystyle \frac{c}{3({\kappa }_{{sg}}^{n+1}+{\sigma }_{{sg}}^{n+1})}\right]{E}_{{sg},j}^{n},\end{array}\end{eqnarray} \tag{ 30 }$$

which is a temporally first-order accurate expression of Fick's law of diffusion. In other words, our operator-split implicit update of the flux reduces to Fick's law in the asymptotic limit of high optical depth. The corresponding explicit material momentum and energy updates  are then given by

$$\begin{eqnarray}&&{\left(\rho {v}_{j}\right)}^{n+1}={\left(\rho {v}_{j}\right)}^{-}-\displaystyle \frac{1}{{c}^{2}}\displaystyle \sum _{s}\displaystyle \sum _{g}\left({F}_{{sgj}}^{n+1}-{F}_{{sgj}}^{-}\right),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{\left(\rho e\right)}^{n+1}={\left(\rho e\right)}^{-}-\displaystyle \frac{1}{{c}^{2}}\displaystyle \sum _{s}\displaystyle \sum _{g}{\left({v}^{j}\right)}^{-}\left({F}_{{sgj}}^{n+1}-{F}_{{sgj}}^{-}\right).\end{eqnarray} \tag{ 32 }$$

Note that Equations (31) and (32) are not momentum- and energy-conservative, respectively, because the equations themselves are not momentum- and energy-conservative as written in the comoving frame. In Section 9.1, we demonstrate how well Fornax conserves total energy in a realistic core-collapse problem for which the code was primarily intended. This completes the description of our treatment of the interaction terms.

Advance the explicit transport subsystem (including advective flux terms and geometrical and gravitational source terms) by Δt:

$$\begin{eqnarray}&&{\rho }_{,t}+{(\rho {v}^{i})}_{;i}=0,\end{eqnarray} \tag{ 33a }$$

$$\begin{eqnarray}&&{(\rho {v}_{j})}_{,t}+{(\rho {v}^{i}{v}_{j}+P{{\delta }^{i}}_{j})}_{;i}=-\rho {\phi }_{,j},\end{eqnarray} \tag{ 33b }$$

$$\begin{eqnarray}&&{(\rho e)}_{,t}+{\left[\rho {v}^{i}\left(e+\displaystyle \frac{P}{\rho }\right)\right]}_{;i}=-\rho {v}^{i}{\phi }_{,i},\end{eqnarray} \tag{ 33c }$$

$$\begin{eqnarray}&&{(\rho {Y}_{e})}_{,t}+{(\rho {Y}_{e}{v}^{i})}_{;i}=0,\end{eqnarray} \tag{ 33d }$$

$$\begin{eqnarray}&&{E}_{s\nu ,t}+{({F}_{s\nu }^{i}+{v}^{i}{E}_{s\nu })}_{;i}=0,\end{eqnarray} \tag{ 33e }$$

$$\begin{eqnarray}&&\begin{array}{l}{F}_{s\nu j,t}+{({c}^{2}{P}_{s\nu j}^{i}+{v}^{i}{F}_{s\nu j})}_{;i}\\ +\ {v}_{;j}^{i}{F}_{s\nu i}-{v}_{;k}^{i}{Q}_{s\nu {ji}}^{k}=0.\end{array}\end{eqnarray} \tag{ 33f }$$

Advance the frequency advection subsystem by Δt:

$$\begin{eqnarray}&&{E}_{s\nu ,t}+{v}_{;j}^{i}\left[{P}_{s\nu i}^{j}-{\partial }_{\nu }\left(\nu \,{P}_{s\nu i}^{j}\right)\right]=0,\end{eqnarray} \tag{ 34a }$$

$$\begin{eqnarray}&&{F}_{s\nu j,t}-{v}_{;k}^{i}\,{\partial }_{\nu }\left(\nu \,{Q}_{s\nu {ji}}^{k}\right)=0.\end{eqnarray} \tag{ 34b }$$

Advance the radiation–matter energy interaction subsystem by Δt:

$$\begin{eqnarray}&&{(\rho e)}_{,t}=-\displaystyle \sum _{s}{\int }_{0}^{\infty }\left({j}_{s\nu }-c{\kappa }_{s\nu }{E}_{s\nu }\right)\,d\nu ,\end{eqnarray} \tag{ 35a }$$

$$\begin{eqnarray}&&{(\rho {Y}_{e})}_{,t}=\displaystyle \sum _{s}{\int }_{0}^{\infty }{\xi }_{s\nu }({j}_{s\nu }-c{\kappa }_{s\nu }{E}_{s\nu })\,d\nu ,\end{eqnarray} \tag{ 35b }$$

$$\begin{eqnarray}&&{E}_{s\varepsilon ,t}={j}_{s\nu }-c{\kappa }_{s\nu }{E}_{s\nu }.\end{eqnarray} \tag{ 35c }$$

Advance the radiation–matter momentum interaction subsystem by Δt:

$$\begin{eqnarray}&&{(\rho {v}_{j})}_{,t}=\displaystyle \frac{1}{c}\displaystyle \sum _{s}{\int }_{0}^{\infty }({\kappa }_{s\nu }+{\sigma }_{s\nu }){F}_{s\nu j}\,d\nu ,\end{eqnarray} \tag{ 36a }$$

$$\begin{eqnarray}&&{(\rho e)}_{,t}=\displaystyle \frac{{v}^{i}}{c}\displaystyle \sum _{s}{\int }_{0}^{\infty }({\kappa }_{s\nu }+{\sigma }_{s\nu }){F}_{s\nu i}\,d\nu ,\end{eqnarray} \tag{ 36b }$$

$$\begin{eqnarray}&&{F}_{s\nu j,t}=-c({\kappa }_{s\nu }+{\sigma }_{s\nu }){F}_{s\nu j}.\end{eqnarray} \tag{ 36c }$$

The purpose of the spherical dendritic grid is to avoid the narrowing of zones in polar and azimuthal angular extent approaching the origin and the narrowing of zones in the azimuthal extent approaching the polar axes that is characteristic of the spherical polar coordinate mapping of a logically Cartesian grid. For example, a unigrid in ordinary spherical polar coordinates would have zones of size ${\rm{\Delta }}r\times r\,{\rm{\Delta }}\theta \times r\,\sin \theta \,{\rm{\Delta }}\phi$. Even at moderate angular resolution, this narrowing of zones can place a severe constraint on the CFL time step. One approach to circumvent this constraint is to coarsen the angular resolution as needed in order to maintain a roughly constant zone aspect ratio. Thus, we coarsen the polar angle in order to keep ${\rm{\Delta }}\theta \sim {\rm{\Delta }}r/r$ as $r\to 0$, and similarly, we coarsen the azimuthal angle in order to keep ${\rm{\Delta }}\phi \sim {\rm{\Delta }}r/r\,\sin \theta$ as $\theta \to 0$ or $\theta \to \pi$. As a result, all zones on the grid have roughly equal aspect ratios, and near the origin where Δr is approximately constant, all zones have roughly the same size. This implies that the CFL time step is essentially constrained only by the minimum radial resolution, Δr_min, i.e., the time step is not constrained by the resolution in either angular direction. A sample arrangement of the dendritic grid in 3D is given in Figure 1.

Zoom In Zoom Out Reset image size

In Fornax, angular resolution is always coarsened by a factor of 2. Thus, wherever the zone aspect ratio would otherwise deviate too far from 1, either Δθ or Δϕ is coarsened as needed,¹⁰ and the aspect ratio of the neighboring zone changes discontinuously by a factor of ∼2. It is optimal to coarsen in angle whenever a zone's aspect ratio would fall outside the range $[1/\sqrt{2},\sqrt{2}]$, ensuring that any zone is at most ∼41% longer in any direction than in the others. It follows that the CFL time-step constraint may be given by

$$\begin{eqnarray}&&{\rm{\Delta }}t\leqslant \displaystyle \frac{{\rm{\Delta }}{r}_{\min }}{3\sqrt{2}{\max }_{i}\{v+{c}_{s}\}},\end{eqnarray} \tag{ 37 }$$

where the max is taken over all zones i. Note that the notion of coarsening does not here imply that the overall zone size is increasing. On the contrary, in the portion of the grid where Δr is constant, the zones are all roughly of equal size. Therein, the resolution can be controlled simply by the value of Δr_min.

To continue, let ${{\ell }}_{\theta }$ and ${{\ell }}_{\phi }$ denote the refinement levels in the θ- and ϕ-directions, respectively, where the level ${{\ell }}_{d}=0$ denotes the coarsest resolution level in the d-direction. We require a minimum of two zones in the θ-direction at the origin and a minimum of four zones in the ϕ-direction at the polar axes, hence, the finest levels in each of these directions are given by ${{\ell }}_{\theta }={\mathrm{log}}_{2}{N}_{\theta }-1$ and ${{\ell }}_{\phi }={\mathrm{log}}_{2}{N}_{\phi }-2$, where N_θ and N_ϕ are the number of zones on the finest level in the θ- and ϕ-directions, respectively. For example, with a resolution of ${N}_{\theta }\times {N}_{\phi }=256\times 512$ zones, we would have maximum levels ${{\ell }}_{\theta }={{\ell }}_{\phi }=7$.

The indexing scheme can be rather subtle. A given zone with global coordinates (i, j, k) at levels ${{\ell }}_{\theta }(i)$ and ${{\ell }}_{\phi }(i,j)$ has a neighbor (or neighbors) in the r-direction at coordinates (i', j', k') at levels ℓ_θ(i') and ${{\ell }}_{\phi }(i^{\prime} ,j^{\prime} )$. The global coordinates are related by $j^{\prime} =\lfloor {n}_{\theta }j\rfloor$ and $k^{\prime} =\lfloor {n}_{\phi }k\rfloor$, where ${n}_{\theta }(i,i^{\prime} )\equiv {2}^{{{\ell }}_{\theta }(i^{\prime} )-{{\ell }}_{\theta }(i)}$ and ${n}_{\phi }(i,i^{\prime} ,j,j^{\prime} )\equiv {2}^{{{\ell }}_{\phi }(i^{\prime} ,j^{\prime} )-{{\ell }}_{\phi }(i,j)}$, and $\lfloor \cdot \rfloor$ denotes the integer floor operation.

The general approach of Fornax toward reconstruction is described in Section 5. However, with our dendritic grid, there are numerous special considerations. When performing reconstruction on the spherical dendritic grid in the r-direction in 2D (and in both the r- and θ-directions in 3D), we must take special care at refinement boundaries to use a pencil of zones of constant coordinate volume. For a given zone at refinement level ℓ, if its neighbor is at refinement level ${{\ell }}_{d}^{{\prime} }\gt {{\ell }}_{d}$ in a given angular direction d, then we perform a fine-to-coarse restriction of the neighbor data to level ℓ_d using an appropriate volume-weighted average. On the contrary, if its neighbor is at level ${{\ell }}_{d}^{{\prime} }\lt {{\ell }}_{d}$, then we perform a coarse-to-fine prolongation of the data to level ℓ_d using an interpolation of the neighbor data based on the monotonized central difference (MC)-limited linear slopes in the transverse direction(s).

The restriction from level ${{\ell }}_{d}^{{\prime} }$ down to level ℓ_d in the r-direction is then given by

$$\begin{eqnarray}&&{{\boldsymbol{U}}}_{{ijk}}=\displaystyle \frac{{\displaystyle \sum }_{p=j^{\prime} }^{j^{\prime} +{n}_{\theta }}{\displaystyle \sum }_{q=k^{\prime} }^{k^{\prime} +{n}_{\phi }}{\rm{\Delta }}{V}_{{ipq}}{{\boldsymbol{U}}}_{{ipq}}}{{\displaystyle \sum }_{p=j^{\prime} }^{j^{\prime} +{n}_{\theta }}{\displaystyle \sum }_{q=k^{\prime} }^{k^{\prime} +{n}_{\phi }}{\rm{\Delta }}{V}_{{ipq}}}.\end{eqnarray} \tag{ 38 }$$

The analogous prolongation from level ℓ_d up to level ${{\ell }}_{d}^{{\prime} }$ in the r-direction can be messy, but straightforward. Prolongation and restriction are done in order to obtain a pencil of constant coordinate volume, which is required to perform a volume-coordinate-based reconstruction.

Fornax is parallelized with non-blocking point-to-point communication and global all-to-all reduction using an implementation of the MPI-3.0 standard. Non-blocking communication allows us to hide most—if not all—of the communication latency behind a significant amount of local computation. Due to the structure of the dendritic grid, domain decomposition becomes a non-trivial issue. Unlike the case for a unigrid or for patch-based AMR/SMR, there is no known optimal domain decomposition vis á vis load balancing. Therefore, we rely on a type of greedy algorithm to map the grid onto processors. In our decomposition algorithm, we attempt to drive the maximum processor load as close to the average load as possible, within reason. At the same time, we attempt to minimize the total communication overhead required to exchange ghost zone data among neighboring grids. Since our decomposition algorithm scales as a power law in the number of processors, N_proc, it can take minutes to hours to compute a single decomposition for a 3D grid with a given resolution and spatial extent. However, we need only to compute this once; the resulting decomposition can be saved to a file and read in at startup.

We typically require the maximum processor load be no more than 50% larger than the average load. This is not possible for every conceivable N_proc, but for a given target, there is typically some nearby N_proc for which our algorithm yields sufficient load-balancing efficiency. We need only search within some range of N_proc around the target—a process that can be performed in parallel—and select the value of N_proc giving the most evenly balanced decomposition. Although our algorithm is likely suboptimal, we still achieve near-perfect strong scaling efficiency for ${N}_{\mathrm{proc}}\gtrsim \mathrm{few}\times {10}^{5}$ (Appendix E). This is due to the fact that the majority of zones belong to a single logically Cartesian block that can be optimally decomposed in a trivial manner. All that remains is to decompose the coarsened dendritic portion of the grid in a way that optimizes the overall load and communication overhead.

In spherical coordinates, certain unit vectors change sign discontinuously across coordinate singularities, i.e., the $\hat{r}$ and $\hat{\phi }$ unit vectors change sign across the origin in the radial direction, and the $\hat{\theta }$ and $\hat{\phi }$ unit vectors change sign across the axis in the polar direction. Following Baumgarte et al. (2013), we enforce parity conditions on the signs of the corresponding vector components in these ghost zones after copying the data from active zones. In the azimuthal direction, the grid at ϕ = 2π is mapped periodically to the grid at ϕ = 0, so that the data in the active zone (i, j, N_ϕ−1) are mapped to ghost zone (i, j, −1) and, similarly, the data in the active zone (i, j, 0) are mapped to ghost zone (i, j, N_ϕ). In the polar direction at the axes, the data are periodically shifted in the azimuthal direction by N_ϕ/2 zones as they are copied from the active zones into the ghost zones. Thus, the active zone (i, 0, k) is mapped to the ghost zone (i, −1, k'), where $k^{\prime} =(k+{N}_{\phi }/2)\mathrm{mod}\ {N}_{\phi }$. Finally, in the radial direction at the origin, the active-zone data are reflected in the polar direction and then periodically shifted as they are copied into the ghost zones. Thus, the active zone (0, j, k) is mapped to the ghost zone (−1, j', k), where j' = N_θ − 1 − j and k' is as before. At the outer radial boundary, Fornax can employ a variety of boundary conditions, including outflow conditions with or without a diode restriction on the mass or radiation flux, Dirichlet and Neumann conditions, and arbitrary user-specified conditions.

We have taken special care in treating the reconstruction of quantities on the mesh and have formulated our reconstruction methods to be agnostic to choices of coordinates and mesh equidistance. As in Mignone (2014), we do this by forming our reconstruction expressions in terms of volume-averaged quantities and moments of the coordinates in each cell. The moments are calculated at start up for the given mesh and coordinates. In the radial direction, we have a uniform mesh in coordinate x₁, where the normal spherical radius $r\,={r}_{t}\sinh ({x}_{1}/{r}_{t})$, with r_t a constant parameter. We reconstruct quantities in x₁, i.e., we locally compute f(x₁) to specify f on the faces of each cell. The mesh in x₁ is entirely uniform. We have no need of, and do not construct, f(r), where then the samples of f would be non-uniformly distributed in r. The same is true in the angular direction. Thus, we have no "non-equidistant" reconstruction, and therefore, our reconstructions yield the same results as Mignone (2014). For any choice of coordinates in any geometry, the method works the same and forms a consistent parabolic reconstruction from volume-averaged inputs.

To demonstrate the scaling of Fornax under this reconstruction regime, we have performed at various resolutions the L₁ error radial/spherical wind test discussed in Mignone (2014). This hydrodynamic problem has an analytic solution (Mignone, his Equation (89)), and we have conducted the test for the Gaussian density profile provided in his Equation (73), with his parameters (a = 10, b = 0) and (a = 16, b = 0.5). Figure 2 depicts the associated convergence scaling for the two parameter sets, compared with quadratic $\left(1/{N}^{2}\right)$ behavior, and should be compared with Mignone's Figure 10. Despite the fact that this study was performed for a non-uniform radial grid, the expected scaling with resolution still emerges.

Zoom In Zoom Out Reset image size

A standard test is the classical self-similar blast wave of Sedov (1959) describing, e.g., the initial energy-conserving phase of a supernova remnant. The problem consists of a stationary background medium with uniform density ${\rho }_{1}$ and near-zero pressure P₁ into which a prescribed blast energy E₁ is deposited. Since P₁ ≈ 0, the solution depends on only the two parameters ρ₁ and E₁; it follows that the radius, r, and the time since explosion, t, are related by a dimensionless similarity variable, ξ. Thus, it can be shown that the position of the shock radius is given by

$$\begin{eqnarray}&&{r}_{\mathrm{sh}}(t)={\xi }_{0}{\left(\displaystyle \frac{{E}_{1}{t}^{2}}{{\rho }_{1}}\right)}^{1/5},\end{eqnarray} \tag{ 39 }$$

where ξ₀ ≈ 1.033 for an adiabatic equation of state with γ = 1.4. Following Kamm & Timmes (2007), we set the dimensionless inputs to γ = 1.4, ρ₁ = 1, P₁ = 4 × 10⁻¹³, and E₁ = 0.851072, such that the outer blast radius is at ${r}_{\mathrm{sh}}=1$ at time t = 1. We use a three-dimensional spherical dendritic grid with radial extent out to ${r}_{\max }=1.2$ and a resolution of 64 × 64 × 128 zones. The radial grid spacing is constant in the interior with Δr ≈ 0.01 out to r ≈ 0.3, then smoothly transitions to logarithmic spacing out to r_max.

For our first test, we deposit all of E₁ in the zones surrounding the origin. As expected, the subsequent blast is aligned with the grid and maintains perfect spherical symmetry. Next, we deposit the same energy in the zones closest to the x-axis at r = 0.12 so that the blast is off-axis. This time, the blast is not aligned with the grid and must cross the polar axis as it expands. Figures 3 and 4 compare the density and internal energy of the on- and off-axis blast waves, respectively. There are small-scale axis artifacts in the density in the off-axis version of the blast wave, but the blast wave remains fairly spherical despite being nowhere aligned with the spherical grid. The scale of the artifacts is smaller in the internal energy, since pressure variations tend to be smoothed out behind the shock where the pressure is large.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall, this test demonstrates that the internal geometric boundary conditions at the origin and polar axes are correctly applied. The artifacts at the axes are likely due to small errors in reconstruction in the azimuthal direction with our dendritic grid. The direction of the polar and azimuthal coordinate vectors is divergent at the polar axis, hence angular variations in the data there can only be resolved by a correspondingly convergent grid. In the limit of a locally "planar" flow across the poles, the rapid azimuthal variation comes entirely from the coordinates themselves, not from the data, although these variations must average to zero. Our dendritic grid deresolves this rapid angular variation, introducing first-order errors due to monotonicity contraints. However, these artifacts are confined to a narrow polar axis region and do not adversely contaminate the solution elsewhere.

To demonstrate the shock-capturing capabilities of Fornax, we consider the classical shock-tube test of Sod (1978). The initial data consist of two states:

$$\begin{eqnarray}(\rho ,P)=\left\{\begin{array}{cc}(1.0,1.0), & {\rm{if}}\,x\lt 0,\\ (0.125,0.1), & {\rm{otherwise}}.\end{array}\right.\end{eqnarray} \tag{ 40 }$$

The initial velocity is set to zero. For this test, we adopt an adiabatic equation of state with index γ = 1.4.

This test is not particularly challenging for modern HRSC codes, but it is nevertheless interesting because it exhibits all fundamental hydrodynamical waves. The exact solution consists of a left propagating rarefaction wave, a right propagating contact discontinuity, and a right propagating shock.

Figure 5 shows the analytic solution and the results from Fornax using a coarse mesh of 64 points. The code correctly captured all of the waves. The shock wave is sharply captured within ∼2 grid points, while the contact wave is smeared over ∼4 grid points.

Zoom In Zoom Out Reset image size

To test the sensitivity of the Fornax code to the numerical diffusion of contact waves, we perform the classic double Mach reflection test of Woodward & Colella (1984). This test consists of a Mach 10 oblique shock through air (γ = 1.4), inclined with respect to a reflecting boundary in which the incident and reflected shocks then interact to produce a triple point. In the space between the reflected shock and the contact discontinuity, an upward-directed jet should form along the slip surface, but if the numerical dissipation of the contact wave is too high, the formation of this jet will be suppressed. We use a two-dimensional grid of 520 × 160 zones over the domain $(x,y)\in [0,3.25]\times [0,1]$ and position an oblique shock at x₀ = 1/6 inclined at an angle α = π/3 with respect to the x-axis. The pre-shock state is given by ${(\rho ,P,{v}_{x},{v}_{y})}_{{\rm{R}}}=(1.4,1,0,0)$ and the post-shock state by ${(\rho ,P,{v}_{x},{v}_{y})}_{{\rm{L}}}=(8,116.5,8.25\sin \alpha ,-8.25\cos \alpha$). The left x-boundary is held fixed at the post-shock state and outflow conditions are imposed at the right x-boundary. At the lower y-boundary, the post-shock state is fixed for x ≤ x₀, and reflecting boundary conditions are imposed for $x\gt {x}_{0}$. Meanwhile, at the upper y-boundary, the time-dependent shock position, ${x}_{s}\,={x}_{0}+y/\tan \alpha +10t/\sin \alpha$, is used to set either the post-shock state (x ≤ x_s) or the pre-shock state (x > x_s).

Figure 6 shows the evolved system at t_final = 0.2 using 30 linearly spaced contours of the density. The features at the various shock fronts remain sharp, and in the region below the triple point, a small upward-directed jet is indeed produced along the slip surface. This indicates that numerical dissipation of the contact wave is well controlled in Fornax.

Zoom In Zoom Out Reset image size

8.5.1. Two-dimensional Rayleigh–Taylor Test

Here, we study the linear and nonlinear development of the Rayleigh–Taylor instability in 2D. The initial data describe an unstably stratified fluid with

$$\begin{eqnarray}\rho =\left\{\begin{array}{cc}{\rho }_{h} & {\rm{if}}\,z\geqslant 0,\\ {\rho }_{l} & {\rm{otherwise}}.\end{array}\right.\end{eqnarray} \tag{ 41 }$$

The gravitational acceleration g = 1/2 is parallel to the z axis, and ρ_h = 2, ρ_l = 1, so that the Atwood number $A=({\rho }_{h}-{\rho }_{l})/({\rho }_{h}+{\rho }_{l})$ is 1/3. The boundary conditions are reflective for z = ±1 and periodic for x = ±1/2. The initial pressure is chosen to ensure hydrostatic equilibrium:

$$\begin{eqnarray}&&P(z)=P(-1)+{\int }_{-1}^{z}\rho \,g\,{dz},\end{eqnarray} \tag{ 42 }$$

with P(−1) = 10/7 + 1/4. For this test, we adopt an adiabatic equation of state with index γ = 1.4.

At time t = 0, the interface between the two fluid components is perturbed to be

$$\begin{eqnarray}&&h(x)={h}_{0}\cos (\kappa x),\end{eqnarray} \tag{ 43 }$$

with h₀ = 0.01 and $\kappa =4\pi$. The sharp interface is then smoothed into a hyperbolic tangent profile with characteristic length 0.005.

According to analytic theory, h should evolve according to the following equation:

$$\begin{eqnarray}&&h(x,t)={h}_{0}\,\cosh (\sqrt{A\kappa g}\,t)\cos (\kappa x).\end{eqnarray} \tag{ 44 }$$

We perform simulations with resolutions ranging from 64 ×128 to 1024 × 2048. We track the mixing of the two fluid components by means of a passive scalar tracer. This is initialized to be $\rho -{\rho }_{l}$, so that the value 1 corresponds to the heavy fluid, while the value 0 corresponds to the light fluid. We track the growth of the initial perturbation by locating the 0.99 isocontour of the tracer. Figure 7 shows the perturbation amplitude, normalized by its initial value, and the prediction from analytic theory (Equation (44)).

Zoom In Zoom Out Reset image size

At low resolution, the perturbations grow exponentially from t = 0, but the growth rate is somewhat smaller than that predicted by analytic theory. The slower growth rate at low resolution is not unexpected and is due to the numerical viscosity, which modifies the Rayleigh–Taylor dispersion relation (Dimonte et al. 2004; Murphy & Burrows 2008a). We speculate that the reason why the initial perturbations grow exponentially at low resolution is that they are not well resolved. At higher resolution, the correct cosh time dependence is recovered and the agreement with linear theory is excellent up to time t ≃ 1.75, when the Rayleigh–Taylor plumes start to break and rolls develop.

During the nonlinear phase of the evolution, secondary Rayleigh–Taylor and Kelvin–Helmholtz instabilities appear (Figure 8). These are seeded by the numerical noise, and their detailed morphology, in the absence of explicit dissipation, is known to be dependent on the details of the numerical scheme (Liska & Wendroff 2003a). Despite the presence of these features, Fornax is able to preserve the sharp discontinuity in the fluid tracer. Artificial mixing between the two fluid components is only present in the Rayleigh–Taylor rolls, where the flow develops features on scales comparable to those of the grid.

Zoom In Zoom Out Reset image size

8.5.2. Three-dimensional Rayleigh–Taylor Test

Neutrino-driven convection plays a central role in the explosion mechanism of CCSNe (Burrows 2013; Burrows et al. 2018; Radice et al. 2018; Vartanyan et al. 2018). In this section, we benchmark Fornax for the modeling of convective flows by studying the nonlinear development of the Rayleigh–Taylor instability in 3D. We consider the setup introduced by Dimonte et al. (2004), for which the dynamics is dominated by mode couplings.

The computational domain is a box with $-L/2\lt x,y\,\lt L/2$, $-L\lt z\lt L$, where L = 10 in arbitrary units. We consider a constant vertical gravitational acceleration $g\,:= -{g}_{z}=2$, and we prepare initial conditions with density stratification

$$\begin{eqnarray}&&\rho (z)={\rho }_{0}{\left(1-\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{{\rho }_{0}{gz}}{{P}_{0}}\right)}^{\tfrac{1}{\gamma -1}}\,,\end{eqnarray} \tag{ 45 }$$

where γ = 5/3,

$$\begin{eqnarray}{\rho }_{0}=\left\{\begin{array}{cc}{\rho }_{h} & {\rm{if}}\,z\geqslant 0,\\ {\rho }_{l} & {\rm{otherwise}},\end{array}\right.\end{eqnarray} \tag{ 46 }$$

and ${p}_{0}=2\pi ({\rho }_{h}+{\rho }_{l}){gL}$. We set ${\rho }_{h}=3$ and ${\rho }_{l}=1$, so that the Atwood number $A=({\rho }_{h}-{\rho }_{l})/({\rho }_{h}+{\rho }_{l})$ is 1/2. The initial pressure is set to

$$\begin{eqnarray}&&P={P}_{0}{\left(\displaystyle \frac{\rho }{{\rho }_{0}}\right)}^{\gamma },\end{eqnarray} \tag{ 47 }$$

which ensures that the initial conditions are in hydrostatic equilibrium. We track the concentration of the "heavy" fluid by evolving a passive scalar field f_h initialized to be one for positive z and zero otherwise. We assume periodicity in the x and y directions and hydrostatic boundary conditions at z = ±L. We use uniform grids and label our simulations by the number of grid points in the x-direction so that, e.g., the N128 run has a resolution of 128 × 128 × 256 points.

At time t = 0, the interface between the two fluid components is perturbed to be

$$\begin{eqnarray}&&\begin{array}{l}h(x,y)=\displaystyle \frac{1}{H}\displaystyle \sum _{8\leqslant {k}_{x}^{2}+{k}_{y}^{2}\leqslant 16}\left[{a}_{{\boldsymbol{k}}}\cos ({k}_{x}X)\cos ({k}_{y}Y)\right.\\ \quad +\ {b}_{{\boldsymbol{k}}}\cos ({k}_{x}X)\sin ({k}_{y}Y)+{c}_{{\boldsymbol{k}}}\sin ({k}_{x}X)\cos ({k}_{y}Y)\\ \quad \left.+\ {d}_{{\boldsymbol{k}}}\sin ({k}_{x}X)\sin ({k}_{y}Y)\right],\end{array}\end{eqnarray} \tag{ 48 }$$

where $X=2\pi x/L$ and $Y=2\pi y/L$. The coefficients ${a}_{{\boldsymbol{k}}}$, ${b}_{{\boldsymbol{k}}}$, ${c}_{{\boldsymbol{k}}}$, and ${d}_{{\boldsymbol{k}}}$ are sampled from a uniform distribution taking values between −1 and 1, boundary excluded. H is a normalization coefficient that we adjust after having sampled ${a}_{{\boldsymbol{k}}}$, ${b}_{{\boldsymbol{k}}}$, ${c}_{{\boldsymbol{k}}}$, and ${d}_{{\boldsymbol{k}}}$ so that ${h}_{\mathrm{rms}}=3\times {10}^{-4}L$. We remark that we use the same initial conditions for all of our simulations, while Dimonte et al. (2004) used a different number of modes in their initial perturbation depending on the resolution. We also use the same realization of the coefficients in Equation (48) for all the simulations presented in this section. In this way, the convergence of our numerical results can be better assessed.

When constructing the perturbed initial data, we compute the cell-averaged density, internal energy, and heavy-fluid concentration, i.e., the quantities evolved by the finite-volume scheme in Fornax, by numerically integrating the profiles discussed above on an auxiliary mesh with 50 times higher resolution than the base grid in each direction.

The initial perturbations generate small buoyant bubbles that interact and grow into larger bubbles as the system evolves, as shown in Figure 9. We find that there is significant mixing between the two fluids. This is evidenced by the intermediate values taken by $\langle {f}_{h}\rangle$ (see Equation (49) below), especially in the buoyant bubbles, which entrain a significant amount of the heavy fluid. At late times, the topology of the interface between the two fluids becomes non-trivial, and we observe the detachment of some of the bubbles.

Zoom In Zoom Out Reset image size

The entrainment of the heavy fluid in the buoyant bubbles and of the light fluid in the sinking plumes is not strongly dependent on the resolution, but they appear to grow slightly with resolution, as shown in Figure 10. This is presumably due to the higher effective Reynolds numbers achieved in the best resolved simulations. Note that we do not include an explicit viscosity in these simulations; consequently, the dissipation scale is a multiple of the grid scale. At high resolutions, secondary instabilities with smaller wavenumbers and faster growth rates are allowed to develop. Their presence might explain the increase in the entrainment with resolution. Despite these quantitative differences, we find a remarkable qualitative agreement between the different resolutions with the gross flow feature captured even at the lowest resolution of 64 × 64 ×128.

Zoom In Zoom Out Reset image size

Following Dimonte et al. (2004), we compute the horizontally averaged concentration of the heavy fluid,

$$\begin{eqnarray}&&\langle {f}_{h}\rangle =\displaystyle \frac{1}{{L}^{2}}\iint {f}_{h}\,{dx}\,{dy},\end{eqnarray} \tag{ 49 }$$

and we define as the bubble penetration depth h_b the z-value at which $\langle {f}_{h}\rangle$ reaches 99%. The profile of $\langle {f}_{h}\rangle$ is observed to be self-similar during the nonlinear development of the Rayleigh–Taylor instability. We show the profile of $\langle {f}_{h}\rangle$ obtained with Fornax at two representative times in Figure 11. Also shown are the experimental results obtained by Dimonte et al. (2004), as reported in that paper. This figure should be contrasted with Figure 7 of Dimonte et al. (2004). We find good agreement between the Fornax results and the experimental data. The agreement improves as the resolution increases, with the exception of a spike appearing at late times in the N512 data around $z/{h}_{b}\simeq -1$. This is caused by the interaction of some of the down-falling plumes with the boundary at the bottom of the simulation box.

Zoom In Zoom Out Reset image size

Following Andrews & Spalding (1990) and Dimonte et al. (2004), we quantify the mixing degree between the heavy and the light fluid with the metric

$$\begin{eqnarray}&&W=\int \langle {f}_{h}\rangle \left(1-\langle {f}_{h}\rangle \right){dz}.\end{eqnarray} \tag{ 50 }$$

Our results are shown in Figure 12. The growth rate of the instability in the nonlinear phase is proportional to ${{Agt}}^{2}/L$, as predicted by theoretical models and confirmed by experimental results (Dimonte et al. 2004). The progression of the mixing between the two fluids as computed by Fornax agrees well with the experimental results reported in Dimonte et al. (2004). We find that Fornax compares favorably to the other Eulerian codes presented in Dimonte et al. (2004), such as FLASH. Fornax yields results comparable in quality to those of arbitrary Lagrangian–Eulerian codes with interface reconstruction methods (see Figure 8 of Dimonte et al. 2004). That said, we want to remark that, as documented in Dimonte et al. (2004) and as we have confirmed in preliminary calculations, the exact rate of mixing between the two fluids depends to some extent on the spectrum and the character of the initial perturbations. We have not fine-tuned the initial conditions of our simulations to match the experimental data.¹¹ Nevertheless, we cannot exclude the possibility that the better agreement with the experimental data of Fornax compared to the results from the other Eulerian codes presented in Dimonte et al. (2004) is due partly to the different choice of initial conditions.

Zoom In Zoom Out Reset image size

We analyze the character of the fluctuations in f_h on the original interface at z = 0 at an advanced time ${{Agt}}^{2}/L\simeq 14$ (Figure 13). We find good qualitative agreement between the two highest resolution simulations. However, there are quantitative differences. Particularly noticeable is the appearance of smaller scale fluid features as the resolution is increased. This is expected since we do not include an explicit viscosity in our simulations. Instead, the dissipation scale is related to the grid resolution, so that higher resolution simulations have a larger effective Reynolds number.

Zoom In Zoom Out Reset image size

We quantify this observation by computing the power spectrum of f_h on the original interface at z = 0 at ${{Agt}}^{2}/L\simeq 14$. We define the 2D spectrum of f_h as

$$\begin{eqnarray}&&{\hat{f}}_{h}({\boldsymbol{k}})=\iint {f}_{h}({\boldsymbol{x}}){| }_{z=0}\exp \left(2\pi i{\boldsymbol{k}}\cdot \displaystyle \frac{{\boldsymbol{x}}}{L}\right){d}^{2}{\boldsymbol{x}}\end{eqnarray} \tag{ 51 }$$

and the 1D spectrum as

$$\begin{eqnarray}&&{\hat{f}}_{h}(k)=\iint \delta ({\boldsymbol{k}}-k){\hat{f}}_{h}({\boldsymbol{k}}){d}^{2}{\boldsymbol{k}}.\end{eqnarray} \tag{ 52 }$$

Because of the periodic boundary conditions, the 2D spectrum is nontrivial only for ${\boldsymbol{k}}=({k}_{x},{k}_{y})$ with k_x, k_y integers. When computing the 1D spectrum, we convert the integral into a weighted summation following Eswaran & Pope (1988), i.e., we set

$$\begin{eqnarray}&&{\hat{f}}_{h}(k)=\displaystyle \frac{2\pi k}{{N}_{k}}\displaystyle \sum _{\parallel {\boldsymbol{k}}\parallel =k}{\hat{f}}_{h}({\boldsymbol{k}}),\end{eqnarray} \tag{ 53 }$$

where N_k is the number of modes with $\kappa -1/2\lt \parallel {\boldsymbol{k}}\parallel \lt \kappa \,+1/2$. Our results are shown in Figure 14.

Zoom In Zoom Out Reset image size

We find that the spectrum appears converged for k ∼ 2, as could have been anticipated by looking at the dipolar component of the oscillations in Figure 13. However, we find that fluctuations with scales k ∼ 10 are progressively suppressed as the resolution increases. The lack of convergence at these intermediate scales can be attributed to nonlinear mode coupling becoming stronger as the effective Reynolds number of the simulations increases. At smaller scales, the flow becomes self-similar and develops an inertial range. We find that the power spectrum follows the scaling expected from Kolmogorov's theory of turbulence for a passively advected scalar field, i.e., $| {\hat{f}}_{h}| \propto {k}^{-5/3}$ (Pope 2000). We also find that, as the resolution increases, a progressively larger part of the inertial range is uncovered. At scales of less than ∼5 grid points, the power spectrum drops due to the direct effect of numerical viscosity.

The Kelvin–Helmholtz (KH) instability is a two-dimensional instability arising from a velocity shear within a fluid. The instability results in vorticity, which can cascade into turbulence. As such, it is of particular interest in the study of CCSNe, where cascading turbulence can contribute to shock revival.

We consider a fluid over the rectangular domain within $x\in [0,4]$ and $z\in [0,2]$, with a setup similar to that found in Lecoanet et al. (2016) and with a density jump characterized by the initial conditions below. We conduct tests at low, medium, high resolutions, defined as 512 × 256, 1024 × 512, and 2048 × 1024, respectively.

$$\begin{eqnarray}&&\rho =1+\displaystyle \frac{\delta \rho }{2{\rho }_{0}}\left[\tanh \left(\displaystyle \frac{z-{z}_{1}}{a}\right)-\tanh \left(\displaystyle \frac{z-{z}_{2}}{a}\right)\right]\end{eqnarray} \tag{ 54a }$$

$$\begin{eqnarray}&&v={u}_{0}\times \left[\tanh \left(\displaystyle \frac{z-{z}_{1}}{a}\right)-\tanh \left(\displaystyle \frac{z-{z}_{2}}{a}\right)\right]\end{eqnarray} \tag{ 54b }$$

$$\begin{eqnarray}&&w=\mathrm{Asin}(\pi x)\times \left[\exp \left(-\displaystyle \frac{{\left(z-{z}_{1}\right)}^{2}}{{\sigma }^{2}}\right)+\exp \left(-\displaystyle \frac{{\left(z-{z}_{2}\right)}^{2}}{{\sigma }^{2}}\right)\right],\end{eqnarray} \tag{ 54c }$$

with a = 0.05, σ = 0.2, z₁ = 0.5, and z₂ = 1.5. We define v as the velocity in the x-direction and w as the velocity in the z-direction.

The initial amplitude, A, of the perturbation in the vertical (z-direction) velocity (w) is set to 0.01, and the initial density perturbation, δρ/ρ₀, is set to 1.0. The initial pressure P₀ is 10, the initial lateral velocity is constant, u₀, is set to 1.0, and the initial density ρ is set to 1.0. We assume an adiabatic gas with adiabatic index γ = 5/3. The resulting Mach number (M) is ∼0.25, consistent with quasi-incompressible flow. As is often done in KH instability tests, we ignore gravity.

In Figure 15, we plot as a function of time the growth rate of the maximum of w, the velocity in the z-direction. The linear growth phase of the KH instability lasts until t ∼ 2.6 (in our dimensionless units) for all models, which evolve virtually identically, independent of resolution, until the nonlinear regime.

Zoom In Zoom Out Reset image size

In Figure 16, we plot the 2D density evolution at three times (t = 0, 3, and 6 in dimensionless units) over our grid. Initial density perturbations evolve to form characteristic vortices at later times.

Zoom In Zoom Out Reset image size

The two-dimensional implosion test of Liska & Wendroff (2003b) consists of Sod-like initial conditions, but rotated by 45° in a two-dimensional reflecting box. Specifically, the domain extends from 0 to 0.3 along both the x- and y-directions, with (ρ, P) set to (1, 1) for x + y > 0.15 and (0.125, 0.14) otherwise. The gas is initially at rest and obeys an ideal gas equation of state with γ = 1.4. The solution is evolved until t = 2.5 on a uniform 400 × 400 mesh. This setup results in a shock, contact, and rarefaction moving from the initial discontinuity that reflect off the boundaries and interact. As highlighted by Stone et al. (2008), the correct solution to this problem includes a low-density jet that travels along the symmetry axis y = x. This problem is exquisitely sensitive to the preservation of this reflection symmetry, with the jet failing to form or wandering off y = x if symmetry is violated even at the level of round off in the discretized update. The directionally unsplit update in fornax is able to maintain this symmetry, as shown in Figure 17. Additionally, the distance traveled by the jet can be a useful measure of the numerical diffusion of contacts. This is illustrated in Figure 17 in the comparison between the results using linear and parabolic reconstructions. In fact, at much higher resolutions, the jet propagates into and interacts with the top-right corner of the domain, as shown in Figure 18. Our results are comparable, though perhaps characteristic of slightly higher numerical diffusion, to the results presented in Stone et al. (2008).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we investigate the nearly pressureless and homologous collapse of a uniformly dense cloud of dust under its own self-gravity. This problem has a solution described in Colgate & White (1966) and was used by Mönchmeyer & Müller (1989) to argue for the reconstruction of volume-averaged grid variables over the corresponding volume centroids, rather than the zone-centered coordinates. At a time t, the outer radius of the cloud, ${r}_{\mathrm{cl}}(t)$, is given by the solution to the equation

$$\begin{eqnarray}&&\left(\displaystyle \frac{8\pi G}{3}{\rho }_{0}\right)t={\left[\displaystyle \frac{{r}_{\mathrm{cl}}}{{r}_{0}}\left(1-\displaystyle \frac{{r}_{\mathrm{cl}}}{{r}_{0}}\right)\right]}^{1/2}+{\sin }^{-1}{\left(1-\displaystyle \frac{{r}_{\mathrm{cl}}}{{r}_{0}}\right)}^{1/2},\end{eqnarray} \tag{ 55 }$$

where ρ₀ and r₀ are the cloud's initial density and radius, respectively. Once r_cl is obtained, the cloud density ${\rho }_{\mathrm{cl}}(t)$ is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{cl}}={\rho }_{0}{\left(\displaystyle \frac{{r}_{\mathrm{cl}}}{{r}_{0}}\right)}^{-3},\end{eqnarray} \tag{ 56 }$$

and the radial velocity profile inside the cloud, ${v}_{\mathrm{cl}}(r,t)$, is given by

$$\begin{eqnarray}&&{v}_{\mathrm{cl}}=-r{\left[\displaystyle \frac{8\pi G}{3}{\rho }_{0}\left(\displaystyle \frac{{r}_{0}}{{r}_{\mathrm{cl}}}-1\right)\right]}^{1/2}.\end{eqnarray} \tag{ 57 }$$

We use ${\rho }_{0}={10}^{9}\ {\rm{g}}\ {\mathrm{cm}}^{-3}$, ${r}_{0}=6500\,\mathrm{km}$, and an adiabatic equation of state with γ = 5/3, and evolve to ${t}_{\mathrm{final}}=0.065\ {\rm{s}}$. The pressure is set to a negligibly small, constant value, P₀, such that ${P}_{0}\ll (4\pi G/\gamma ){\rho }_{0}^{2}{r}_{0}^{2}$, hence at t_final, the outer edge of the cloud is collapsing at a Mach number of ∼80. Finally, we use a three-dimensional spherical dendritic grid with radial extent out to ${r}_{\max }=7000\,\mathrm{km}$ and a resolution of 200 × 64 ×128 zones. The radial grid spacing is constant in the interior with ${\rm{\Delta }}r\approx 0.5\,\mathrm{km}$ out to approximately 15 km, then smoothly transitions to logarithmic spacing out to r_max.

Figure 19 shows the density, ρ, and radial velocity, v, at time t_final. We scale these by the cloud density and radius, ρ_cl and v_cl, respectively, as obtained from Equations (56) and (57), which in turn depend on the numerical solution for r_cl obtained from Equation (55) at time t_final. As evident from the figure, the cloud density remains constant throughout the cloud, and the velocity remains linear. There is some smearing at the outer edge of the cloud due to finite pressure gradients, but otherwise, the results are comparable to those obtained by other authors (Mönchmeyer & Müller 1989; Mignone 2014). This test demonstrates the code's ability to obtain accurate fluxes from reconstruction in the volume coordinate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We perform several spherical 1D core-collapse hydrodynamical simulations without transport or GR corrections and with Newtonian gravity to study energy conservation in fornax. For this test, we use the SFHo equation of state (Steiner et al. 2013). The total mechanical energy is defined as the sum of the kinetic, internal, and gravitational potential energies on the grid, extending out to 20,000 km. The latter is calculated as

$$\begin{eqnarray}&&{E}_{\mathrm{grav}}=-\displaystyle \frac{1}{8\pi G}\int | {\rm{\nabla }}{\rm{\Phi }}{| }^{2}{d}^{3}x,\end{eqnarray} \tag{ 58 }$$

where Φ is the gravitational potential. We take the jump in total energy around bounce as an apt measure, since gravitation is not implemented in automatically conservative form (the other terms are), and bounce is the most problematic phase of core-collapse simulations vis á vis total energy conservation. Moreover, we continue the calculations to ∼100 ms after bounce for three resolutions: Lo (304 radial cells, green), Default (608 radial cells, blue), and Hi (1216 radial cells, red). We find that energy is conserved from just before bounce to just after bounce to approximately 10⁵⁰ erg, $2\times {10}^{49}\,\mathrm{erg}$, and $\sim 3\times {10}^{48}\ \mathrm{erg}$, respectively, for the three resolutions (Figure 20).

Earlier generations of codes, including AGILE-BOLTZTRAN (Liebendörfer et al. 2004) and BETHE-hydro (Murphy & Burrows 2008b), saw total energy shifts from just before bounce to over 100 mas after bounce of on the order of 10⁵¹ erg. For comparison, CHIMERA (Bruenn et al. 2009, 2016) conserves energy to 0.5 Bethe ($1\ {\rm{B}}={10}^{51}\,\mathrm{erg}$) over ∼1 s post-bounce. Our default resolution from before bounce to ∼1 s after bounce conserves total energy to 0.05 B, an order of magnitude better. Importantly, we find that the total energy is conserved to $\sim {10}^{47}\,\mathrm{erg}$ from 10 ms post-bounce to 100 ms post-bounce. By comparison, using CoCoNUT for a Newtonian core-collapse simulation and subtracting out neutrino losses, Müller et al. (2010) report energy conservation of $\sim 2\times {10}^{49}\,\mathrm{erg}$ from just before bounce to just after bounce, and subsequent non-conservation for the following tens of milliseconds at the ∼10⁴⁸ level. In Section 9.1, we perform an analogous core-collapse test using multidimensional, multispecies transport.

As a first test of our full radiation-hydrodynamics code, we test the conservation of total energy by performing a 2D core-collapse calculation on a grid with a resolution of 608 radial cells by 128 polar-angle cells. We use Newtonian monopolar gravity, the SFHo equation of state, 12 energy groups per species, and our dendritic grid decomposed over 1024 processors. We track the gravitational potential energy (defined as in Equation (58)), internal plus kinetic energy, and the total laboratory-frame radiation energy at each time step, accounting for the flux of each energy component through the outer radial boundary.

The results of this test are shown in Figure 21. As previously demonstrated in Section 8.9, we experience a glitch in the total energy at bounce, since gravity is treated in non-conservation form. With radiation transport, there are other additional sources of energy non-conservation, e.g., due to the finite tolerance of the implicit solver or due to the truncation of the comoving-to-lab frame transformation in powers of v/c. However, as Figure 21 demonstrates, the total energy remains well controlled to within a few $\times \sim {10}^{50}\,\mathrm{erg}$ (or to within 0.5% relative to E_grav) throughout the 500 ms of the simulation run.¹² Energy conservation within such low levels is especially important considering the high-velocity energy fluxes crossing the internal refinement boundaries of our dendritic grid in both the radial and angular directions and across a distribution of many processors.

Zoom In Zoom Out Reset image size

To test the advection of radiation energy in frequency space caused by strong spatial variations in the gas velocity, we perform the Doppler-shift test of Vaytet et al. (2011). We use a one-dimensional domain with 50 equally spaced zones between x = 0 cm and x = 10 cm with constant background density of ρ = 1 g cm⁻³ and zero opacity. The velocity varies smoothly between $v=0\ \mathrm{cm}\,{{\rm{s}}}^{-1}$ at the left edge of the domain and v = A in the center, where $A=5\times {10}^{7}\ \mathrm{cm}\ {{\rm{s}}}^{-1}$. At the left boundary, we introduce a radiation field having a Planck spectrum at a temperature T = 1000 K, which we resolve using 20 equally spaced frequency groups from ν = 0 Hz to ν = 2 × 10¹⁴ Hz, with two additional groups to account for radiation in the frequency range $2\times {10}^{14}\ \mathrm{Hz}\to \infty$. We evolve the system for 10 light-crossing times in order to reach a steady state.

In Figure 22, we plot the difference between the Doppler-shifted spectral energy distribution at the center of the domain, ${{ \mathcal E }}_{\nu }(v=A)$, and the unshifted distribution at the left boundary, ${{ \mathcal E }}_{\nu }(v=0)$. We also plot the analytic solution, obtained from the difference between the relativistic Doppler shift of the incident spectrum into the frame comoving at velocity v = A and the unshifted spectrum. To do this, we first multiply the laboratory-frame frequencies by the relativistic Doppler factor,

$$\begin{eqnarray}&&\sqrt{\displaystyle \frac{1+\beta }{1-\beta }},\end{eqnarray} \tag{ 59 }$$

where $\beta =A/c$, to obtain the corresponding comoving-frame frequencies. Then, we set ${{ \mathcal E }}_{\nu }=(4\pi /c){B}_{\nu }$, with frequencies evaluated in the comoving frame, and finally divide by the relativistic Doppler factor to transform the resulting spectrum back into the laboratory frame.

Zoom In Zoom Out Reset image size

The computed solution agrees well with the analytic solution at almost all frequencies. The relative error is slightly larger near $\nu \sim 5\times {10}^{13}\ \mathrm{Hz}$, the peak of the Planck spectrum at T = 1000 K, where TVD slope limiting in our frequency advection step results in a local first-order error. With 20 groups, the relative error is less than ∼6.5% away from the spectral peak, where it jumps up to ∼23%. Running the same test with 40 groups reduces the maximum relative error to 0.85% away from the spectral peak, where it jumps to 5.6%.

9.2.1. Continuity of Laboratory-frame Flux at Transparent Shocks

To demonstrate further that fornax accounts properly for frame effects, we show a comparison of the radial profile of the total radiation energy density and flux calculated in the steady state for a constant point luminosity source at the spherical grid center. We set the material opacity to zero and impose a core-collapse-like radial velocity field given by

$$\begin{eqnarray}v(r)=\left\{\begin{array}{ll}-0.1\,{v}_{\mathrm{shock}}\displaystyle \frac{r}{{r}_{\mathrm{shock}}}, & r\lt {r}_{\mathrm{shock}}\\ -{v}_{\mathrm{shock}}\displaystyle \frac{\sqrt{{r}_{\max }/r}-1}{\sqrt{{r}_{\max }/{r}_{\mathrm{shock}}}-1}, & r\geqslant {r}_{\mathrm{shock}}\end{array}\right.,\end{eqnarray} \tag{ 60 }$$

where r_shock = 100 km is the location of a strong standing shock of maximum velocity v_shock = 0.1c and r_max = 1000 km is the maximum radius.

Since fornax calculates radiation quantities in the comoving frame, the corresponding energy density and flux profiles will have discontinuities at shocks. However, since the radiation is uncoupled from the matter, the radiation quantities should possess a smooth profile that goes as 1/r² in the laboratory frame. Thus, performing a simple Lorentz transformation back to the laboratory frame should eliminate the discontinuity. Figure 23 depicts the radial profiles of the radiation energy density and flux, each scaled to their values in the first radial zone, with the left and right panels depicting the profile in the comoving and laboratory frames, respectively. The right-hand side panel shows that the Lorentz transformation has indeed removed the discontinuity. This simple test is one way of demonstrating that fornax handles the velocity-dependent terms in radiation Equations 3(a) and (b) properly.

Zoom In Zoom Out Reset image size

To test the effect of a strong velocity gradient on a frequency-dependent opacity function, we perform a modified version of a test by Vaytet et al. (2011). We use a one-dimensional domain with $x\in [0,1]\ \mathrm{cm}$ resolved over 100 zones and set a fixed velocity profile $v={ \mathcal D }x$, where ${ \mathcal D }$ is either 0 s⁻¹ or 10⁷ s⁻¹. We use a density profile of $\rho =1/({ \mathcal C }x)$ with ${ \mathcal C }=0.2\ {\mathrm{cm}}^{2}\ {{\rm{g}}}^{-1}$ and set the gas temperature to T₀ = 3 K everywhere. We use a frequency-dependent opacity with ${\kappa }_{\nu }\,=100\ {\mathrm{cm}}^{2}\ {{\rm{g}}}^{-1}$ for $\nu \lt 2\times {10}^{13}\ \mathrm{Hz}$, transitioning smoothly to ${\kappa }_{\nu }=1\ {\mathrm{cm}}^{2}\ {{\rm{g}}}^{-1}$ over a region of width ${\rm{\Delta }}\nu =4.5\times {10}^{9}\ \mathrm{Hz}$ as shown in Figure 24. At the left boundary, we inject a Gaussian radiation spectrum that is normalized to have total energy ${a}_{{\rm{R}}}{T}_{1}^{4}$ with ${T}_{1}=1000\ {\rm{K}}$, is peaked at $\nu =2\times {10}^{13}\ \mathrm{Hz}$, has FWHM equal to two-thirds the width of the transition region, and is resolved over 20 equally spaced frequency groups (also shown in Figure 24). For each case, we evolve only the radiation system for one light-crossing time, keeping the hydrodynamics frozen.

Zoom In Zoom Out Reset image size

Figure 25 shows selected group temperatures, ${T}_{g}\equiv {({{ \mathcal E }}_{g}/{a}_{{\rm{R}}})}^{1/4}$, for the cases with ${ \mathcal D }=0\ {{\rm{s}}}^{-1}$ (top) and ${ \mathcal D }={10}^{7}\ {{\rm{s}}}^{-1}$ (bottom). In the zero-velocity case, the group energies are unshifted, and only those groups with low-enough optical depth are able to stream freely across the domain. By contrast, in the case with a strong velocity gradient, the comoving-frame frequencies are Doppler-shifted to the right where the opacities are lower. Thus, the energy in groups 9–12 in the transition region of the opacity function is able to stream more freely once optically thin conditions are reached.

Zoom In Zoom Out Reset image size

To study the behavior of fornax in the strong diffusion regime, where radiation is transported both by diffusion through the gas and advection along with it, we perform a variation of the pulse advection test of Krumholz et al. (2007; see also Zhang et al. 2013), modified for our explicit transport solver. We initialize the gas temperature and density profiles as

$$\begin{eqnarray}&&T={T}_{0}+({T}_{1}-{T}_{0})\exp \left(-\displaystyle \frac{{x}^{2}}{2{w}^{2}}\right),\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&\rho ={\rho }_{0}\displaystyle \frac{{T}_{0}}{{T}_{1}}+\displaystyle \frac{{a}_{{\rm{R}}}\mu }{3{k}_{{\rm{B}}}}\left(\displaystyle \frac{{T}_{0}^{4}}{T}-{T}^{3}\right),\end{eqnarray} \tag{ 62 }$$

where T₀ and ρ₀ are the background temperature and density, respectively; T₁ = 2T₀ is the peak of a Gaussian pulse of thermal energy of width w; and $\mu =({m}_{e}+{m}_{p})/2$ is the mean particle mass for ionized hydrogen. Following Zhang et al. (2013), we use a temperature- and frequency-dependent absorption opacity of the form

$$\begin{eqnarray}&&{\kappa }_{\nu }(T)={\kappa }_{0}{\left(\displaystyle \frac{T}{{T}_{1}}\right)}^{-1/2}{\left(\displaystyle \frac{\nu }{{\nu }_{1}}\right)}^{-3}\left[1-\exp \left(-\displaystyle \frac{h\nu }{{k}_{{\rm{B}}}T}\right)\right],\end{eqnarray} \tag{ 63 }$$

where κ₀ is some normalization constant, and we take ${\nu }_{1}=2.821{k}_{{\rm{B}}}{T}_{1}/h$ to approximate the spectral peak of a Planck distribution at temperature T₁. The problem dimensions are set by the choice of an arbitrary length scale, w, and by the values of the non-dimensional parameters β ≡ v/c, ${ \mathcal M }\equiv v/{a}_{0}$, ${ \mathcal P }\,\equiv {a}_{{\rm{R}}}{T}_{0}^{4}/({\rho }_{0}{a}_{0}^{2})$, and $\tau \equiv {\kappa }_{{\rm{P}}}({T}_{1})w$, where ${\kappa }_{{\rm{P}}}=3.457\,{\kappa }_{0}{(T/{T}_{1})}^{-3.5}$ is the Planck mean opacity such that τ represents the optical depth of the pulse. For our version of this problem, we choose $w=1\ \mathrm{cm}$, β = 0.01, ${ \mathcal M }=0.1$, ${ \mathcal P }=0.1$, and ${\kappa }_{0}=0.2892\,\tau \ {\mathrm{cm}}^{-1}$, where τ is chosen to control which diffusion regime characterizes the flow.

The product βτ is equal to the ratio of the radiation diffusion and advection timescales. For βτ ≲ 1, radiation diffusion dominates, and the system is in the so-called "static diffusion" regime, but for βτ ≫ 1, advection dominates and the system is in the "dynamic diffusion" regime. In either case, although the velocity-dependent radiation work and advection terms will be very different depending on the frame in which the equations are solved, the results should be frame-independent. Therefore, for both the static and dynamic diffusion regimes, we compare runs in which the radiation pulse is initially at rest to runs in which the radiation pulse is advected over twice its initial width and shifted back to lie on top of the unadvected results. In either case, with τ ≫ 1, the radiation is in thermal equilibrium with the gas, and the total pressure is initially constant.

We use a computational grid of length L = 20w resolved over 512 zones, with eight radiation groups logarithmically spaced over the frequency range $[4\times {10}^{18},4\times {10}^{22}]\ \mathrm{Hz}$. In Figure 26, we show the resulting temperature, velocity, and density profiles for both the advected and unadvected runs with τ = 100 such that βτ = 1, putting the system in the static diffusion regime. Over the course of the run, the radiation has diffused through the gas, decreasing the temperature and thermal pressure support in the pulse, hence the density increases correspondingly. Since the advected and unadvected runs are visually indistinguishable, we plot the relative errors of the density and temperature in Figure 27, which are bounded by 2.3 × 10⁻⁵ and 1.4 × 10⁻⁵, respectively, in absolute value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Similarly, in Figure 28, we plot the results for the case with τ = 10⁴ such that βτ = 100, which puts the system in the dynamic diffusion regime. This time, there is hardly any diffusion as the pulse is advected, since the radiation is effectively trapped within the gas. In Figure 29, we plot the relative errors of the density and temperature, which are bounded by 1.4 × 10⁻³ and 1.5 × 10⁻⁴, respectively, in absolute value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Importantly, these results demonstrate the ability of fornax to reproduce the correct behavior to a high degree of accuracy in both the static and dynamic diffusion regimes and for both advected and unadvected flows. This represents a very sensitive test of the code's ability to model radiative trapping in the CCSN. In the CCSN context, trapping of ${\nu }_{e}{\rm{s}}$ by infalling matter and their subsequent compression produces a spectrum of degenerate ν_e whose energies can reach the ∼100–300 MeV range while preserving a relatively high lepton fraction, which is a critically important aspect that allows the core to reach nuclear densities at bounce.

To test the fornax code's ability to couple the multigroup radiation subsystem to the hydrodynamics in a regime where the fluid and radiation are out of equilibrium, we present here the results of the classical gray non-equilibrium radiation shock described by Zel'Dovich & Raizer (1969). Following that work, we adopt the Eddington approximation and set the radiation pressure tensor to $({{ \mathcal P }}_{\epsilon }{)}_{j}^{i}=({{ \mathcal E }}_{\epsilon }/3){\delta }_{j}^{i}$, where ${\delta }_{j}^{i}$ is the Kronecker symbol. The shock structure has a semi-analytic solution described by Lowrie & Edwards (2008). Using their parameters with an upstream Mach number of ${ \mathcal M }=3$, we obtain a subcritical radiation shock whose temperature jumps discontinuously at the shock interface. We set the constant, gray absorption opacity to $\kappa =577\ {\mathrm{cm}}^{-1}$, the mean particle mass to $\mu ={m}_{{\rm{H}}}$, and use an adiabatic EOS with γ = 5/3. We use a one-dimensional Cartesian domain with $x\,\in [-0.0132,0.00255]$ resolved over N_x = 512 zones. The problem is set in the rest frame of the shock, which we initialize at x = 0. In the upstream state (x < 0), we set the gas temperature ${T}_{0}=2.18\times {10}^{6}\ {\rm{K}}$, density ${\rho }_{0}=5.69\ {\rm{g}}\ {\mathrm{cm}}^{-3}$, and velocity ${v}_{0}=5.19\times {10}^{7}\ \mathrm{cm}\ {{\rm{s}}}^{-1}$, and using the Rankine–Hugoniot jump conditions to determine the downstream state (x < 0), we set ${T}_{1}=7.98\times {10}^{6}\ {\rm{K}}$, density ${\rho }_{1}=17.1\ {\rm{g}}\ {\mathrm{cm}}^{-3}$, and velocity ${v}_{1}=1.73\times {10}^{7}\ \mathrm{cm}\ {{\rm{s}}}^{-1}$. Similar to Vaytet et al. (2011), we use 8 radiation groups logarithmically spaced over the frequency range $\nu \in [{10}^{15}\ \mathrm{Hz},{10}^{19}\ \mathrm{Hz}]$, initialize the group radiation energy densities using a Planck spectrum at the local gas temperature with zero radiation flux, and evolve for 3 shock-crossing times to ${t}_{\mathrm{final}}=9.08\times {10}^{-10}\ {\rm{s}}$. Finally, since the structure of the steady-state shock solution is independent of the radiation propagation speed, $\hat{c}$, and since we must solve the semi-explicit radiation subsystem on a hydrodynamic timescale, we adopt a reduced speed of light $\hat{c}\,=10({v}_{0}+{a}_{0})$, where ${a}_{0}=1.73\times {10}^{7}\ \mathrm{cm}\ {{\rm{s}}}^{-1}$ is the upstream sound speed.

Figure 30 shows the gas temperature structure at time t_final with the semi-analytic solution from Lowrie & Edwards (2008) overplotted and an inset showing the detail of the Zel'Dovich spike in the upstream temperature near the shock interface. There is very good agreement between the two solutions, including in the non-equilibrium spike region and in the radiatively heated shock precursor in the downstream state. The relative error between the solutions is bounded by 0.8% everywhere; the agreement in the remaining variables is similarly good. Since this test simultaneously exercises the multigroup, velocity-dependent, Doppler-shift, and matter–radiation coupling features of the fornax code, it also simultaneously demonstrates the code's ability to compute the most problematic aspects of full radiation–hydrodynamics accurately and consistently.

Zoom In Zoom Out Reset image size