CASTRO: A NEW COMPRESSIBLE ASTROPHYSICAL SOLVER. III. MULTIGROUP RADIATION HYDRODYNAMICS

W. Zhang; L. Howell; A. Almgren; A. Burrows; J. Dolence; J. Bell

doi:10.1088/0067-0049/204/1/7

1. INTRODUCTION

Numerical simulations are a useful tool for many radiation hydrodynamic problems in astrophysics. However, numerical modeling of radiation hydrodynamic phenomena is very challenging for a number of reasons. To achieve a highly accurate description of the system usually requires a multidimensional treatment with high spatial resolution. It is also desirable to handle the radiation frequency variable properly because radiation transport is a frequency-dependent process. Moreover, the numerical algorithms for radiation hydrodynamics must be stable and efficient.

Core-collapse supernovae are such complex phenomena (Bethe 1990; Kotake et al. 2006; Janka et al. 2007; Janka 2012) that neutrino radiation hydrodynamics simulations have been indispensable in helping to understand the explosion mechanism. Various algorithms have been developed to solve neutrino radiation hydrodynamics equations. Because of the complexities involved in the problem and the limited computer resources available, various tradeoffs have to be made. Many of these algorithms are one-dimensional (1D; e.g., Bowers & Wilson 1982; Bruenn 1985; Mezzacappa & Bruenn 1993; Burrows et al. 2000; Liebendörfer et al. 2004; O'Connor & Ott 2010), or have used the so-called "ray-by-ray" approach, which does not truly perform multidimensional transport (e.g., Burrows et al. 1995; Rampp & Janka 2002; Buras et al. 2006; Takiwaki et al. 2012). The two-dimensional (2D) code of Miller et al. (1993) is a gray flux-limited diffusion (FLD) code that does not include order O(v/c) terms, where v is the characteristic velocity of the system. The 2D multigroup multiangle code of Livne et al. (2004) does not include all order O(v/c) terms, and in particular it lacks direct coupling among radiation groups. The 2D algorithm of Hubeny & Burrows (2007) is based on a mixed-frame formulation of a multigroup two-moment system that is correct to order O(v/c), but it remains to be implemented in a working code. Swesty & Myra (2009) have developed a fully coupled 2D multigroup FLD code based on a comoving-frame formulation that is correct to order O(v/c). Besides the grid-based methods, smoothed particle hydrodynamics (SPH) methods have also been applied to multidimensional neutrino radiation hydrodynamics simulations of core-collapse supernovae (e.g., Herant et al. 1994; Fryer & Warren 2002; Fryer et al. 2006), but the SPH scheme has not yet been extended to a multigroup. Recently, Abdikamalov et al. (2012) developed a Monte Carlo method for neutrino transport, but hydrodynamics and radiation transport are not yet coupled in the scheme.

In this paper, we present a numerical algorithm for solving multidimensional multigroup photon and neutrino radiation hydrodynamics equations in the comoving frames that are correct to order O(v/c). Radiation quantities that are of interest are the radiation energy density, radiation flux, and radiation pressure tensor. It is more accurate to evolve both the radiation energy density and the radiation flux especially when the radiation is highly anisotropic and optically thin. However, the two-moment approach is very expensive in terms of computer time and memory for multidimensional multigroup radiation hydrodynamics simulations. Thus, we have adopted the FLD approach (Alme & Wilson 1973) for computational efficiency. In the FLD approach, only the radiation energy density needs to be evolved over time, and the radiation flux is derived through a diffusion approximation. Furthermore, we use an analytic closure for the radiation pressure tensor. We have adopted a multigroup method, in which the radiation frequency is discretized into multiple groups. The equations that we solve are similar to those in Swesty & Myra (2009). However, we have neglected inelastic scattering of neutrinos on matter, which is currently treated as elastic scattering. Since the total cross sections for the latter at the neutrino energies near and outside the neutrinospheres are ∼20–100 times smaller than the cross sections for scattering off nucleons and nuclei, not including inelastic effects at this stage has a minor effect on the qualitative behavior of collapse, bounce, and shock dynamics. Nevertheless, the approach of Thompson et al. (2003) to handle inelasticity and energy redistribution can be incorporated into our scheme. Another difference between our scheme and the scheme of Swesty & Myra (2009) is that we have developed a Godunov scheme in which radiation is fully coupled into a characteristic-based Riemann solver for the hyperbolic part of the system.

Many multigroup neutrino transport codes use a fixed Eulerian grid. However, it is difficult to achieve sufficient resolution in multidimensional simulations with this simple approach. To efficiently achieve high resolution, we use a Eulerian grid with block-structured adaptive mesh refinement (AMR). The AMR algorithm we use has a nested hierarchy of logically rectangular variable-sized grids and refines simultaneously in both space and time. AMR techniques have been successfully used in multidimensional multigroup photon radiation transport (Gittings et al. 2008; van der Holst et al. 2011). However, to the best of our knowledge, this is the first time that AMR has been employed in a multidimensional multigroup neutrino solver. Besides AMR, the arbitrary Lagrange–Eulerian (ALE) method is an alternative approach for mesh refinement. We note that ALE has been successfully employed in multidimensional multigroup neutrino transport (Livne et al. 2004; Burrows et al. 2007; Ott et al. 2008) and photon transport (Marinak et al. 2001).

We have implemented our algorithm in the compressible astrophysics code, CASTRO. This is the third paper in a series on the CASTRO code and its numerical algorithms. In our previous papers, we describe our treatment of hydrodynamics, including gravity and nuclear reactions (Almgren et al. 2010, henceforth Paper I), and our algorithm for flux-limited gray radiation hydrodynamics based on a mixed-frame formulation (Zhang et al. 2011, henceforth Paper II). Here, we describe our algorithm for flux-limited multigroup radiation hydrodynamics based on a comoving-frame formulation.

In Section 2, we present the multigroup radiation hydrodynamics equations that CASTRO solves. We show that the system can be divided into a hyperbolic subsystem (Section 2.3), a subsystem of advection in frequency space (Section 2.4), and a parabolic subsystem for radiation diffusion (Section 2.5). Analytic results for the mathematical characteristics of the hyperbolic subsystem are also presented in Section 2.3. The hyperbolic and the frequency space advection steps are solved by an explicit solver, whereas the parabolic subsystems are solved by an iterative implicit solver. We describe the explicit solvers in Section 3.1, followed by a description of the implicit solver in Section 3.2. We also discuss acceleration schemes that can greatly improve the convergence rate of the iterative implicit solver (Section 3.2.4). In Section 4, we present results from a series of test problems. A scaling test is shown in Section 5 to demonstrate the scaling behavior of the multigroup radiation solver. Finally, the results of the paper are summarized in Section 6.

2. MULTIGROUP RADIATION HYDRODYNAMICS

In this section, we present the multigroup radiation hydrodynamics equations that CASTRO solves. Our formulation is based on a comoving-frame approach, unlike our gray radiation hydrodynamics solver presented in Paper II. In the comoving-frame approach, the radiation quantities, and the opacities are measured in the comoving frames. The set of equations are correct to order O(v/c). We then derive the multigroup radiation hydrodynamics equations using frequency space discretization. We also analyze the mathematical characteristics of the system. Here we present the results in terms of neutrino radiation and indicate modifications needed for photon radiation.

2.1. Frequency-dependent Radiation Hydrodynamics Equations

The system of neutrino radiation hydrodynamics that is correct to order O(v/c) can be described by

$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \boldsymbol{u}) = { } & 0, \end{eqnarray} \tag{ 1 }$

$\begin{eqnarray} \frac{\partial (\rho \boldsymbol{u})}{\partial t} + \nabla \cdot (\rho \boldsymbol{u} \boldsymbol{u}) + \nabla p = { } & \int _0^{\infty } \frac{\chi }{c} \boldsymbol{F}_\nu d\nu + \boldsymbol{F}_G, \end{eqnarray} \tag{ 2 }$

$\begin{eqnarray} \frac{\partial (\rho E)}{\partial t} + \nabla \cdot (\rho E \boldsymbol{u} + p \boldsymbol{u}) = { } & \boldsymbol{u} \cdot \int _0^{\infty } \frac{\chi }{c} \boldsymbol{F}_\nu d\nu + \int _0^{\infty } (c\kappa E_{\nu } - 4\pi \eta) d\nu + \boldsymbol{u} \cdot \boldsymbol{F}_G, \end{eqnarray} \tag{ 3 }$

$\begin{eqnarray} \frac{\partial (\rho Y_e)}{\partial t} + \nabla \cdot (\rho Y_e \boldsymbol{u}) = { } & \int _0^{\infty } \xi (c\kappa E_{\nu } - 4\pi \eta) d\nu , \end{eqnarray} \tag{ 4 }$

$\begin{eqnarray} \frac{\partial E_{\nu }}{\partial t} + \nabla \cdot (E_{\nu } \boldsymbol{u}) = { } & -\nabla \cdot \boldsymbol{F}_\nu - (c\kappa E_{\nu } - 4\pi \eta) + \frac{\partial }{\partial \ln {\nu }} (\mathbf{\mathsf {P}}_{\nu } : \nabla \boldsymbol{u}) \end{eqnarray} \tag{ 5 }$

(Buchler 1983; Mihalas & Mihalas 1999; Castor 2004; Swesty & Myra 2009). Here, ρ, u, p, and E are the mass density, velocity, pressure, and total matter energy per unit mass (internal energy, e, plus kinetic energy, $\boldsymbol{u}^2/2$ ), respectively. Y_e is the electron fraction defined to be the net electron number (electrons minus positrons) per baryon. F_G is the external force on the fluid (e.g., gravitational force). E_ν, $\boldsymbol{F}_\nu$ , and $\mathbf{\mathsf {P}}_\nu$ are the monochromatic radiation energy density, radiation flux, and radiation pressure tensor at frequency ν, measured in the comoving fluid frame, respectively. The speed of light is denoted by c. The absorption coefficient is denoted by κ, and the total interaction coefficient including conservative scattering is given by χ. The emission coefficient is denoted by η. For photon radiation ξ is zero, whereas for neutrinos ξ is given by

$\begin{equation} \xi = s \frac{m_{B}}{h \nu }, \end{equation} \tag{ 6 }$

where m_B is the baryon mass, h is the Planck constant, and s = 1 for ν_e neutrinos, s = −1 for $\bar{\nu }_e$ neutrinos, and s = 0 for other flavor neutrinos. Also note that for neutrino radiation with multiple species, the integrals over frequency imply summation over the species. The colon product : in $\mathbf{\mathsf{P}}_{\nu } : \nabla \boldsymbol{u}$ indicates contraction over two indices as follows:

$\begin{equation} \mathbf{\mathsf {P}}_{\nu } : \nabla \boldsymbol{u} = \left(\mathbf {\mathsf {P}}_{\nu }\right)_{ij} \left(\nabla \boldsymbol{u}\right)_{ij}. \end{equation} \tag{ 7 }$

The system of the radiation hydrodynamics (Equations (1)–(5)) is closed by an equation of state (EOS) for matter, the diffusion approximation for radiation flux, and a closure relation for radiation that computes the radiation pressure tensor from radiation energy density and radiation flux. In the FLD approximation (Alme & Wilson 1973), the monochromatic radiation flux is written in the form of Fick's law of diffusion,

$\begin{equation} \boldsymbol{F}_{\nu } = - \frac{c \lambda }{\chi } \nabla E_{\nu }, \end{equation} \tag{ 8 }$

where λ is the flux limiter at frequency ν. We have implemented a variety of choices for the flux limiter (e.g., Minerbo 1978; Bruenn et al. 1978; Levermore & Pomraning 1981). For example, one form of the flux limiter in Levermore & Pomraning (1981) is

$\begin{equation} \lambda = \frac{2+R}{6+3R+R^2}, \end{equation} \tag{ 9 }$

where

$\begin{equation} R = \frac{|\nabla E_{\nu }|}{\chi E_{\nu }}. \end{equation} \tag{ 10 }$

In the optically thick limit, R → 0, λ → 1/3, and the flux approaches the classical Eddington approximation $\boldsymbol{F}_\nu \approx - (c/3\chi) \nabla E_\nu$ . In the optically thin limit, R → ∞, λ ≈ 1/R, and $\boldsymbol{F}_\nu \approx c E_\nu$ as expected for streaming radiation. For the radiation pressure tensor, we use an approximate closure (Minerbo 1978; Levermore 1984)

$\begin{equation} \mathbf {\mathsf {P}}_{\nu } = \frac{1}{2}(1-f) E_{\nu } \mathbf {\mathsf {I}} + \frac{1}{2}(3 f-1) E_{\nu } \hat{\boldsymbol{n}}\hat{\boldsymbol{n}}, \end{equation} \tag{ 11 }$

where

$\begin{equation} f = \lambda + \lambda ^2 R^2 \end{equation} \tag{ 12 }$

is the scalar Eddington factor, $\mathbf {\mathsf {I}}$ is the identity tensor of rank 2, and $\hat{\boldsymbol{n}} = \nabla E_{\nu } / |\nabla E_{\nu }|$ . There are other choices for the scalar Eddington factor. For example, Kershaw (1976) has suggested

$\begin{equation} f = \frac{1}{3} + \frac{2}{3} \lambda ^2 R^2. \end{equation} \tag{ 13 }$

Note that in the optically thick limit, f → 1/3, and in the optically thin limit, f → 1.

Using the FLD approximation and the approximate closure, the equations of the radiation hydrodynamics now become

$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \boldsymbol{u}) = { } & 0, \end{eqnarray} \tag{ 14 }$

$\begin{eqnarray} \frac{\partial (\rho \boldsymbol{u})}{\partial t} + \nabla \cdot (\rho \boldsymbol{u} \boldsymbol{u}) + \nabla p = { } & -\int _0^{\infty } \lambda \nabla E_{\nu } d\nu + \boldsymbol{F}_G, \end{eqnarray} \tag{ 15 }$

$\begin{eqnarray} \frac{\partial (\rho E)}{\partial t} + \nabla \cdot (\rho E \boldsymbol{u} + p \boldsymbol{u}) = { } & -\boldsymbol{u} \cdot \int _0^{\infty } \lambda \nabla E_{\nu } d\nu + \int _0^{\infty } (c\kappa E_{\nu } - 4 \pi \eta) d\nu +\boldsymbol{u} \cdot \boldsymbol{F}_G, \end{eqnarray} \tag{ 16 }$

$\begin{eqnarray} \frac{\partial (\rho Y_e)}{\partial t} + \nabla \cdot (\rho Y_e \boldsymbol{u}) = { } & \int _0^{\infty } \xi (c\kappa E_{\nu } - 4\pi \eta) d\nu , \end{eqnarray} \tag{ 17 }$

$\begin{eqnarray} \frac{\partial E_{\nu }}{\partial t} + \nabla \cdot (E_{\nu } \boldsymbol{u}) = { } & \nabla \cdot \left(\frac{c\lambda }{\chi } \nabla E_{\nu } \right) - (c\kappa E_{\nu } - 4\pi \eta) \nonumber \\ { } & + \frac{\partial }{\partial \ln {\nu }} \left(\frac{1-f}{2} E_\nu \nabla \cdot \boldsymbol{u} + \frac{3f-1}{2} E_\nu \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}\right). \end{eqnarray} \tag{ 18 }$

2.2. Multigroup Radiation Hydrodynamics Equations

To numerically solve Equations (14)–(18), we first discretize the system in frequency space by dividing the radiation into N energy groups. We define the radiation energy density in each group as

$\begin{equation} E_g = \int _{\nu _{g-1/2}}^{\nu _{g+1/2}} E_\nu d\nu , \end{equation} \tag{ 19 }$

where ν_{g − 1/2} and ν_{g + 1/2} are the frequency at the lower and upper boundaries, respectively, for the gth group. For clarity, $\int _{\nu _{g-1/2}}^{\nu _{g+1/2}} d\nu$ is denoted by ∫_gdν hereafter.

In the multigroup methods, various group mean coefficients are introduced to the system. For example, the Rosseland mean interaction coefficient is defined as

$\begin{equation} \chi _g^{-1} = \frac{\int _g \chi ^{-1} (\partial B_\nu /\partial T) d\nu }{\int _g (\partial B_\nu /\partial T) d\nu }, \end{equation} \tag{ 20 }$

where B_ν = η/κ. In this paper, we will simply assume that

$\begin{eqnarray} \kappa _g = { } & \kappa (\nu _g), \end{eqnarray} \tag{ 21 }$

$\begin{eqnarray} \chi _g = { } & \chi (\nu _g), \end{eqnarray} \tag{ 22 }$

where ν_g is the representative frequency for the gth group. Our simple treatment of group mean coefficients is sufficient for the transport of continuum radiation. In particular, it is sufficient for neutrino transport in the core-collapse supernova problem because the monochromatic opacities and emissivities of neutrinos are smooth functions of neutrino frequency. It should, however, be emphasized that our algorithm for radiation hydrodynamics does not rely on the simple treatment of group mean coefficients used in this paper. We also assume that the flux limiter is approximately independent of frequency within each group, and denote the flux limiter and the scalar Eddington factor for the gth group as λ_g and f_g, respectively. As for the emissivity, we define

$\begin{equation} j_g = \frac{4\pi }{c}\eta (\nu _g) \Delta \nu _g, \end{equation} \tag{ 23 }$

where Δν_g = ν_{g + 1/2} − ν_{g − 1/2} is the group width. However, in the case of photon radiation, we have adopted the polylogarithm function-based method of Clark (1987) for computing the group-integrated Planck function. Thus, assuming local thermodynamic equilibrium, we have for photon radiation

$\begin{equation} j_g = \frac{4\pi }{c}\kappa _g B_g, \end{equation} \tag{ 24 }$

where B_g = ∫_gB_νdν.

Using these definitions and approximations, we obtain the multigroup radiation hydrodynamics equations

$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \boldsymbol{u}) = { } & 0, \end{eqnarray} \tag{ 25 }$

$\begin{eqnarray} \frac{\partial (\rho \boldsymbol{u})}{\partial t} + \nabla \cdot (\rho \boldsymbol{u} \boldsymbol{u}) + \nabla p + \sum _{g} \lambda _g \nabla E_g = { } & \boldsymbol{F}_G, \end{eqnarray} \tag{ 26 }$

$\begin{eqnarray} \frac{\partial (\rho E)}{\partial t} + \nabla \cdot (\rho E \boldsymbol{u} + p \boldsymbol{u}) + \boldsymbol{u} \cdot \sum _{g}\lambda _g \nabla E_g = { } & \sum _{g} c (\kappa _gE_{g}-j_g) + \boldsymbol{u}\cdot \boldsymbol{F}_G, \end{eqnarray} \tag{ 27 }$

$\begin{eqnarray} \frac{\partial (\rho Y_e)}{\partial t} + \nabla \cdot (\rho Y_e \boldsymbol{u}) = { } & \sum _{g} c \xi _g (\kappa _gE_{g}-j_g) , \end{eqnarray} \tag{ 28 }$

$\begin{eqnarray} \frac{\partial E_g}{\partial t} + \nabla \cdot \left(\frac{3-f_g}{2} E_g \boldsymbol{u}\right) - \boldsymbol{u} \cdot \nabla \left(\frac{1-f_g}{2} E_g\right) = { } & {-} c (\kappa _gE_{g}-j_g) + \nabla \cdot \left(\frac{c\lambda _g}{\chi _g} \nabla E_g \right), \nonumber\\ + \int _g \frac{\partial }{\partial \nu } \Bigg[\left(\frac{1-f}{2} \nabla \cdot \boldsymbol{u} + \frac{3f-1}{2} \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u} \right) & \nu E_\nu \Bigg] d\nu - \frac{3f_g-1}{2} E_g \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}, \end{eqnarray} \tag{ 29 }$

where g = 1, 2, ..., N, ξ_g = ξ(ν_g), ∑_g denotes ∑^N_{g = 1}.

2.3. Mathematical Characteristics of the Hyperbolic Subsystem

The system of multigroup radiation hydrodynamics (Equations (25)–(29)), similar to the system of gray radiation hydrodynamics we have studied in Paper II, has a hyperbolic subsystem,

$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \boldsymbol{u}) = { } & 0, \end{eqnarray} \tag{ 30 }$

$\begin{eqnarray} \frac{\partial (\rho \boldsymbol{u})}{\partial t} + \nabla \cdot (\rho \boldsymbol{u} \boldsymbol{u}) + \nabla p + \sum _{g} \lambda _g \nabla E_g = { } & \boldsymbol{F}_G, \end{eqnarray} \tag{ 31 }$

$\begin{eqnarray} \frac{\partial (\rho E)}{\partial t} + \nabla \cdot (\rho E \boldsymbol{u} + p \boldsymbol{u}) +\boldsymbol{u} \cdot \sum _{g} \lambda _g \nabla E_g = { } & \boldsymbol{u} \cdot \boldsymbol{F}_G, \end{eqnarray} \tag{ 32 }$

$\begin{eqnarray} \frac{\partial (\rho Y_e)}{\partial t} + \nabla \cdot (\rho Y_e \boldsymbol{u}) = { } & 0, \end{eqnarray} \tag{ 33 }$

$\begin{eqnarray} \frac{\partial E_g}{\partial t} + \nabla \cdot \left(\frac{3-f_g}{2} E_g \boldsymbol{u}\right) - \boldsymbol{u} \cdot \nabla \left(\frac{1-f_g}{2} E_g\right) = { } & 0, \end{eqnarray} \tag{ 34 }$

where g = 1, 2, ..., N. Since Y_e is just an auxiliary variable that does not alter the characteristic wave speeds of the hyperbolic subsystem, we will neglect it in the remainder of this subsection for simplicity. Gravitational force terms are source terms, and thus they are not included in the following analysis. The term $-[(3f_g-1)/2]E_g \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}$ in Equation (29) is not included in the analysis of the hyperbolic subsystem here even though it becomes comparable to the term $\nabla \cdot \lbrace [(3-f_g)/2] E_g \boldsymbol{u}\rbrace$ in Equation (34) in the streaming limit where f_g → 1. The first reason for not including the term is that the main purpose of including radiation in the hyperbolic subsystem is to improve the accuracy in the optically thick limit where the radiation can affect the hydrodynamics significantly. The term $-[(3f_g-1)/2]E_g \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}$ is not significant in the optically thick limit where f_g → 1/3. The second reason is that in the streaming limit, where λ → 0, radiation force has little impact on the motion of matter. Thus, radiation is essentially decoupled from the hydrodynamics of matter in the streaming limit. Moreover, in the streaming limit the diffusion term ∇ · [(cλ_g/χ_g)∇E_g] (Equation (29)) will dominate the term $-[(3f_g-1)/2]E_g \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}$ and all the terms in Equation (34). Therefore, the term $-[(3f_g-1)/2]E_g \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}$ is insignificant in both the diffusion limit and the streaming limit.

In CASTRO, we solve the hyperbolic subsystem with a Godunov method, which utilizes a characteristic-based Riemann solver (see Section 3.1 for details). The Godunov method requires that we analyze the mathematical characteristics of the hyperbolic subsystem. For simplicity, let us consider the system in one dimension, which can be written in terms of primitive variables as

$\begin{equation} \frac{\partial {Q}}{\partial {t}} + A \frac{\partial {Q}}{\partial {x}} = 0, \end{equation} \tag{ 35 }$

where the primitive variables are

$\begin{equation} Q = \left(\begin{array}{@{}c@{}}\rho \\ u \\ p \\ E_1 \\ E_2 \\ \vdots \\ E_N\end{array}\right), \end{equation} \tag{ 36 }$

and the Jacobian matrix is

$\begin{equation} A = \left(\begin{array}{@{}ccccccc@{}}u & \rho & 0 & 0 & 0 & \ldots & 0 \\ 0 & u & {1}/{\rho } & {\lambda _1}/{\rho } & {\lambda _2}/{\rho } & \ldots & {\lambda _N}/{\rho } \\ 0 & \gamma p & u & 0 & 0 & \ldots & 0 \\ 0 & E_1 (3-f_1)/2 & 0 & u & 0 & \ldots & 0 \\ 0 & E_2 (3-f_2)/2 & 0 & 0 & u & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & E_N (3-f_N)/2 & 0 & 0 & 0 & \ldots & u \end{array} \right), \end{equation} \tag{ 37 }$

where γ is the adiabatic index of the matter. This system is hyperbolic because the Jacobian matrix is diagonalizable with N + 3 real eigenvalues,

$\begin{equation} u-c_s, \ u+c_s, \ u, \ u, \ u, \ldots , \ u. \end{equation} \tag{ 38 }$

Here

$\begin{equation} c_s = \sqrt{\gamma \frac{p}{\rho } + \sum _{g}\left(\frac{3-f_g}{2}\right) \left(\frac{\lambda _g E_g}{\rho }\right)} = \sqrt{c_m^2 + \sum _{g} c_g^2} \end{equation} \tag{ 39 }$

is the radiation modified sound speed, where $c_m = \sqrt{\gamma p / \rho }$ is the sound speed without radiation and $c_g = \sqrt{[(3-f_g)/2] (\lambda _g E_g / \rho)}$ is the contribution to the sound speed from each group. The corresponding right eigenvectors are

$\begin{equation} \left(\begin{array}{@{}c@{}}1 \\ -c_s/\rho \\ \gamma p/\rho \\ c_1^2/\lambda _1 \\ c_2^2/\lambda _2 \\ \vdots \\ c_N^2/\lambda _N\end{array}\right), \ \left(\begin{array}{@{}c@{}}1 \\ c_s/\rho \\ \gamma p/\rho \\ c_1^2/\lambda _1 \\ c_2^2/\lambda _2 \\ \vdots \\ c_N^2/\lambda _N\end{array}\right), \ \left(\begin{array}{@{}c@{}}1 \\ 0 \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0\end{array} \right), \ \left(\begin{array}{@{}c@{}}0 \\ 0 \\ -\lambda _1 \\ 1 \\ 0 \\ \vdots \\ 0\end{array} \right), \ \left(\begin{array}{@{}c@{}}0 \\ 0 \\ -\lambda _2 \\ 0 \\ 1 \\ \vdots \\ 0\end{array} \right), \ \ldots , \left(\begin{array}{@{}c@{}}0 \\ 0 \\ -\lambda _N \\ 0 \\ 0 \\ \vdots \\ 1\end{array} \right), \ \end{equation} \tag{ 40 }$

and the corresponding left eigenvectors are

$\begin{equation} \begin{array}{rrrrrrr}(0,& -\rho /2c_s,& 1/2{c_s}^2,& \lambda _1/2 {c_s}^2,& \lambda _2/2 {c_s}^2,& \ldots ,& \lambda _N/2 {c_s}^2), \\ (0,& \rho /2c_s,& 1/2{c_s}^2,& \lambda _1/2 {c_s}^2,& \lambda _2/2 {c_s}^2,& \ldots ,& \lambda _N/2 {c_s}^2), \\ (1,& 0,& -1 / {c_s}^{2},& -\lambda _1 / {c_s}^2,& -\lambda _2 / {c_s}^2,& \ldots , & -\lambda _N / {c_s}^2), \\ (0,& 0,& -\omega _1/\lambda _1,& 1- \omega _1,& -\omega _1\lambda _2/\lambda _1, & \ldots ,& -\omega _1 \lambda _N/\lambda _1), \\ (0,& 0,& -\omega _2/\lambda _2,& -\omega _2\lambda _1/\lambda _2,& 1-\omega _2,& \ldots ,& -\omega _2\lambda _N/\lambda _2), \\ & & & \vdots & & &, \\ (0,& 0,& -\omega _N/\lambda _N,& - \omega _N\lambda _1/\lambda _N,& -\omega _N\lambda _2/\lambda _N,& \ldots ,& 1-\omega _N), \end{array} \end{equation} \tag{ 41 }$

where

$\begin{equation} \omega _g = \frac{c_g^2}{c_s^2}. \end{equation} \tag{ 42 }$

These N + 3 eigenvectors define the characteristic fields for the 1D system. By computing the product of the right eigenvectors and the gradients of their corresponding eigenvalues (LeVeque 2002), we find that the first and second fields are genuinely nonlinear corresponding to either a shock wave or a rarefaction wave. The rest of the fields are linearly degenerate corresponding to a contact discontinuity in either density or radiation energy density of one group moving at a speed of u. Note that there is also a jump in fluid pressure such that the total pressure, p_tot = p + ∑_gλ_gE_g, is constant across the contact discontinuity. Obviously, in three dimensions, there are two additional contact discontinuities for transverse velocities just like the case of pure hydrodynamics.

2.4. Advection in Frequency Space

The integrals in Equation (29) form the following system:

$\begin{equation} \frac{\partial E_g}{\partial t} = \int _g \frac{\partial }{\partial \nu } \left[\left(\frac{1-f}{2} \nabla \cdot \boldsymbol{u} + \frac{3f-1}{2} \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u} \right) \nu E_\nu \right] d\nu , \end{equation} \tag{ 43 }$

which can also be written in the conservation law form

$\begin{equation} \frac{\partial q_1}{\partial t} + \frac{\partial (a_qq_1)}{\partial \ln {\nu }} = 0, \end{equation} \tag{ 44 }$

where q₁ = νE_ν and $a_q = -[(1-f)/2] (\nabla \cdot \boldsymbol{u}) - [(3f-1)/2] \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}$ . Here we have used Equation (19). Equation (44) describes the advection of q₁ in logarithmic frequency space with an advection speed of a_q. Because ∫q₁dln ν is the integrated radiation energy, the system conserves energy.

Various terms in Equation (29) have been split between the hyperbolic subsystem (Section 2.3) and the system for the advection in the frequency space (Section 2.4). There is, however, another way to look at Equation (29). The evolution of the radiation energy density (Equation (29)) without the diffusion term ∇ · [(cλ_g/χ_g)∇E_g] and the source–sink term −c(κ_gE_g − j_g) can also be split as

$\begin{eqnarray} \frac{\partial E_g}{\partial t} = { } & -\nabla \cdot (E_g \boldsymbol{u}), \end{eqnarray} \tag{ 45 }$

$\begin{eqnarray} \frac{\partial E_\nu }{\partial t} = { } & \frac{\partial }{\partial \ln {\nu }} \left[ \left(\frac{1-f}{2} \nabla \cdot \boldsymbol{u} + \frac{3f-1}{2} \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}\right) E_\nu \right]. \end{eqnarray} \tag{ 46 }$

Equation (46) also represents an advection in logarithmic frequency space and can be written in the form

$\begin{equation} \frac{\partial q_2}{\partial t} + \frac{\partial (a_q q_2)}{\partial \ln {\nu }} = 0, \end{equation} \tag{ 47 }$

where q₂ = E_ν. Here the conserved quantity is the number of radiation particles (i.e., ∫E_ν/ν dν). This new way of splitting will be further discussed later in Section 3.1.

2.5. Multigroup Radiation Diffusion

The parabolic part of the system consists of the radiation diffusion and source–sink terms, which were omitted from the discussion of the hyperbolic subsystem (Section 2.3) and advection in frequency space (Section 2.4),

$\begin{eqnarray} \frac{\partial (\rho e)}{\partial t} = { } & \sum _{g} c (\kappa _gE_{g}-j_g), \end{eqnarray} \tag{ 48 }$

$\begin{eqnarray} \frac{\partial (\rho Y_e)}{\partial t} = { } & \sum _{g} c \xi _g (\kappa _gE_{g}-j_g), \end{eqnarray} \tag{ 49 }$

$\begin{eqnarray} \frac{\partial E_g}{\partial t} = { } & -c (\kappa _gE_{g}-j_g) + \nabla \cdot \left(\frac{c\lambda _g}{\chi _g} \nabla E_g \right), \end{eqnarray} \tag{ 50 }$

where e is the specific internal energy of the fluid, and g = 1, 2, ..., N. The term c(κ_gE_g − j_g) represents the energy exchange in the comoving frame between the matter and the gth radiation group through absorption and emission of radiation. Note that the terms on the right-hand side (RHS) of Equations (48)–(50) all contain the speed of light c. Thus, these are order O(1) terms in v/c. An implicit treatment for these equations is usually necessary because of their stiffness. For photon radiation, Equation (49) is not needed.

3. NUMERICAL METHODS

CASTRO uses a nested hierarchy of logically rectangular, variable-sized grids with simultaneous refinement in both space and time. The AMR algorithms for hyperbolic and radiation diffusion equations in CASTRO have been described in detail in Papers I and II, respectively. The multigroup solver does not require any changes in the AMR algorithms. Thus, we will not describe them again in this paper.

In this section, we describe the single-level integration scheme. For each step at a single level of refinement, the state is first evolved using an explicit method for the hyperbolic subsystem (Equations (30)–(34)) and advection in frequency space (Section 3.1). Then an implicit update for radiation diffusion and source–sink terms is performed (Section 3.2). In Section 3.2.4, we describe two acceleration schemes for speeding up the convergence rate of the implicit solver for multigroup radiation diffusion equations.

3.1. Explicit Solver for Hyperbolic Subsystem and Advection in Frequency Space

The hyperbolic subsystem (Section 2.3) is treated explicitly. Our Godunov algorithm for the hyperbolic subsystem of the multigroup radiation hydrodynamics is an extension of the hyperbolic solver for the gray radiation hydrodynamics presented in Paper II, which in turn is based on the hydrodynamics algorithm presented in Paper I of this series. We refer the reader to Papers I and II for a detailed description of the integration scheme, which supports a general equation of state, self-gravity, nuclear reactions, the advection of nuclear species and electron fraction. Here we will only present a brief overview of the Godunov scheme.

The time evolution of the hyperbolic subsystem can be written in the form

$\begin{equation} \frac{\partial \boldsymbol{U}}{\partial {t}} = - \nabla \cdot \boldsymbol{F}, \end{equation} \tag{ 51 }$

where $\boldsymbol{U} = (\rho , \rho \boldsymbol{u}, \rho E, E_1, E_2, \ldots , E_N)^{T}$ denotes the conserved variables and $\boldsymbol{F}$ is their flux. The conserved variables are defined at cell centers. We predict the primitive variables, including ρ, $\boldsymbol{u}$ , p, ρe, E_g, where g = 1, 2, ..., N, from cell centers at time tⁿ to edges at time t^{n + 1/2} and use an approximate Riemann solver to construct fluxes, $\boldsymbol{F}^{n+1/2}$ , on cell faces. This algorithm is formally second order in both space and time. The time step is computed using the standard Courant–Friedrichs–Lewy (CFL) condition for explicit methods, and the sound speed used in the computation is now the radiation-modified sound speed c_s (Equation (39)). The time step is also limited by maximally allowed changes in temperature and electron fraction in the implicit solver. Additional constraints for the time step are applied if additional physics, such as burning, is included. CASTRO solves the hyperbolic subsystem of radiation hydrodynamics with an unsplit piecewise parabolic method (PPM) with characteristic tracing and full corner coupling (Miller & Colella 2002). An approximate Riemann solver based on the Riemann solver of Bell et al. (1989) and P. Colella et al. (unpublished), which utilizes the eigenvectors of the system, is used to compute the Godunov states and fluxes at the cell interface.

Advection in frequency space (Section 2.4) can also be treated explicitly because its timescale is comparable to that of the hyperbolic subsystem. One approach is to adopt an operator-splitting method and to evolve Equation (44) with a Godunov method. This is the approach taken by van der Holst et al. (2011) in the CRASH code for photon radiation. The term $-[(3f_g-1)/2]E_\nu \hat{\boldsymbol{n}}\hat{\boldsymbol{n}} : \nabla \boldsymbol{u}$ can be treated as a source term of Equation (44). The success of this approach for a number of test problems has been demonstrated by van der Holst et al. (2011). We have also found that this approach can produce satisfactory results for all the photon test problems in Section 4. However, we have found that erroneous results were obtained for core-collapse supernova simulations (Section 4.2), probably due to operator-splitting errors. One difficulty in core-collapse supernova simulations is that the divergence of velocity can be significant because there is a strong shock and the typical velocity in front of the shock is ∼0.1–0.2 c. In the case of core-collapse supernovae, the negative divergence of velocity will shift neutrinos toward high frequency (Equation (44)). The stronger conservative scattering at higher frequency can further amplify the neutrino energy density due to the work term $\boldsymbol{u} \cdot \nabla \lbrace [(1-f_g)/2]E_g\rbrace$ (Equation (34)). These two processes can combine to cause unphysical results due to splitting errors. However, the origin of the various terms involved in these two processes is the term $\partial (\mathbf{\mathsf {P}}_\nu : \nabla \boldsymbol{u})/\partial \ln {\nu }$ in Equation (5). To avoid operator-splitting errors, it is therefore desirable to treat the term $\partial (\mathbf{\mathsf {P}}_\nu : \nabla \boldsymbol{u})/\partial \ln {\nu }$ without splitting it. On the other hand, the term can also make a significant contribution to the hyperbolic subsystem in the optically thick limit, and we do not wish to sacrifice accuracy of the hyperbolic subsystem by taking the traditional approach of leaving out radiation from the Riemann solver for the hydrodynamics.

To avoid the operator-splitting without sacrificing the accuracy of the hyperbolic subsystem, we have adopted a new approach based on splitting the radiation energy density equation into Equations (45) and (46). Note that the radiation diffusion and source–sink terms will be treated separately by an implicit solver (Section 3.2). Our new approach still utilizes the Godunov scheme for the hyperbolic subsystem governed by Equations (30)–(34). After the Godunov states and fluxes at the cell interfaces are obtained, all variables in the hyperbolic system except the radiation energy density are updated as usual. For the radiation energy density, two steps are taken. The radiation energy density is first updated according to Equation (45) with its RHS computed using the Godunov states. Then another Godunov solver is used to evolve Equation (47), which is equivalent to Equation (46). The conversion between E_g and E_ν is performed according to E_g = E_ν(ν_g)ν_gΔln ν_g. The Godunov states at time t^{n + 1/2} of the first Godunov solver are used in computing $\nabla \cdot \boldsymbol{u}$ and $\nabla \boldsymbol{u}$ to obtain the wave speed a_q in the frequency space advection. The second Godunov scheme for the frequency space advection is a standard approach based on the method of lines. For time integration we use the third-order total variation diminish Runge–Kutta scheme of Shu & Osher (1988). We use the piecewise linear method with a monotonized central slope limiter (van Leer 1977) to reconstruct q₂ and a_q in order to achieve high order in (logarithmic frequency) space. The reconstructed variables are used to compute the HLLE flux (Harten et al. 1983; Eisenstat 1981). The frequency space advection has its own CFL condition for stability. Since CASTRO uses the hydrodynamic CFL condition with other additional constraints for the time steps, if necessary, subcycling in time is employed in order to satisfy the frequency space CFL condition.

3.2. Implicit Solver for Radiation Diffusion and Source–Sink Terms

The implicit solver evolves the radiation and matter according to Equations (48)–(50) with the results of the explicit update (Section 3.1) as the initial conditions. Our scheme is based on a first-order backward Euler method.

3.2.1. Outer Newton Iteration

We solve the parabolic subsystem (Equations (48)–(50)) iteratively via Newton's method. We define

$\begin{eqnarray} F_e = { } & \rho e - \rho e^{-} - \Delta t \sum _g c (\kappa _g E_g - j_g), \end{eqnarray} \tag{ 52 }$

$\begin{eqnarray} F_Y = { } & \rho Y_e - \rho Y_e^{-} - \Delta t \sum _g c \xi _g (\kappa _g E_g - j_g), \end{eqnarray} \tag{ 53 }$

$\begin{eqnarray} F_g = { } & E_g - E_g^{-} - \Delta t\, \nabla \cdot \left(\frac{c\lambda _g}{\chi _g} \nabla E_g \right) + \Delta t\, c (\kappa _g E_g - j_g), \end{eqnarray} \tag{ 54 }$

where the − superscript denotes the state following the explicit update, and g = 1, 2, ..., N. The desired solution is for F_g, F_e, and F_Y to all be zero. We approach the solution by the Newton iteration, and in each iteration step we solve the following system:

$\begin{equation} \left[\begin{array}{ccc}({\partial F_e}/{\partial T})^{(k)} & ({\partial F_e}/{\partial Y_e})^{(k)} & ({\partial F_e}/{\partial E_g})^{(k)} \\ [3pt] ({\partial F_Y}/{\partial T})^{(k)} & ({\partial F_Y}/{\partial Y_e})^{(k)} & ({\partial F_Y}/{\partial E_g})^{(k)} \\ [3pt] ({\partial F_g}/{\partial T})^{(k)} & ({\partial F_g}/{\partial Y_e})^{(k)} & ({\partial F_g}/{\partial E_g})^{(k)} \end{array} \right] \left[\begin{array}{c}\delta T^{(k+1)}\\ \delta Y_e^{(k+1)}\\ \delta E_g^{(k+1)} \end{array}\right] = \left[\begin{array}{c}- F_e^{(k)}\\ - F_Y^{(k)}\\ - F_g^{(k)} \end{array}\right]. \end{equation} \tag{ 55 }$

Here, δT^{(k + 1)} = T^{(k + 1)} − T^(k), δY^{(k + 1)}_e = Y_e^{(k + 1)} − Y^(k)_e, and δE^{(k + 1)}_g = E_g^{(k + 1)} − E^(k)_g, where the (k) superscript denotes the stage of the Newton iteration. The exact Jacobian matrix in Equation (55) is very complicated. For simplicity, we neglect the derivatives of the diffusion coefficient cλ_g/χ_g. This is a common practice in radiation transport and does not appear to significantly degrade the convergence rate. With the approximate Jacobian matrix, we can eliminate δT and δY_e from Equation (55) and obtain the following equation for E^{(k + 1)}_g:

$\begin{eqnarray} \left(c\kappa _g + \frac{1}{\Delta t}\right) E_g^{(k+1)} &-& \nabla \cdot \left(\frac{c\lambda _g}{\chi _g} \nabla E_g^{(k+1)} \right) \\ = c j_g &+& \frac{E_g^{-}}{\Delta t}+ H_g \left[c\sum _{g^{\prime }}\left(\kappa _{g^{\prime }} E_{g^{\prime }}^{(k+1)} - j_{g^{\prime }} \right) - \frac{1}{\Delta t} (\rho e^{(k)} - \rho e^{-}) \right]\nonumber\\ &+& \Theta _g \left[c \sum _{g^{\prime }}\xi _{g^{\prime }}\left(\kappa _{g^{\prime }} E_{g^{\prime }}^{(k+1)} - j_{g^{\prime }} \right) - \frac{1}{\Delta t} \big(\rho Y_e^{(k)} - \rho Y_e^{-}\big) \right].\end{eqnarray} \tag{ 56 }$

To reduce clutter, we have dropped the (k) superscript for λ_g, κ_g, χ_g, and j_g. The newly introduced variables in Equation (56) are given by

$\begin{eqnarray} H_g = { } & \left(\frac{\partial j_g}{\partial T} - \frac{\partial \kappa _g}{\partial T} E_g^{(k)} \right) \eta _T - \left(\frac{\partial j_g}{\partial Y_e} - \frac{\partial \kappa _g}{\partial Y_e} E_g^{(k)}\right) \eta _Y, \end{eqnarray} \tag{ 57 }$

$\begin{eqnarray} \Theta _g = { } & \left(\frac{\partial j_g}{\partial Y_e} - \frac{\partial \kappa _g}{\partial Y_e} E_g^{(k)}\right) \theta _Y - \left(\frac{\partial j_g}{\partial T} - \frac{\partial \kappa _g}{\partial T} E_g^{(k)} \right) \theta _T, \end{eqnarray} \tag{ 58 }$

where

$\begin{eqnarray} \eta _T = { } & \frac{c\Delta t}{\Omega } \left[\rho + c \Delta t \sum _g \xi _g \left(\frac{\partial j_g}{\partial Y_e} - \frac{\partial \kappa _g}{\partial Y_e} E_g^{(k)}\right) \right], \end{eqnarray} \tag{ 59 }$

$\begin{eqnarray} \eta _Y = { } & \frac{c\Delta t}{\Omega } \left[c \Delta t \sum _g \xi _g \left(\frac{\partial j_g}{\partial T} - \frac{\partial \kappa _g}{\partial T} E_g^{(k)} \right) \right], \end{eqnarray} \tag{ 60 }$

$\begin{eqnarray} \theta _T = { } & \frac{c\Delta t}{\Omega } \left[ \rho \frac{\partial e}{\partial Y_e} + c \Delta t \sum _g \left(\frac{\partial j_g}{\partial Y_e} - \frac{\partial \kappa _g}{\partial Y_e} E_g^{(k)}\right) \right], \end{eqnarray} \tag{ 61 }$

$\begin{eqnarray} \theta _Y = { } & \frac{c\Delta t}{\Omega } \left[\rho \frac{\partial e}{\partial T} + c\Delta t \sum _g \left(\frac{\partial j_g}{\partial T} - \frac{\partial \kappa _g}{\partial T} E_g^{(k)} \right) \right], \end{eqnarray} \tag{ 62 }$

$\begin{eqnarray} \Omega = { } & \left[\rho \frac{\partial e}{\partial T} + c\Delta t \sum _g \left(\frac{\partial j_g}{\partial T} - \frac{\partial \kappa _g}{\partial T} E_g^{(k)} \right) \right] \left[\rho + c \Delta t \sum _g \xi _g \left(\frac{\partial j_g}{\partial Y_e} - \frac{\partial \kappa _g}{\partial Y_e} E_g^{(k)}\right) \right] \nonumber \\ & - \left[ \rho \frac{\partial e}{\partial Y_e} + c \Delta t \sum _g \left(\frac{\partial j_g}{\partial Y_e} - \frac{\partial \kappa _g}{\partial Y_e} E_g^{(k)}\right) \right] \left[c \Delta t \sum _g \xi _g \left(\frac{\partial j_g}{\partial T} - \frac{\partial \kappa _g}{\partial T} E_g^{(k)} \right) \right]. \end{eqnarray} \tag{ 63 }$

Here all the derivatives are computed from the results of the previous Newton iteration (i.e., the (k) state).

3.2.2. Inner Iteration

We now present our algorithm for solving Equation (56), which is actually a set of N equations coupled through the summation terms on the RHS of Equation (56). The system is usually too large to be solved directly. Instead, we solve it by decoupling the groups and utilizing an iterative method. Note that the iteration here for solving Equation (56) is different from the Newton iteration for solving the entire parabolic subsystem. The inner iteration here is embedded inside the outer Newton iteration step. In our algorithm, we avoid the coupling by replacing the radiation energy density in the coupling terms with its state from the previous inner iteration. Thus, each step of the inner iteration solves

$\begin{eqnarray} \left(c\kappa _g + \frac{1}{\Delta t}\right) E_g^{\ell +1} &-& \nabla \cdot \left(\frac{c\lambda _g}{\chi _g} \nabla E_g^{\ell +1} \right) \\ = c j_g &+& \frac{E_g^{-}}{\Delta t}+ H_g \left[c\sum _{g^{\prime }}\left(\kappa _{g^{\prime }} E_{g^{\prime }}^{\ell } - j_{g^{\prime }} \right) - \frac{1}{\Delta t} (\rho e^{(k)} - \rho e^{-}) \right] \nonumber\\ &+& \Theta _g \left[c \sum _{g^{\prime }}\xi _{g^{\prime }}\left(\kappa _{g^{\prime }} E_{g^{\prime }}^{\ell } - j_{g^{\prime }} \right) - \frac{1}{\Delta t} \big(\rho Y_e^{(k)} - \rho Y_e^{-}\big) \right], \end{eqnarray} \tag{ 64 }$

where ℓ is the inner iteration index. Here, we have dropped the (k + 1) superscript for E^{ℓ + 1}_g and E^ℓ_g to reduce clutter. For the first inner iteration, the state of the previous inner iteration is set to the result of the previous outer iteration. The inner iteration for Equation (64) is stopped when the maximum of ∑_g|(E^{ℓ + 1}_g − E^ℓ_g)|/∑_gE^{ℓ + 1}_g on the computational domain is below a preset tolerance (e.g., 10⁻⁶).

CASTRO can solve equations in the following canonical form (Paper II):

$\begin{equation} A E_g^{\ell +1} - \sum _i \frac{\partial }{\partial x^i} \left(B_i \frac{\partial E_g^{\ell +1}}{\partial x^i} \right) + \sum _i \frac{\partial }{\partial x^i} \left(C_i E_g^{\ell +1} \right) + \sum _i D_i \frac{\partial E_g^{\ell +1}}{\partial x^i} = \mathrm{RHS}, \end{equation} \tag{ 65 }$

where A are cell-centered coefficients and B_i, C_i, and D_i are centered at cell faces. We use the same approach for spatial discretization at the AMR coarse–fine interface as in Almgren et al. (1998). For Equation (64) the C_i and D_i coefficients are not used. The same canonical form works for Cartesian, cylindrical, and spherical coordinates so long as appropriate metric factors are included in the coefficients and the RHS. CASTRO uses the hypre library (Falgout & Yang 2002; hypre Code Project 2012) to solve the linear system.

3.2.3. Matter Update

When the numerical solution to Equation (56) is obtained via the inner iteration (Equation (64)), we update the matter and continue with the next outer iteration until convergence is achieved. Using Equation (55), we can obtain

$\begin{eqnarray} \delta T^{(k+1)} &=& { } \eta _T \left[\sum _g\big(\kappa _gE_g^{(k+1)} - j_g\big) - \frac{\rho e^{(k)} - \rho e^-}{c\Delta t} \right]\\ &&- \theta _T\left[ \sum _g\xi _g\big(\kappa _gE_g^{(k+1)} - j_g\big) - \frac{\rho Y_e^{(k)} - \rho Y_e^-}{c\Delta t} \right], \end{eqnarray} \tag{ 66 }$

$\begin{eqnarray} \delta Y_e^{(k+1)} &=& { } -\eta _Y \left[\sum _g\big(\kappa _gE_g^{(k+1)} - j_g\big) - \frac{\rho e^{(k)} - \rho e^-}{c\Delta t} \right]\\ &&+ \theta _Y \left[ \sum _g\xi _g\big(\kappa _gE_g^{(k+1)} - j_g\big) - \frac{\rho Y_e^{(k)} - \rho Y_e^-}{c\Delta t} \right]. \end{eqnarray} \tag{ 67 }$

Thus, we could update the matter temperature as T^{(k + 1)} = T^(k) + δT^{(k + 1)} and electron fraction as Y^{(k + 1)}_e = Y_e^(k) + δY^{(k + 1)}_e, and then compute the matter specific internal energy e^{(k + 1)} from T^{(k + 1)} and Y^{(k + 1)}_e. However, the total energy (radiation energy plus matter internal energy) and lepton number would not be conserved in this type of updates. Furthermore, our experience with this approach is that sometimes the update gives T and Y_e that are out of the bounds of opacity and EOS tables used in our core-collapse simulations and hence causes convergence issues. Thus, a different approach is taken in CASTRO. To conserve the total energy and lepton number, we update the matter as follows:

$\begin{eqnarray} \rho e^{(k+1)} &=& { } H \rho e^{(k)} + (1-H) \rho e^- + \Theta \big(\rho Y_e^{(k)} - \rho Y_e^-\big) \\ &&+ c \Delta t \sum _g\left[\big(\kappa _gE_g^{\ell +1} - j_g \big) - (H + \Theta \xi _g) \big(\kappa _gE_g^{\ell } - j_g\big)\right], \end{eqnarray} \tag{ 68 }$

$\begin{eqnarray} \rho Y_e^{(k+1)} &=& { }\bar{\Theta } \rho Y_e^{(k)} + (1-\bar{\Theta }) \rho Y_e^- + \bar{H} (\rho e^{(k)} - \rho e^-) \\ &&+ c\Delta t \sum _g\left[\xi _g\big(\kappa _gE_g^{\ell +1} - j_g\big) - (\bar{H} + \bar{\Theta } \xi _g) \big(\kappa _gE_g^{\ell } - j_g\big)\right], \end{eqnarray} \tag{ 69 }$

where H = ∑_gH_g, Θ = ∑_gΘ_g, $\bar{H} = \sum _g \xi _g H_g$ , and $\bar{\Theta } = \sum _g \xi _g \Theta _g$ . Then we compute temperature T^{(k + 1)} and electron fraction Y^{(k + 1)}_e from the newly updated ρe^{(k + 1)} and ρY^{(k + 1)}_e. We also compute the new absorption, scattering, and emission coefficients using the new temperature and electron fraction. By updating the matter using Equations (68) and (69), the implicit algorithm conserves the total energy and electron lepton number regardless of the accuracy of the solution to the parabolic subsystem.

The outer iteration for solving Equations (48)–(50) uses the following convergence criteria:

$\begin{eqnarray} |T^{(k+1)}-T^{(k)}| < {} & \epsilon \, T^{(k+1)}, \end{eqnarray} \tag{ 70 }$

$\begin{eqnarray} \big|Y_e^{(k+1)}-Y_e^{(k)}\big| < {} & \epsilon \, Y_e^{(k+1)}, \end{eqnarray} \tag{ 71 }$

$\begin{eqnarray} \big|F_e^{(k+1)}\big| < {} & \epsilon \, \rho \left[\left|\left(\frac{\partial e}{\partial T}\right)^{(k+1)}\right| |T^{(k+1)}-T^{(k)}| + \left|\left(\frac{\partial e}{\partial Y_e}\right)^{(k+1)}\right| \big|Y_e^{(k+1)}-Y_e^{(k)}\big| \right], \end{eqnarray} \tag{ 72 }$

$\begin{eqnarray} \big|F_Y^{(k+1)}\big| < {} & \epsilon \, \rho Y_e^{(k+1)}, \end{eqnarray} \tag{ 73 }$

where epsilon is a small number such as epsilon = 10⁻⁶. We do not use the relative error of the matter-specific internal energy e because it may not accurately reflect the error since that e can be shifted by an arbitrary amount. Typically, the outer iteration is stopped only when all four conditions (inequalities (70)–(73)) are satisfied. But for problems with extremely high opacities, inequalities (72) and (73) can be difficult to achieve, and only inequalities (70) and (71) are used.

3.2.4. Acceleration Scheme

We solve the multigroup diffusion equation (Equation (56)) using an iterative method. The system of the multigroup radiation diffusion equations is coupled. Each individual radiation group depends on the matter, which in turn depends on all radiation groups. When the coupling is strong (e.g., large opacity and small heat capacity) and/or the time step is large, the convergence rate of the iterative method we have presented in Section 3.2.2 can be very slow. To improve the convergence rate, we have extended the synthetic acceleration scheme of Morel et al. (1985, 2007) for photon radiation to neutrino radiation, and have developed two types of acceleration schemes.

Local acceleration scheme. We define the error of radiation energy density after ℓ inner iterations as

$\begin{equation} \epsilon _g^{\ell +1} = E_g^{e} - E_g^{\ell +1}, \end{equation} \tag{ 74 }$

where E^e_g is the exact solution of Equation (56) and E^{ℓ + 1}_g is the result of the ℓth inner iteration via Equation (64). In our local acceleration scheme, the spatial derivatives of epsilon ^{ℓ + 1}_g are assumed to be zero. This is justified because the spatially constant limit corresponds to the most slowly convergent mode (Morel et al. 1985) and our goal is to improve the convergence rate by curing the slowest decaying mode of the errors. Using Equations (56) and (64), we derive an equation for the error in the gth group,

$\begin{equation} \kappa _g^t \epsilon _g^{\ell +1} - H_g \sum _{g^{\prime }} \kappa _{g^{\prime }} \epsilon _{g^{\prime }}^{\ell +1} - \Theta _g \sum _{g^{\prime }} \xi _{g^\prime } \kappa _{g^{\prime }} \epsilon _{g^{\prime }}^{\ell +1} = H_g r_T^{\ell +1} + \Theta _g r_Y^{\ell +1}, \end{equation} \tag{ 75 }$

where

$\begin{eqnarray} \kappa _g^t = { } & \kappa _g + \frac{1}{c\Delta t}, \end{eqnarray} \tag{ 76 }$

$\begin{eqnarray} r_T^{\ell +1} = { } &\sum _{g^{\prime }} \kappa _{g^{\prime }} \big(E_{g^{\prime }}^{\ell +1} -E_{g^{\prime }}^{\ell }\big), \end{eqnarray} \tag{ 77 }$

$\begin{eqnarray} r_Y^{\ell +1} = { } &\sum _{g^{\prime }} \xi _{g^\prime } \kappa _{g^{\prime }} \big(E_{g^{\prime }}^{\ell +1} -E_{g^{\prime }}^{\ell }\big). \end{eqnarray} \tag{ 78 }$

Equation (75) can be solved analytically and gives

$\begin{equation} \epsilon _g^{\ell +1} = \frac{(qH_g + r \Theta _g) r_T^{\ell +1}+(p \Theta _g + sH_g) r_Y^{\ell +1}}{\kappa _g^t (pq-rs)}, \end{equation} \tag{ 79 }$

where

$\begin{eqnarray} p = { } & 1 - \sum _g \frac{H_g \kappa _g}{\kappa _g^t}, \end{eqnarray} \tag{ 80 }$

$\begin{eqnarray} q = { } & 1 - \sum _g \frac{\Theta _g \xi _g \kappa _g}{\kappa _g^t}, \end{eqnarray} \tag{ 81 }$

$\begin{eqnarray} r = { } & \sum \frac{H_g \xi _g \kappa _g}{\kappa _g^t}, \end{eqnarray} \tag{ 82 }$

$\begin{eqnarray} s = { } & \sum \frac{\Theta _g \kappa _g}{\kappa _g^t}. \end{eqnarray} \tag{ 83 }$

The improved solution is therefore

$\begin{equation} E_{g}^{\ell +1} := E_{g}^{\ell +1} + \epsilon _g^{\ell +1}. \end{equation} \tag{ 84 }$

The acceleration is placed after the convergence check for the inner iteration so that no unnecessary acceleration is applied to converged solutions.

Gray acceleration scheme. The original acceleration scheme of Morel et al. (1985) utilizes a gray diffusion equation for errors. It is postulated that the normalized multigroup spectrum of the errors is given by the normalized equilibrium spectrum. Note that our local acceleration scheme is based on the assumption of local equilibrium. Therefore, the spectrum defined as

$\begin{equation} \psi _g = \frac{{\epsilon }_g^{\ell +1}}{\sum _{g^{\prime }} {\epsilon }_{g^{\prime }}^{\ell +1}} \end{equation} \tag{ 85 }$

is postulated to be

$\begin{equation} \psi _g = \frac{\tilde{\epsilon }_g^{\ell +1}}{\sum _{g^{\prime }} \tilde{\epsilon }_{g^{\prime }}^{\ell +1}}, \end{equation} \tag{ 86 }$

where $\tilde{\epsilon }_g^{\ell +1}$ is the solution of the local acceleration scheme (Equation (79)). Substituting E^{ℓ + 1}_g = E_g^e − epsilon ^{ℓ + 1}_g into Equation (64) and summing over groups, we obtain

$\begin{equation} A \epsilon ^{\ell +1} - \nabla \cdot (\bar{d} \nabla \epsilon ^{\ell +1}) + \nabla \cdot (\boldsymbol{\phi } \epsilon ^{\ell +1}) = c H r_T^{\ell +1} + c \Theta r_Y^{\ell +1}, \end{equation} \tag{ 87 }$

where

$\begin{eqnarray} \epsilon ^{\ell +1} = {} & {\sum _{g} {\epsilon }_{g}^{\ell +1}}, \end{eqnarray} \tag{ 88 }$

$\begin{eqnarray} A = { } & c(1-H)\bar{\kappa } - c \Theta \bar{\kappa }_Y + \frac{1}{\Delta t}, \end{eqnarray} \tag{ 89 }$

$\begin{eqnarray} \bar{d} = { } & \sum _{g} \psi _g \frac{c\lambda _g}{\chi _g}, \end{eqnarray} \tag{ 90 }$

$\begin{eqnarray} \boldsymbol{\phi } = { } & - \sum _{g} \frac{c\lambda _g}{\chi _g} \nabla \psi _g, \end{eqnarray} \tag{ 91 }$

$\begin{eqnarray} \bar{\kappa } = {} & \sum _g \psi _g \kappa _g, \end{eqnarray} \tag{ 92 }$

$\begin{eqnarray} \bar{\kappa }_Y = {} & \sum _g \psi _g \xi _g \kappa _g. \end{eqnarray} \tag{ 93 }$

Equation (87) is in the canonical form that CASTRO solves (Equation (65)). The solution of the gray diffusion equation for error, epsilon ^{ℓ + 1}, is used to obtain the improved multigroup solution

$\begin{equation} E_{g}^{\ell +1} := E_{g}^{\ell +1} + \psi _g \epsilon ^{\ell +1}. \end{equation} \tag{ 94 }$

We note again that the acceleration operation is placed after the convergence check for the inner iteration so that no unnecessary acceleration is applied to converged solutions.

4. TEST PROBLEMS

In this section, we present detailed tests of the code demonstrating its ability to handle a wide range of radiation hydrodynamics problems. Photon radiation problems are presented in Section 4.1, whereas Section 4.2 shows neutrino radiation hydrodynamics simulations of core-collapse supernovae. Note that not every term or equation in the full system of radiation hydrodynamics is included in every test. This allows various parts of the code to be tested separately.

A CFL number of 0.8 is used for these tests unless stated otherwise or a fixed time step is used. For photon radiation test problems, no additional constraint is posed by the implicit solver. For core-collapse supernova simulations, the time step is further constrained by limiting the estimated relative changes in temperature to less than 1% and the estimated absolute changes in electron fraction to less than 0.01 using the change rates from the previous time step, and in our simulations this only happens around the time of core bounce. The relative tolerances for both the outer and inner iterations in the implicit update are 10⁻⁶. All convergence criteria for the outer iteration (inequalities (70)–(73)) are used except for the tests in Sections 4.1.4 and 4.1.7 with extremely high opacities. By default, the radiation groups are uniform in logarithmic frequency space and the representative frequency for the gth group is $\nu _g = \sqrt{\nu _{g-1/2}\nu _{g+1/2}}$ . For AMR simulations, the "deferred sync" scheme (Paper II) is used for flux synchronization between coarse and fine levels in the implicit solver.

4.1. Photon Radiation Test Problems

4.1.1. Linear Multigroup Diffusion

In this problem, we simulate a linear multigroup diffusion test problem introduced by Shestakov & Bolstad (2005). This can test the multigroup diffusion part of the system

$\begin{eqnarray} \frac{\partial (\rho e)}{\partial t} = { } & \sum _{g} c (\kappa _gE_{g}-j_g), \end{eqnarray} \tag{ 95 }$

$\begin{eqnarray} \frac{\partial E_g}{\partial t} = { } & {-}c (\kappa _gE_{g}-j_g) + \nabla \cdot \left(\frac{c\lambda _g}{\chi _g} \nabla E_g \right). \end{eqnarray} \tag{ 96 }$

To make the nonlinear system mathematically tractable, Shestakov & Bolstad (2005) linearize the system using the following assumptions. First, the absorption coefficient is assumed to be proportional to ν⁻³, and the group-averaged absorption is given by

$\begin{equation} \kappa _g = C_\kappa \nu _g^{-3}, \end{equation} \tag{ 97 }$

where C_κ is a constant. Second, radiation emission is assumed to obey Wien's law rather than Planck's law, and it is further modified so that

$\begin{equation} B_g = \frac{2\pi h}{c^2} \nu _g^3 \left[\exp \left(-\frac{h\nu _{g-1/2}}{k T_{f}}\right) - \exp \left(-\frac{h\nu _{g+1/2}}{k T_{f}}\right)\right] T, \end{equation} \tag{ 98 }$

where T_f is a fixed temperature. Note that the group-integrated B_g now has linear dependence on T. Finally, there is no scattering and the diffusion is not flux limited. Under these assumptions, Shestakov & Bolstad (2005) have obtained exact solutions for the linearized system.

The test that we have simulated is the first problem in Shestakov & Bolstad (2005). We have set up the test using physical units rather than the dimensionless units of Shestakov & Offner (2008). The 1D computational domain in physical units is (0, 1.0285926 × 10⁶) cm, corresponding to (0, 5.12) in dimensionless units. We use a symmetry boundary condition at the lower boundary, and a Marshak boundary with no incoming flux at the upper boundary (i.e., E_g + (2/3κ_g)∂E_g/∂x = 0). The mass density is 1.8212111 × 10⁻⁵ g cm⁻³. The specific heat capacity of the matter at constant volume is 9.9968637 × 10⁷ erg K⁻¹ g⁻¹. Initially, the temperature is set to 0.1 keV and 0, for x < 0.5 and x > 0.5, respectively. Here, x is a dimensionless coordinate. The radiation energy density is initially set to zero everywhere. The constant C_κ is set to 4.0628337 × 10⁴³ cm⁻¹ Hz³. The fixed temperature in the group-integrated radiation emission (Equation (98)) is T_f = 0.01 keV, corresponding to T_f = 0.1 in dimensionless units. We use 64 radiation groups with the lowest boundary at zero. The width of the first group is set to 1.2089946 × 10¹³ Hz, corresponding to 5 × 10⁻⁴ in dimensionless units, and the width of other groups is set to be 1.1 times the width of its immediately preceding group. The representative frequency for each group is set to $\nu _g = \sqrt{\nu _{g-1/2}\,\nu _{g+1/2}}$ , except that ν₁ is set to 0.5 ν_3/2 for the first group. We have run a simulation with 2048 uniform cells. Thus, the dimensionless cell size is 1/400. The time step is fixed at 5.8034112 × 10⁻⁸ s, corresponding to 1/200 in dimensionless units. No acceleration scheme is used in this test. Figure 1 shows the results after 200 steps. The comparison with the tabular data of Shestakov & Offner (2008) is shown in Table 1. Relative errors are computed as epsilon (f) = |f_n − f_e|/f_e, where f_n and f_e are the numerical and exact results, respectively. Because the tabular data of the exact solution are located on the faces of numerical cells, averaging of the numerical data has been performed for the comparison. The results are in excellent agreement with the exact solution and are comparable to the numerical results of Shestakov & Offner (2008).

**Figure 1.** Linear multigroup diffusion test. Numerical results are shown in lines, whereas the analytic solutions are shown in circle symbols. We show dimensionless total radiation energy density and temperature at the dimensionless time of 1 in panels (a) and (b), respectively. The spatial coordinate x is also dimensionless.
Download figure:
Standard image High-resolution image

Table 1. Relative Errors of the Linear Multigroup Diffusion Test Problem

x^a	(T) × 10³	(E_r) × 10³	x	(T) × 10³	(E_r) × 10³
0.00	0.0017	0.3209	0.51	4.8283	0.2585
0.20	0.0016	0.3220	0.52	1.8030	0.0232
0.40	0.0007	0.3460	0.53	1.0335	0.1495
0.46	0.0077	0.4095	0.54	0.7124	0.2330
0.47	0.0170	0.4443	0.60	0.3104	0.5469
0.48	0.0461	0.5136	0.80	0.5774	1.4067
0.49	0.2195	0.7169	1.00	1.3612	2.2413
0.50	0.0020	0.3712

Notes. The relative errors of temperature and total radiation energy density at the dimensionless time of 1 are shown. ^aDimensionless spatial coordinate.

Download table as: ASCII Typeset image

To assess the convergence behavior of our implicit scheme, we have performed a series of uniform-grid and AMR simulations of the linear multigroup diffusion test problem with various resolutions. The initial setup is the same as above except for spatial resolution and time step. Two levels are used in the AMR runs with the dimensionless domain of (0.25,0.75) covered by the fine (i.e., level 1) grids. In the convergence study, when the spatial resolution changes by a factor of two from one run to another, we change the time step by a factor of four. Because our implicit scheme is first order in time and second order in space, the expected convergence rate in the convergence study is second order with respect to Δx. Tables 2 and 3 show the L₁-norm errors and convergence rate for the total radiation energy density on the dimensionless spatial domain of (0, 1) at the dimensionless time of 1. Here, the L₁-norm error of the total radiation energy density is computed as

$\begin{equation} \epsilon _1(E_r) = \frac{ \sum _{i} |E_{r,i}-E_r(x_i)| \Delta x_{i} }{ \sum _{i} E_r(x_i) \Delta x_{i}}, \end{equation} \tag{ 99 }$

where E_{r, i} = ∑_gE_{g, i} is the numerical result for cell i and E_r(x_i) is the "exact" solution obtained from a high-resolution uniform-grid simulation with Δx = 1/3200 and Δt = 1/12800. The results of both uniform-grid and AMR simulations demonstrate the second-order convergence rate with respect to Δx as expected. Because only part of the domain in the AMR runs is covered by fine grids, the AMR runs have somewhat larger errors than the uniform-grid runs with the same resolution as the fine grids.

Table 2. L₁-norm Errors and Convergence Rate for the Total Radiation Energy Density at the Dimensionless Time of 1 in the Linear Multigroup Diffusion Test Problem

Δx	Δt	L₁ Error	Convergence Rate
1/100	2/25	2.91E−3
1/200	1/50	7.36E−4	2.0
1/400	1/200	1.90E−4	1.9
1/800	1/800	4.87E−5	2.0

Notes. Four one-dimensional uniform-grid runs with various resolutions are shown. The errors are computed on the dimensionless domain of (0, 1).

Download table as: ASCII Typeset image

Table 3. L₁-norm Errors and Convergence Rate for the Total Radiation Energy Density at the Dimensionless Time of 1 in the Linear Multigroup Diffusion Test Problem

Δx^a	Δt^b	L₁ Error	Convergence Rate
1/50, 1/100	4/25, 2/25	3.58E−3
1/100, 1/200	1/25, 1/50	1.05E−3	1.8
1/200, 1/400	1/100, 1/200	2.69E−4	2.0
1/400, 1/800	1/400, 1/800	6.82E−5	2.0

Notes. Four one-dimensional AMR runs with various resolutions are shown. The errors are computed on the dimensionless domain of (0, 1). ^aCell size on levels 0 and 1. ^bTime step on levels 0 and 1.

Download table as: ASCII Typeset image

4.1.2. Nonequilibrium Radiative Transfer with Picket-fence Model

Su & Olson (1999) have studied a nonequilibrium radiative transfer problem with a multigroup model, more specifically the picket-fence model. There are some special assumptions in this test problem. The volumetric heat capacity at constant volume is assumed to be c_v = αT³, where T is the matter temperature and α is a parameter. The group-integrated Planck function is assumed to be

$\begin{equation} B_g = p_g \left(\frac{a c}{4 \pi }\right) T^4, \end{equation} \tag{ 100 }$

where a is the radiation constant and p_g are parameters such that ∑_gp_g = 1. It is also assumed that the absorption coefficient is independent of temperature and there is no scattering. Su & Olson (1999) have also defined dimensionless variables for convenience. The dimensionless spatial coordinate is defined as

$\begin{equation} x = \bar{\kappa } z, \end{equation} \tag{ 101 }$

where z is the coordinate in physical units, and

$\begin{equation} \bar{\kappa } = \sum _{g}p_g \kappa _g. \end{equation} \tag{ 102 }$

The dimensionless time is defined as

$\begin{equation} \tau = \left(\frac{4 a c \bar{\kappa }}{\alpha }\right) t, \end{equation} \tag{ 103 }$

where t is the time in physical units. The dimensionless variables U_g for radiation energy density and V for matter energy density are defined as

$\begin{eqnarray} U_g &=& \frac{E_g}{a T_0^4}, \end{eqnarray} \tag{ 104 }$

$\begin{eqnarray} V &=& \left(\frac{T}{T_0}\right)^4, \end{eqnarray} \tag{ 105 }$

where T₀ is a reference temperature. In this test problem, there is a constant radiation source that exists for a finite period of time (0 ⩽ τ ⩽ τ₀) on a finite range (|x| < x₀). Su & Olson (1999) have found analytic solutions to the problem for both transport and diffusion approaches.

Here, we simulate Case C of Su & Olson (1999) by solving the multigroup diffusion equations (Equations (95) and (96)) and compare the results with the diffusion solution of Su & Olson (1999). Two radiation groups are used in this test, with p₁ = p₂ = 0.5. The absorption coefficients are set to κ₁ = 2/101 cm⁻¹ and κ₂ = 200/101 cm⁻¹. Thus, $\bar{\kappa } = 1\,\mathrm{cm}^{-1}$ . The constant α in the expression for volumetric heat capacity is set to α = 4a. The reference temperature is chosen to be T₀ = 10⁶ K. The parameters for the radiation source are x₀ = 0.5 and τ₀ = 10. When τ < τ₀, the radiation source deposits radiation energy to the region |x| < x₀ and increases the radiation energy density at a rate of $p_gc\bar{\kappa }aT_0^4$ for each group. The computational domain is 0 < x < 102.4 covered by 1024 uniform cells. The lower boundary is symmetric and the radiation energy density outside the upper boundary is set to zero. Initially, there is no radiation and the temperature is zero. A fixed time step Δτ = 0.1 is employed. The explicit solver is turned off in the simulation, and the radiation flux is not limited. No acceleration scheme is used in this test. The numerical results are shown in Figure 2. Our results are in excellent agreement with the analytic solution.

**Figure 2.** Nonequilibrium radiative transfer with the picket-fence model. Numerical results of U₁, U₂, and V are shown for τ = 3 (solid lines) and 30 (dashed lines). Also shown are analytic results for τ = 3 (circles) and 30 (diamonds).
Download figure:
Standard image High-resolution image

4.1.3. Radiating Sphere

Graziani (2008) found an analytic solution for the spectrum of a radiating sphere in a medium with frequency-dependent opacity. This test was used by Swesty & Myra (2009) to test their V2D code. In this test, a hot sphere with a fixed temperature is surrounded by a cold static medium. The absorption coefficient is assumed to be zero and the radiative transfer is due to coherent and isotropic scattering. Thus, there is no radiation–matter coupling or group coupling. The time-dependent spectrum at a given place in the medium is therefore the superposition of the blackbody spectrum of the medium and the spectrum of the radiation originated from the hot sphere.

We perform this test in 1D spherical geometry. The computational domain is 0.02 < r < 0.1 cm and is covered by 128 uniform cells. The hot sphere has a fixed temperature of 1.5 keV and the medium is at 50 eV. We use 60 energy groups covering an energy range of 0.5 eV to 308 keV. The group Rosseland mean coefficient is assumed to be χ_g = 10¹³(ν_g/3.6 × 10¹⁴ Hz)⁻³ cm⁻¹. The explicit solver is turned off. No acceleration scheme is used in this test. There is no outer iteration in the implicit solver because the matter properties do not change. The time step is fixed at 10⁻¹⁵ s. The spectrum at r = 0.06 cm and t = 10⁻¹² s is shown in Figure 3. The numerical results are again in good agreement with the analytic solution.

**Figure 3.** Spectrum from the radiating sphere observed at t = 10⁻¹² cm and r = 0.06 cm. Both numerical (circle symbols) and analytic (solid line) results are shown.
Download figure:
Standard image High-resolution image

4.1.4. Shock Tube Problem in Strong Equilibrium Regime

In this 1D problem, the initial conditions of the matter are ρ_L = 10⁻⁵ g cm⁻³, T_L = 1.5 × 10⁶ K, u_L = 0 and ρ_R = 10⁻⁵ g cm⁻³, T_R = 3 × 10⁵ K, u_R = 0, where L stands for the left half and R is the right half of the computational domain (0 < x < 100 cm). The gas is assumed to be ideal with an adiabatic index of γ = 4/3 and a mean molecular weight of μ = 1. Initially, the radiation is assumed to be in thermal equilibrium with the gas (i.e., E_g = 4πB_g/c). The interaction coefficients are set to κ_g = 10⁶ cm⁻¹ and χ_g = 10⁸ cm⁻¹. Thus, due to the huge opacities, the system is close to the limit of strong equilibrium with no diffusion and is essentially governed by Equations (30)–(34) with λ_g = 1/3 and f_g = 1/3. The parameters of this test problem are chosen such that neither radiation pressure nor gas pressure can be ignored. This test is adapted from the shock tube problem in Paper II. However, the adiabatic index of the gas is now γ = 4/3 rather than 5/3 so that we can compare the numerical results with the exact solution of a pure hydrodynamic Riemann problem. Furthermore, although this is a one-temperature gray problem due to the huge opacity, the numerical calculation is performed with 16 radiation groups with the lowest and uppermost boundaries at 10¹⁴ and 10¹⁸ Hz, respectively. We use 128 uniform zones in this test. The numerical results are shown in Figure 4. Also shown are the exact solutions of a pure hydrodynamic Riemann problem corresponding to the numerical problem. The results show good agreement. We also note that our scheme is stable without using small time steps even though the system is in the "dynamic diffusion" limit (Mihalas & Mihalas 1999) with τu/c ∼ 10⁷, where τ is the optical depth of the system. Either the gray acceleration scheme or the local acceleration scheme (Section 3.2.4) must be employed in this test, otherwise the implicit solver, which is responsible for the thermal equilibrium between radiation and matter, exhibits extremely slow convergence.

**Figure 4.** Shock tube at t = 10⁻⁶ s. Numerical results from a full radiation hydrodynamics calculation with 16 radiation groups and 128 cells are shown in symbols. The exact solutions of a pure hydrodynamic Riemann problem corresponding to the numerical problem are shown in solid lines for comparison. We show mass density (ρ), velocity (v), total pressure (p_tot), and radiation energy density (*E_r*). For the radiation hydrodynamics simulation, the total pressure is the sum of gas pressure and radiation pressure, and the radiation energy density is summed over all groups. In the pure hydrodynamic Riemann problem, there is no radiation energy density. All quantities are in cgs units.
Download figure:
Standard image High-resolution image

4.1.5. Nonequilibrium Radiative Shock

In this test, we simulate the Mach 5 shock problem of Lowrie & Edwards (2008). In Paper II, we showed results computed with the gray radiation hydrodynamics solver. Here, we present the numerical results obtained with the multigroup radiation hydrodynamics solver. In this test, the 1D numerical region (− 4000 cm, 2000 cm) initially consists of two constant states: ρ_L = 5.45969 × 10⁻¹³ g cm⁻³, T_L = 100 K, u_L = 5.88588 × 10⁵ cm s⁻¹ and ρ_R = 1.96435 × 10⁻¹² g cm⁻³, T_R = 855.720 K, u_R = 1.63592 × 10⁵ cm s⁻¹, where L stands for the left and R the right of the point x = 0, respectively. The gas is assumed to be ideal with an adiabatic index of γ = 5/3 and a mean molecular weight of μ = 1. The setup is the same as that in Paper II ⁴ except that 16 radiation groups are used in this paper. The lowest and uppermost boundaries of the groups are at 10¹⁰ and 10¹⁵ Hz, respectively. Initially, the radiation is assumed to be in thermal equilibrium with the gas (i.e., E_g = 4πB_g/c). The interaction coefficients are set to κ_g = 3.92664 × 10⁻⁵ cm⁻¹ and χ_g = 0.848903 cm⁻¹, respectively. The states at spatial boundaries are held at fixed values. Four refinement levels (five total levels) are used with a refinement factor of two with the finest resolution being Δx ≈ 1.5 cm. In this test, a cell is tagged for refinement if the second derivative of either density or temperature exceeds certain thresholds (Paper II). The initial conditions are set according to the pre-shock and post-shock states of the Mach 5 shock. After a brief period of adjustment, the system settles down to a steady shock. The simulation is stopped at t = 0.04 s. The results are shown in Figure 5. The agreement with the analytic solution including the narrow spike in temperature is excellent, except that the numerical results in the figure had to be shifted by −207 cm in space to match the analytic shock position. This discrepancy in shock position is due to the initial transient phase as the initial state relaxes to the correct steady-state profile. No flux limiter is used in this calculation because the analytic solution of Lowrie & Edwards (2008) assumes λ = 1/3. We again found that the implicit solver would fail to converge without an acceleration scheme. Both the gray acceleration scheme and the local acceleration scheme (Section 3.2.4) work in this test.

**Figure 5.** Density and temperature profiles for the Mach 5 supercritical shock. Numerical results are shown in symbols. We show density (circles), gas temperature (plus signs), and radiation temperature (squares). Here radiation temperature is defined as (*E_r*/a)^1/4, where *E_r* is the total radiation energy density. Also shown are the analytic results for density (red solid line), gas temperature (blue solid line), and radiation temperature (green solid line). We show only part of the region near the density jump. The inset shows a blow-up of the spike in temperature. The numerical results have been shifted by −207 cm in space to compensate for the discrepancy in shock position caused by the initial setup.
Download figure:
Standard image High-resolution image

4.1.6. Advecting Radiation Pulse

In Paper II, we performed a test of advecting a radiation pulse by the motion of matter with the gray radiation hydrodynamics solver (see also Krumholz et al. 2007). Here, we repeat the test with small changes to exercise the multigroup aspect of CASTRO. The setup is the same as that in Paper II, except that the absorption coefficient κ in this paper is given by

$\begin{equation} \kappa = 180 \left(\frac{T}{10^7\,\mathrm{K}}\right)^{-0.5} \left(\frac{\nu }{10^{18}\,\mathrm{Hz}}\right)^{-3} \left[1 - \exp \left(-\frac{h\nu }{kT}\right)\right] \,\mathrm{cm}^{-1}. \end{equation} \tag{ 106 }$

As in the test in Paper II, the initial temperature and density profiles are

$\begin{eqnarray} T &=& T_0 + (T_1-T_0) \exp {\left(-\frac{x^2}{2w^2}\right)}, \end{eqnarray} \tag{ 107 }$

$\begin{eqnarray} \rho &=& \rho _0\frac{T_0}{T_1} + \frac{a \mu }{3R} \left(\frac{T_0^4}{T} - T^3\right), \end{eqnarray} \tag{ 108 }$

where T₀ = 10⁷ K, T₁ = 2 × 10⁷ K, ρ₀ = 1.2 g cm⁻³, w = 24 cm, the mean molecular weight of the gas is μ = 2.33, and R is the ideal gas constant. Initially, the radiation is assumed to be in thermal equilibrium with the gas. No flux limiter is used in this test (i.e., λ = 1/3 and f = 1/3). If there were no radiation diffusion, the system would be in an equilibrium with the gas pressure and radiation pressure balancing each other. Because of radiation diffusion, the balance is lost and the gas moves.

The purpose of this test is twofold. First, we use this test to assess the ability of the code to handle radiation being advected by the motion of matter. We solve the problem numerically in different laboratory frames. Then we compare the case in which the gas is initially at rest with the case in which the gas initially moves at a constant velocity. Another purpose is to assess the effect of group resolution. Since the opacity is very large, the system is close to thermal equilibrium between radiation and matter. Hence, the results of multigroup radiation hydrodynamics simulations should become increasingly close to those of a gray radiation hydrodynamics simulation as more groups are used.

We have performed four runs: an unadvected run with 8 radiation groups, an advected run with 8 groups, an advected run with 64 groups, and an advected run using the gray radiation hydrodynamics solver presented in Paper II. The velocity in the unadvected run is initially zero everywhere, whereas in the advected runs it is 10⁶ cm s⁻¹ everywhere. For the multigroup runs, the lowest and uppermost group boundaries are at 10¹⁵ and 10¹⁹ Hz, respectively. In the gray simulation, the single group Planck mean coefficient, κ_P, and Rosseland mean coefficient, κ_R, of the opacity given by Equation (106) are

$\begin{eqnarray} \kappa _{\mathrm{P}} &=& 3063.96 \left(\frac{T}{10^7\,\mathrm{K}}\right)^{-3.5} \,\mathrm{cm}^{-1}, \end{eqnarray} \tag{ 109 }$

$\begin{eqnarray} \kappa _{\mathrm{R}} &=& 101.248 \left(\frac{T}{10^7\,\mathrm{K}}\right)^{-3.5} \,\mathrm{cm}^{-1}. \end{eqnarray} \tag{ 110 }$

The computational domain in all four runs is a 1D region of −512 < x < 512 cm with periodic boundaries. The numerical grid consists of 512 uniform cells. The simulations are stopped at t = 4.8 × 10⁻⁵ s. Figure 6 shows the density, velocity, and temperature profiles for all four runs. The results of the advected runs have been shifted in space for comparison. The profiles from the unadvected eight-group run and the advected eight-group run are almost identical, demonstrating the accuracy of our scheme in handling radiation advection by the motion of matter. This is expected because the velocity is significantly smaller than the speed of light (v/c ∼ 3 × 10⁻⁵). The results from the advected 64 group run and the gray run are also almost identical as expected. It is clear that the eight-group results differ from the gray results. This is not surprising because we did not perform group averaging for interaction coefficients (Equations (21) and (22)). We also show the relative differences in Figure 7. The relative difference of the advected 8 group run with respect to the unadvected 8 group run is computed as (advected − unadvected)/unadvected, and the relative difference of the 64 group run with respect to the gray run is computed as (multigroup − gray)/gray. The figure shows that the difference between the advected and unadvected 8 group runs is less than 0.03% everywhere, and the difference between the 64 group and gray runs is less than 0.1% everywhere.

**Figure 6.** Density, velocity and temperature profiles at t = 4.8 × 10⁻⁵ s for the test of an advecting radiation pulse. We show the results of four runs including the unadvected 8 group run (red solid lines), the advected 8 group run (blue dotted lines), the advected 64 group run (black dotted lines), and the advected gray run (green solid lines). The results of the advected runs have been shifted in space for comparison.
Download figure:
Standard image High-resolution image

**Figure 7.** Relative differences in density (solid lines) and gas temperature (dashed lines) in the test of an advecting radiation pulse. The relative differences between the advected and unadvected 8 group runs are shown in panel (a), whereas the relative differences between the advected 64 group and advected gray runs are shown in panel (b).
Download figure:
Standard image High-resolution image

4.1.7. Static Equilibrium

In Paper II, we have demonstrated that the gray solver can maintain a perfect static equilibrium in multiple dimensions because the radiation pressure and the gas pressure are coupled in the Riemann solver. Here, we repeat the test using the multigroup solver. The setup is the same as that in Paper II. The initial conditions for radiation and matter are the same as those in the test of the advecting radiation pulse (Section 4.1.6) except that the interaction coefficients are set to a very high value, κ_g = χ_g = 10²⁰ cm² g⁻¹ × ρ. Thus, almost no diffusion can happen, and the radiation and the matter are in perfect thermal equilibrium. We have performed a calculation on a 2D Cartesian grid of −512 < x < 512 cm and −512 < y < 512 cm with 512 uniform cells in each direction. The initial velocity is zero everywhere. For the initial setup, the coordinate x in Equation (107) is replaced by $r = \sqrt{x^2 + y^2}$ . We use 16 radiation groups with the lowest and uppermost boundaries at 10¹⁵ and 10¹⁹ Hz, respectively. Figure 8 shows the velocity profile at t = 10⁻⁴ s (after about 6000 steps). The maximal velocity at that time is ∼8 × 10⁻⁸ cm s⁻¹, which is about 10⁻¹⁵ of the sound speed. Such a small gas velocity indicates that the multigroup solver in CASTRO can also maintain a perfect static equilibrium in multiple dimensions, as does the gray solver in CASTRO.

**Figure 8.** Structure of the magnitude of velocity at t = 10⁻⁴ s in the static equilibrium test. The unit of velocity is cm s⁻¹.
Download figure:
Standard image High-resolution image

4.2. Core-collapse Supernova Simulations

In this section, we present core-collapse supernova simulations in 1D spherical and 2D cylindrical coordinates. Newtonian gravity with a monopole approximation is used in the 2D simulation for simplicity, although CASTRO has a multigrid Poisson solver for gravity. In the simulations, there are three neutrino species: electron neutrinos, ν_e, electron antineutrinos, $\bar{\nu }_e$ , and ν_μ, which stands for a combined species representing ν_μ, ν_τ, $\bar{\nu }_\mu$ , and $\bar{\nu }_\tau$ . We use a tabular equation of state that provides thermodynamic quantities as a function of temperature, density, and electron fraction and was constructed from the Shen et al. (1998a, 1998b) mean-field table, augmented with photons and a general electron/positron EOS, and extended down to 10 g cm⁻³. The opacities used are fully described in Burrows et al. (2006). They include (1) elastic scattering off nucleons, alpha particles, and the single representative heavy nucleus at the peak of the binding energy curve calculated by Shen et al. (1998a, 1998b), and (2) charged-current absorption cross sections off nucleons and nuclei (the latter using the approach of Bruenn 1985). The scattering and absorption cross sections off nucleons are corrected for weak magnetism and recoil effects, the form factor, electron screening, and ion–ion correlation effects on Freedman scattering off nuclei are applied, and pair production by nucleon–nucleon bremsstrahlung and electron–positron annihilation (and their inverses) are included as sources (and sinks). The sink terms of the latter do not include Fermi blocking corrections by the neutrinos in the final state. Inelastic neutrino–electron scattering is treated as elastic scattering.

4.2.1. One-dimensional Core-collapse Supernova Simulation

We present here a core-collapse supernova simulation in 1D spherical coordinates. The initial model is based on a 15 M_☉ progenitor of Woosley & Weaver (1995). The computational domain for this run is 0 < r < 5 × 10⁸ cm. The simulation employs four refinement levels (five total levels) with 1024 zones on the coarsest level (i.e., level 0) and a resolution of Δr ≈ 0.305 km on the finest level (i.e., level 4). The AMR refinement factor is two. The regions where the density is greater than 7.2 × 10⁶ g cm⁻³ or the enclosed mass is greater than 1.32 M_☉ are refined to level 2. In addition, cells on level ℓ will be refined to level ℓ + 1 if Y_e < 0.4, s > 5, ∇Y_e > 0.01/Δr_ℓ, or ∇s > 1/Δr_ℓ, where ℓ < 4, s is entropy in units of k baryon⁻¹, and Δr_ℓ is the cell width on level ℓ. We use 40 energy groups for electron neutrinos with the first group at 1 MeV and the last group at 300 MeV. We use 32 energy groups for electron antineutrinos with the first and last groups at 1 and 100 MeV, respectively. For ν_μ, we use 40 energy groups with the first and last groups at 1 and 300 MeV, respectively. The flux limiter of Levermore & Pomraning (1981) is used is this simulation. The local acceleration scheme is employed. The inner boundary is symmetric. An outflow boundary condition is used for the outer boundary in the explicit solver. The modified Marshak boundary of Sanchez & Pomraning (1991) with zero incoming flux is used for the outer boundary in the implicit solver. Thus, the region outside the outer boundary is essentially treated as a vacuum, and the radiation diffuses into the vacuum. Note that in the usual Marshak boundary based on the classical diffusion theory, the flux at the boundary does not have the correct behavior in the streaming limit, and the issue has been fixed in the modified Marshak boundary of Sanchez & Pomraning (1991) by using flux-limited diffusion.

The simulation is carried out to 200 ms after core bounce. Figure 9 shows four snapshots of density, velocity, temperature, entropy, electron fraction, and net electron lepton fraction (Y_L). The infall of the material due to gravity results in the increase of the central density until the core bounces. At bounce, the density reaches ρ^core ≈ 2.9 × 10¹⁴ g cm⁻³ at the core. A shock appears around bounce at M_sh ≈ 0.635 M_☉. The electron and lepton fractions at bounce are Y^core_e ≈ 0.288 and Y^core_L ≈ 0.326, respectively, at the core. The core temperature at bounce is T^core ≈ 9.98 MeV. The core entropy at bounce is s^core ≈ 0.901 k baryon⁻¹. These properties are in agreement with previous studies (e.g., Rampp & Janka 2000; Liebendörfer et al. 2001; Thompson et al. 2003). In a recent study, Lentz et al. (2012) have obtained M_sh ≈ 0.492 M_☉, ρ_core ≈ 4.264 × 10¹⁴ g cm⁻³, Y^core_e ≈ 0.2407, and Y^core_L ≈ 0.2782 for their model N-FullOp at bounce, and M_sh ≈ 0.717 M_☉, ρ_core ≈ 3.336 × 10¹⁴ g cm⁻³, Y^core_e ≈ 0.3046, and Y^core_L ≈ 0.3696 for their model N-ReduceOp at bounce. Except for ρ_core, the properties of our model at bounce are in between those in models N-FullOp and N-ReduceOp of Lentz et al. (2012).

Figure 9 shows a spike near the shock in the Y_e and Y_L profiles shortly after bounce (see also Lentz et al. 2012). The spike of Y_L is caused by electron neutrinos leaking out of the shock. The absorption of some of those neutrinos by the pre-shock material results in a spike of Y_e. In the post-shock region, the minimal Y_e drops rapidly for the first few milliseconds after bounce due to electron capture; it then becomes steady at ∼0.065–0.075. The velocity in some region behind the shock is positive for the first ∼2 ms after bounce; the maximum positive velocity reaches ∼2 × 10⁹ cm s⁻¹ similar to the maximum positive velocity in Thompson et al. (2003). The maximum temperature in the shocked region increases rapidly to 17.5 MeV at ∼0.1 ms after bounce; it then decreases quickly as the shock expands outward in radius; it increases again and reaches ∼22 MeV at 200 ms after bounce due to compression and deleptonization as matter settles onto the protoneutron star. The evolution of the temperature, electron fraction, and velocity profiles is in agreement with that in previous studies (e.g., Thompson et al. 2003).

In Figure 10, we show shock radius as a function of time. For the first few milliseconds after bounce, the shock expands very rapidly in radius and reaches ∼72 km. The peak shock radius reaches ∼140 km at ∼100 ms after bounce. After that the shock slowly shrinks in radius although it continues to grow in mass.

The comoving-frame luminosities of ν_e, $\bar{\nu }_e$ , and ν_μ neutrinos measured at 500 km are shown in Figure 11. In agreement with previous studies (e.g., Thompson et al. 2003; Marek & Janka 2009; Lentz et al. 2012), there is a small dip in $L_{\nu _e}$ after the initial rise. The dip is followed by a large pulse in $L_{\nu _e}$ that reaches 540 Bethe s⁻¹ at 4 ms after bounce. The FWHM for $L_{\nu _e}$ is ∼7 ms. The peak $L_{\nu _e}$ in our model is about 20% larger than that in Marek & Janka (2009) and Lentz et al. (2012). In our results, there is also a noticeable spike in $L_{\nu _\mu }$ . Note that Thompson et al. (2003) and Lentz et al. (2012) have obtained somewhat similar spikes in their simulations that did not include inelastic scattering. The lack of inelastic scattering of neutrinos on electrons results in harder neutrino spectra, especially for ν_μ, as shown in the rms average energies (Figure 12). Here the rms average is defined as

$\begin{equation} \langle \epsilon _\nu \rangle = \left[ \frac{\int \epsilon _\nu ^2 n_\nu d \nu }{\int n_\nu d \nu } \right]^{1/2}, \end{equation} \tag{ 111 }$

where epsilon _ν is the neutrino energy and n_ν is the monochromatic neutrino number density. Nevertheless, we obtain $\langle \epsilon _{\nu _e} \rangle < \langle \epsilon _{\bar{\nu }_e}\rangle < \langle \epsilon _{\nu _\mu } \rangle$ after the initial phases, as expected. At 200 ms after bounce, we obtain $\langle \epsilon _{\nu _e} \rangle \approx 14$ , $\langle \epsilon _{\bar{\nu }_e}\rangle \approx 19$ and $\langle \epsilon _{\nu _\mu } \rangle \approx 26\:\mathrm{MeV}$ . The luminosity spectra measured at 500 km in the comoving are shown Figure 13. The spectra are evidently harder than those in Thompson et al. (2003) because inelastic scattering is not included in our model.

**Figure 11.** Comoving-frame luminosities of ν_e (red solid line), $\bar{\nu }_e$ (green dotted line), and ν_μ neutrinos (blue dashed line) measured at 500 km for the 1D core-collapse supernova simulation. The luminosity of electron neutrinos reaches a peak value of 540 Bethe s⁻¹ at 4 ms after bounce. For $L_{\nu _e}$ , the FWHM is ∼7 ms. There is also a noticeable spike in the luminosity of ν_μ neutrinos with a peak value of 160 Bethe s⁻¹ at 17 ms after bounce.
Download figure:
Standard image High-resolution image

**Figure 11.** Comoving-frame luminosities of ν_e (red solid line), $\bar{\nu }_e$ (green dotted line), and ν_μ neutrinos (blue dashed line) measured at 500 km for the 1D core-collapse supernova simulation. The luminosity of electron neutrinos reaches a peak value of 540 Bethe s⁻¹ at 4 ms after bounce. For $L_{\nu _e}$ , the FWHM is ∼7 ms. There is also a noticeable spike in the luminosity of ν_μ neutrinos with a peak value of 160 Bethe s⁻¹ at 17 ms after bounce.
Download figure:
Standard image High-resolution image

**Figure 12.** Rms average energies of ν_e (red solid line), $\bar{\nu }_e$ (green dotted line), and ν_μ neutrinos (blue dashed line) measured at 500 km for the 1D core-collapse supernova simulation. For $\langle \epsilon _{\nu _\mu }\rangle$ , the peak value of 61 MeV is reached at 13 ms after bounce. The rms average energies at 200 ms after bounce are 14, 19, and 26 MeV for ν_e, $\bar{\nu }_{e}$ , and ν_μ, respectively. We do not show $\langle \epsilon _{\nu _\mu } \rangle$ before bounce because there is very little ν_μ at those early times.
Download figure:
Standard image High-resolution image

**Figure 13.** Luminosity spectrum of ν_e (red solid line), $\bar{\nu }_e$ (green dotted line), and ν_μ neutrinos (blue dashed line) measured at 500 km and 200 ms after bounce for the 1D core-collapse supernova simulation. The actual energy groups are shown in circles.
Download figure:
Standard image High-resolution image

4.2.2. Two-dimensional Core-collapse Supernova Simulation

Here we present a 2D core-collapse supernova simulation. The computational domain for this 2D cylindrical run is 0 < r < 2.5 × 10⁸ cm and −2.5 × 10⁸ cm < z < 2.5 × 10⁸ cm. Two refinement levels (three total levels) with a refinement factor of four are used. The coarsest level has 256 and 512 cells at the r- and z-directions, respectively. Thus, the size of the cells on the finest level is about 0.61 km. The same refinement criteria as those in the 1D simulation in Section 4.2.1 are used. We use 20, 16, and 20 groups with logarithmic spacing for ν_e, $\bar{\nu }_e$ , and ν_μ, respectively. The other aspects of this 2D run such as the outer boundary and flux limiter are the same as those in the 1D simulation in Section 4.2.1. We start the 2D simulation from the result of a 1D simulation rather than the 15 M_☉ presupernova model. The initial model is obtained from the 1D simulation at about 10 ms before bounce. At that time, multidimensional effects are still expected to be small. The 1D simulation that provides the initial model for the 2D run has similar resolution to the 2D run. Note that, other than resolution, the low-resolution 1D simulation here is essentially the same as the 1D high-resolution simulation in Section 4.2.1.

Figure 14 shows the density, entropy, temperature, and electron fraction profiles at 1.5 ms after bounce. The profiles are still highly symmetric. Spherically averaged profiles of density, velocity, temperature, entropy, electron fraction, and lepton fraction at 1.5 ms are shown in Figure 15. Also shown in Figure 15 are the results of the two 1D runs. At this point, the 2D results are only slightly different from the 1D results. However, in some region behind the shock, the gradient of entropy is negative. A negative entropy gradient region from ∼70 to ∼100 km is also clearly visible in the entropy profile at 4.1 ms (left panel of Figure 16). The white and black contour lines in the left panel of Figure 16 denote s = 7.5 and 5 k baryon⁻¹, respectively. Convective instabilities can potentially develop in the region between ∼60 and 90 km. In fact, convection does emerge quickly in our simulation. In the right panel of Figure 16, we show the entropy profile at 9.8 ms after bounce. Also shown in the figure is the velocity field. High entropy plumes rise while some low entropy material falls down. Secondary Kelvin–Helmholtz instabilities also develop in the shearing regions. The convection helps push the shock outward because of the enhanced material energy transport from the inner hot region to the region behind the shock by the convective motion. Figure 17 shows the density, entropy, temperature, and electron fraction profiles at 20 ms after bounce. The shock has moved to ∼200 km at 20 ms in the 2D simulation (see also Figure 18). At 40 ms after bounce, the north pole and the south pole of the shock have moved to ∼450 km. In contrast, in the 1D simulation in Section 4.2.1 (Figure 10), the shock radius is ∼85 km at 20 ms and it never exceeds 141 km. The comoving-frame luminosities of ν_e, $\bar{\nu }_e$ , and ν_μ neutrinos measured at 500 km for the 2D simulation are shown in Figure 19. Also shown in Figure 19 are the results of the 1D low-resolution run, which are nearly identical to those of the 1D high-resolution results (Figure 11). In the 2D results, there is also a large narrow pulse in $L_{\nu _e}$ that rises and drops at roughly the same rate. However, after ∼7 ms past bounce, the 2D results are qualitatively different from the 1D results because of the emergence of the convection. A major difference is that ν_μ neutrinos are much less luminous in two dimensions than in one dimension, because the temperature in two dimensions is cooler than that in one dimension due to the rapid expansion of 2D shock.

**Figure 15.** Density (upper left), radial velocity (lower left), temperature (upper center), entropy per baryon (lower center), electron fraction (upper right), and net lepton fraction (lower right) as a function of enclosed mass at 1.5 ms after bounce. Here, radial velocity is the velocity in the spherical radial direction. We show the results of the 2D cylindrical simulation (red solid lines), 1D low-resolution spherical run (blue dashed lines), and 1D high-resolution spherical run (green dotted lines). The radial profiles of the 2D results are generated via averaging.
Download figure:
Standard image High-resolution image

**Figure 16.** Entropy at 4.1 (left) and 9.8 ms (right) after bounce for the 2D cylindrical core-collapse supernova simulation. The white and black contour lines in the left panel denote s = 7.5 and 5 k baryon⁻¹, respectively. The arrows in the right panel show the velocity field in the shocked region at 9.8 ms after bounce. It is shown that convection has developed in the shocked region at 9.8 ms after bounce.
Download figure:
Standard image High-resolution image

**Figure 17.** Snapshot at 20 ms after bounce for the 2D cylindrical core-collapse supernova simulation. We show (a) logarithmic density, (b) entropy per baryon, (c) temperature, and (d) electron fraction, where the unit for density is g cm⁻³. Only the inner 270 km of the model is shown here.
Download figure:
Standard image High-resolution image

**Figure 18.** Shock radius vs. post-bounce time for the 2D core-collapse supernova simulation. We show the shock position in the northern hemisphere (red solid line), southern hemisphere (green dotted line), and equatorial plane (blue dashed lines) for the 2D simulation. Also shown is the result of the 1D low-resolution simulation (black dash-dotted line).
Download figure:
Standard image High-resolution image

**Figure 19.** Comoving-frame luminosities of ν_e, $\bar{\nu }_e$ , and ν_μ neutrinos measured at 500 km for the 2D cylindrical (red solid lines) and the low-resolution 1D spherical (black dashed lines) core-collapse supernova simulations.
Download figure:
Standard image High-resolution image

Since this paper is not on core-collapse supernova explosions per se, the 2D simulation is stopped at ∼55 ms after bounce. A detailed analysis of the 2D simulation is beyond the scope of this paper on algorithm, implementation, and testing.

5. PARALLEL PERFORMANCE

In Paper II, we demonstrated that the gray radiation solver in CASTRO has very good scaling behavior on up to 32,768 cores via a weak scaling study. The multigroup radiation solver presented in this paper has similar good weak scaling behavior because the same module for solving Equation (65) is used in both solvers. Here we present a strong scaling study that was carried out on Hopper, a petascale Cray XE6 supercomputer at the National Energy Research Scientific Computing Center. We have performed a series of 3D core-collapse supernova simulations in Cartesian coordinates with the same resolution on various numbers of cores. A hybrid MPI/OpenMP approach with six OpenMP threads per MPI task was used for these simulations. PFMG, a parallel multigrid solver from hypre (Falgout & Yang 2002; hypre Code Project 2012), was used to solve the linear systems. In each run, we start the 3D simulation from the result of a 1D simulation at about 10 ms before bounce. We use 12 energy groups for each species (ν_e, $\bar{\nu }_e$ , and ν_μ). The computational domain is a 3D box with a size of 5000 km in each dimension. Three refinement levels (four total levels) are used with a refinement factor of two between levels 0 and 1, and a refinement factor of four for other levels. The coarsest level has 256³ cells. The finest level, which approximately covers the inner region of $\sqrt{x^2+y^2+z^2} < 165\;\mathrm{km}$ , has a cell size of 0.61 km. Figure 20 shows the wall clock time per coarse time step for a series of runs with the number of cores ranging from 1536 to 24,576. Also shown in Figure 20 is the parallel efficiency defined as

$\begin{equation} E(N) = \frac{N_0 T_0}{N T_N}, \end{equation} \tag{ 112 }$

where N is the number of cores, T_N is the wall clock time for the N-core run, N₀ is the number of cores in a baseline run (i.e., 1536-core run), and T₀ is the wall clock time for the baseline run. In this strong scaling study, good scaling behavior is achieved on up to ∼10, 000 cores. However, because this is a strong scaling study, the workload per core diminishes as more cores are used. For example, the average number of cells for each core in the 24,576 core run is about 16³. An unfavorable consequence of small computational boxes is an increased communication overhead. Another unfavorable consequence is that more computation is wasted at domain boundaries due to larger surface area to volume ratios of smaller boxes. Hence, the degradation of parallel performance on more than 10,000 cores is not surprising.

**Figure 20.** Strong scaling behavior of 3D simulations on Hopper at NERSC. We show (a) wall clock time per coarse time step and (b) parallel efficiency in circles for five runs with the number of cores ranging from 1536 to 24,576. The dotted line in panel (a) indicates perfect scaling.
Download figure:
Standard image High-resolution image

6. SUMMARY

In this paper, we have presented a multidimensional multigroup radiation hydrodynamics solver that is part of the CASTRO code. Block-structured AMR is utilized in CASTRO. In this paper, we focus attention on the single-level algorithms because the AMR algorithms have been presented in Papers I and II and references therein.

Our multigroup radiation hydrodynamics solver is based on a comoving-frame formulation that is correct to order O(v/c) and uses an FLD approximation. Our mathematical analysis shows that the system we solve contains a hyperbolic subsystem that is the usual hydrodynamics system modified by radiation, a frequency space advection part that accounts for the Doppler shift of radiation due to the motion of matter, and a parabolic diffusion part. We have presented the mathematical characteristics of the hyperbolic subsystem. The eigenvalues and eigenvectors we have obtained are useful for characteristic-based Godunov schemes. We also described our treatment of the frequency space advection part. An implicit solver involving two nested iterations is presented. We have also presented two acceleration schemes that improve the convergence rate of the implicit solver.

Extensive testing is presented demonstrating the accuracy and stability of CASTRO to solve multigroup radiation hydrodynamics problems. We have presented a series of photon radiation test problems that cover a wide range of regimes. In addition to photon radiation problems, we have demonstrated the applicability of CASTRO in core-collapse supernova simulations.

Recently, Lentz et al. (2012) have argued, via a series of 1D simulations, that general relativistic effects and inelastic scattering of neutrinos on matter have significant effects on the numerical modeling of core-collapse supernova explosions. In contrast, Nordhaus et al. (2010) have concluded that spatial dimensionality has far more impact than general relativistic effects or microphysics such as inelastic scattering. Our 1D simulations are consistent with the assessment of Lentz et al. (2012) on the effects of inelastic scattering. However, our 2D simulation supports the assessment of Nordhaus et al. (2010) on the importance of spatial dimensionality. Our simulations have shown that multidimensional effects (e.g., convective motion in the region behind the shock and exterior to the nascent neutron star) appear less than 10 ms after bounce, and 2D results are profoundly different from 1D results. Thus, one should be cautious about the assessment of various effects such as inelastic scattering that are based on 1D simulations that last several hundred milliseconds past bounce like the ones presented in Lentz et al. (2012).

The main advantages of the CASTRO code are the efficiency due to the use of AMR and the accuracy due to the coupling of radiation force into the Riemann solver. Three-dimensional multigroup neutrino radiation hydrodynamics simulations of core-collapse supernovae with a good resolution using CASTRO can be carried out, if not now, in the near future with faster computers. In the strong scaling study with a total of 36 radiation groups and a cell size of 0.6 km on the finest level (Section 5), it took about 0.8 hr to advance one coarse time step (about 0.1 ms) on 12,288 cores of Hopper. Thus, it would take 800 hr on 12,288 cores of Hopper to evolve to 100 ms. Admittedly, it is still very expensive to perform long-duration 3D simulations, but it is possible to study the first 100 ms after bounce and investigate the development of the turbulent convection region via 3D simulations.

The work at LBNL was supported by the Office of High Energy Physics and the Office of Advanced Scientific Computing Research of the U.S. Department of Energy under contract No. DE-AC02-05CH11231. The work performed at LLNL was supported by the SciDAC program of the U.S. Department of Energy under the auspices of contract No. DE-AC52-07NA27344. Adam Burrows was supported by the SciDAC program of DOE under grant number DE-FG02-08ER41544, the NSF under subaward No. ND201387 to the Joint Institute for Nuclear Astrophysics, and the NSF PetaApps program, under award OCI-0905046 via a subaward No. 44592 from Louisiana State University to Princeton University. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC02-05CH11231. The authors would like to thank S.E. Woosley for many useful comments on the manuscript.

CASTRO: A NEW COMPRESSIBLE ASTROPHYSICAL SOLVER. III. MULTIGROUP RADIATION HYDRODYNAMICS

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. MULTIGROUP RADIATION HYDRODYNAMICS

2.1. Frequency-dependent Radiation Hydrodynamics Equations

2.2. Multigroup Radiation Hydrodynamics Equations

2.3. Mathematical Characteristics of the Hyperbolic Subsystem

2.4. Advection in Frequency Space

2.5. Multigroup Radiation Diffusion

3. NUMERICAL METHODS

3.1. Explicit Solver for Hyperbolic Subsystem and Advection in Frequency Space

3.2. Implicit Solver for Radiation Diffusion and Source–Sink Terms

3.2.1. Outer Newton Iteration

3.2.2. Inner Iteration

3.2.3. Matter Update

3.2.4. Acceleration Scheme

4. TEST PROBLEMS

4.1. Photon Radiation Test Problems

4.1.1. Linear Multigroup Diffusion

4.1.2. Nonequilibrium Radiative Transfer with Picket-fence Model

4.1.3. Radiating Sphere

4.1.4. Shock Tube Problem in Strong Equilibrium Regime

4.1.5. Nonequilibrium Radiative Shock

4.1.6. Advecting Radiation Pulse

4.1.7. Static Equilibrium

4.2. Core-collapse Supernova Simulations

4.2.1. One-dimensional Core-collapse Supernova Simulation

4.2.2. Two-dimensional Core-collapse Supernova Simulation

5. PARALLEL PERFORMANCE

6. SUMMARY

Footnotes

CASTRO: A NEW COMPRESSIBLE ASTROPHYSICAL SOLVER. III. MULTIGROUP RADIATION HYDRODYNAMICS

Article metrics

Permissions

Share this article

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. MULTIGROUP RADIATION HYDRODYNAMICS

2.1. Frequency-dependent Radiation Hydrodynamics Equations

2.2. Multigroup Radiation Hydrodynamics Equations

2.3. Mathematical Characteristics of the Hyperbolic Subsystem

2.4. Advection in Frequency Space

2.5. Multigroup Radiation Diffusion

3. NUMERICAL METHODS

3.1. Explicit Solver for Hyperbolic Subsystem and Advection in Frequency Space

3.2. Implicit Solver for Radiation Diffusion and Source–Sink Terms

3.2.1. Outer Newton Iteration

3.2.2. Inner Iteration

3.2.3. Matter Update

3.2.4. Acceleration Scheme

4. TEST PROBLEMS

4.1. Photon Radiation Test Problems

4.1.1. Linear Multigroup Diffusion

4.1.2. Nonequilibrium Radiative Transfer with Picket-fence Model

4.1.3. Radiating Sphere

4.1.4. Shock Tube Problem in Strong Equilibrium Regime

4.1.5. Nonequilibrium Radiative Shock

4.1.6. Advecting Radiation Pulse

4.1.7. Static Equilibrium

4.2. Core-collapse Supernova Simulations

4.2.1. One-dimensional Core-collapse Supernova Simulation

4.2.2. Two-dimensional Core-collapse Supernova Simulation

5. PARALLEL PERFORMANCE

6. SUMMARY

Footnotes