ABSTRACT
We present a formulation for multigroup radiation hydrodynamics that is correct to order O(v/c) using the comoving-frame approach and the flux-limited diffusion approximation. We describe a numerical algorithm for solving the system, implemented in the compressible astrophysics code, CASTRO. CASTRO uses a Eulerian grid with block-structured adaptive mesh refinement based on a nested hierarchy of logically rectangular variable-sized grids with simultaneous refinement in both space and time. In our multigroup radiation solver, the system is split into three parts: one part that couples the radiation and fluid in a hyperbolic subsystem, another part that advects the radiation in frequency space, and a parabolic part that evolves radiation diffusion and source–sink terms. The hyperbolic subsystem and the frequency space advection are solved explicitly with high-order Godunov schemes, whereas the parabolic part is solved implicitly with a first-order backward Euler method. Our multigroup radiation solver works for both neutrino and photon radiation.
Export citation and abstract BibTeX RIS
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
(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, ), respectively. Ye is the electron fraction defined to be the net electron number (electrons minus positrons) per baryon. FG is the external force on the fluid (e.g., gravitational force). Eν, , and 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
where mB is the baryon mass, h is the Planck constant, and s = 1 for νe neutrinos, s = −1 for 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 indicates contraction over two indices as follows:
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,
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
where
In the optically thick limit, R → 0, λ → 1/3, and the flux approaches the classical Eddington approximation . In the optically thin limit, R → ∞, λ ≈ 1/R, and as expected for streaming radiation. For the radiation pressure tensor, we use an approximate closure (Minerbo 1978; Levermore 1984)
where
is the scalar Eddington factor, is the identity tensor of rank 2, and . There are other choices for the scalar Eddington factor. For example, Kershaw (1976) has suggested
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
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
where νg − 1/2 and νg + 1/2 are the frequency at the lower and upper boundaries, respectively, for the gth group. For clarity, 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
where Bν = η/κ. In this paper, we will simply assume that
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 fg, respectively. As for the emissivity, we define
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
where Bg = ∫gBνdν.
Using these definitions and approximations, we obtain the multigroup radiation hydrodynamics equations
where g = 1, 2, ..., N, ξg = ξ(νg), ∑g denotes ∑Ng = 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,
where g = 1, 2, ..., N. Since Ye 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 in Equation (29) is not included in the analysis of the hyperbolic subsystem here even though it becomes comparable to the term in Equation (34) in the streaming limit where fg → 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 is not significant in the optically thick limit where fg → 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)∇Eg] (Equation (29)) will dominate the term and all the terms in Equation (34). Therefore, the term 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
where the primitive variables are
and the Jacobian matrix is
where γ is the adiabatic index of the matter. This system is hyperbolic because the Jacobian matrix is diagonalizable with N + 3 real eigenvalues,
Here
is the radiation modified sound speed, where is the sound speed without radiation and is the contribution to the sound speed from each group. The corresponding right eigenvectors are
and the corresponding left eigenvectors are
where
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, ptot = p + ∑gλgEg, 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:
which can also be written in the conservation law form
where q1 = νEν and . Here we have used Equation (19). Equation (44) describes the advection of q1 in logarithmic frequency space with an advection speed of aq. Because ∫q1dln ν 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)∇Eg] and the source–sink term −c(κgEg − jg) can also be split as
Equation (46) also represents an advection in logarithmic frequency space and can be written in the form
where q2 = 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),
where e is the specific internal energy of the fluid, and g = 1, 2, ..., N. The term c(κgEg − jg) 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
where denotes the conserved variables and is their flux. The conserved variables are defined at cell centers. We predict the primitive variables, including ρ, , p, ρe, Eg, where g = 1, 2, ..., N, from cell centers at time tn to edges at time tn + 1/2 and use an approximate Riemann solver to construct fluxes, , 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 cs (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 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 (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 in Equation (5). To avoid operator-splitting errors, it is therefore desirable to treat the term 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 Eg and Eν is performed according to Eg = Eν(νg)νgΔln νg. The Godunov states at time tn + 1/2 of the first Godunov solver are used in computing and to obtain the wave speed aq 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 q2 and aq 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
where the − superscript denotes the state following the explicit update, and g = 1, 2, ..., N. The desired solution is for Fg, Fe, and FY to all be zero. We approach the solution by the Newton iteration, and in each iteration step we solve the following system:
Here, δT(k + 1) = T(k + 1) − T(k), δY(k + 1)e = Ye(k + 1) − Y(k)e, and δE(k + 1)g = Eg(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 δYe from Equation (55) and obtain the following equation for E(k + 1)g:
To reduce clutter, we have dropped the (k) superscript for λg, κg, χg, and jg. The newly introduced variables in Equation (56) are given by
where
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
where ℓ is the inner iteration index. Here, we have dropped the (k + 1) superscript for Eℓ + 1g 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ℓ + 1g − Eℓg)|/∑gEℓ + 1g on the computational domain is below a preset tolerance (e.g., 10−6).
CASTRO can solve equations in the following canonical form (Paper II):
where A are cell-centered coefficients and Bi, Ci, and Di 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 Ci and Di 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
Thus, we could update the matter temperature as T(k + 1) = T(k) + δT(k + 1) and electron fraction as Y(k + 1)e = Ye(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 Ye 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:
where H = ∑gHg, Θ = ∑gΘg, , and . 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:
where is a small number such as = 10−6. 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
where Eeg is the exact solution of Equation (56) and Eℓ + 1g is the result of the ℓth inner iteration via Equation (64). In our local acceleration scheme, the spatial derivatives of ℓ + 1g 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,
where
Equation (75) can be solved analytically and gives
where
The improved solution is therefore
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
is postulated to be
where is the solution of the local acceleration scheme (Equation (79)). Substituting Eℓ + 1g = Ege − ℓ + 1g into Equation (64) and summing over groups, we obtain
where
Equation (87) is in the canonical form that CASTRO solves (Equation (65)). The solution of the gray diffusion equation for error, ℓ + 1, is used to obtain the improved multigroup solution
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−6. 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 . 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
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 ν−3, and the group-averaged absorption is given by
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
where Tf is a fixed temperature. Note that the group-integrated Bg 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 × 106) 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., Eg + (2/3κg)∂Eg/∂x = 0). The mass density is 1.8212111 × 10−5 g cm−3. The specific heat capacity of the matter at constant volume is 9.9968637 × 107 erg K−1 g−1. 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 × 1043 cm−1 Hz3. The fixed temperature in the group-integrated radiation emission (Equation (98)) is Tf = 0.01 keV, corresponding to Tf = 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 × 1013 Hz, corresponding to 5 × 10−4 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 , except that ν1 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−8 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 (f) = |fn − fe|/fe, where fn and fe 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).
Table 1. Relative Errors of the Linear Multigroup Diffusion Test Problem
xa | (T) × 103 | (Er) × 103 | x | (T) × 103 | (Er) × 103 |
---|---|---|---|---|---|
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: ASCIITypeset 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 L1-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 L1-norm error of the total radiation energy density is computed as
where Er, i = ∑gEg, i is the numerical result for cell i and Er(xi) 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. L1-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 | L1 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: ASCIITypeset image
Table 3. L1-norm Errors and Convergence Rate for the Total Radiation Energy Density at the Dimensionless Time of 1 in the Linear Multigroup Diffusion Test Problem
Δxa | Δtb | L1 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: ASCIITypeset 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 cv = αT3, where T is the matter temperature and α is a parameter. The group-integrated Planck function is assumed to be
where a is the radiation constant and pg are parameters such that ∑gpg = 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
where z is the coordinate in physical units, and
The dimensionless time is defined as
where t is the time in physical units. The dimensionless variables Ug for radiation energy density and V for matter energy density are defined as
where T0 is a reference temperature. In this test problem, there is a constant radiation source that exists for a finite period of time (0 ⩽ τ ⩽ τ0) on a finite range (|x| < x0). 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 p1 = p2 = 0.5. The absorption coefficients are set to κ1 = 2/101 cm−1 and κ2 = 200/101 cm−1. Thus, . The constant α in the expression for volumetric heat capacity is set to α = 4a. The reference temperature is chosen to be T0 = 106 K. The parameters for the radiation source are x0 = 0.5 and τ0 = 10. When τ < τ0, the radiation source deposits radiation energy to the region |x| < x0 and increases the radiation energy density at a rate of 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.
Download figure:
Standard image High-resolution image4.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 = 1013(νg/3.6 × 1014 Hz)−3 cm−1. 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−15 s. The spectrum at r = 0.06 cm and t = 10−12 s is shown in Figure 3. The numerical results are again in good agreement with the analytic solution.
Download figure:
Standard image High-resolution image4.1.4. Shock Tube Problem in Strong Equilibrium Regime
In this 1D problem, the initial conditions of the matter are ρL = 10−5 g cm−3, TL = 1.5 × 106 K, uL = 0 and ρR = 10−5 g cm−3, TR = 3 × 105 K, uR = 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., Eg = 4πBg/c). The interaction coefficients are set to κg = 106 cm−1 and χg = 108 cm−1. 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 fg = 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 1014 and 1018 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 ∼ 107, 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.
Download figure:
Standard image High-resolution image4.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−13 g cm−3, TL = 100 K, uL = 5.88588 × 105 cm s−1 and ρR = 1.96435 × 10−12 g cm−3, TR = 855.720 K, uR = 1.63592 × 105 cm s−1, 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 4 except that 16 radiation groups are used in this paper. The lowest and uppermost boundaries of the groups are at 1010 and 1015 Hz, respectively. Initially, the radiation is assumed to be in thermal equilibrium with the gas (i.e., Eg = 4πBg/c). The interaction coefficients are set to κg = 3.92664 × 10−5 cm−1 and χg = 0.848903 cm−1, 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.
Download figure:
Standard image High-resolution image4.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
As in the test in Paper II, the initial temperature and density profiles are
where T0 = 107 K, T1 = 2 × 107 K, ρ0 = 1.2 g cm−3, 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 106 cm s−1 everywhere. For the multigroup runs, the lowest and uppermost group boundaries are at 1015 and 1019 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
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−5 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−5). 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.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution image4.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 = 1020 cm2 g−1 × ρ. 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 . We use 16 radiation groups with the lowest and uppermost boundaries at 1015 and 1019 Hz, respectively. Figure 8 shows the velocity profile at t = 10−4 s (after about 6000 steps). The maximal velocity at that time is ∼8 × 10−8 cm s−1, which is about 10−15 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.
Download figure:
Standard image High-resolution image4.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, , and νμ, which stands for a combined species representing νμ, ντ, , and . 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−3. 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 × 108 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 × 106 g cm−3 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 Ye < 0.4, s > 5, ∇Ye > 0.01/Δrℓ, or ∇s > 1/Δrℓ, where ℓ < 4, s is entropy in units of k baryon−1, 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 (YL). 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 × 1014 g cm−3 at the core. A shock appears around bounce at Msh ≈ 0.635 M☉. The electron and lepton fractions at bounce are Ycoree ≈ 0.288 and YcoreL ≈ 0.326, respectively, at the core. The core temperature at bounce is Tcore ≈ 9.98 MeV. The core entropy at bounce is score ≈ 0.901 k baryon−1. 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 Msh ≈ 0.492 M☉, ρcore ≈ 4.264 × 1014 g cm−3, Ycoree ≈ 0.2407, and YcoreL ≈ 0.2782 for their model N-FullOp at bounce, and Msh ≈ 0.717 M☉, ρcore ≈ 3.336 × 1014 g cm−3, Ycoree ≈ 0.3046, and YcoreL ≈ 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).
Download figure:
Standard image High-resolution imageFigure 9 shows a spike near the shock in the Ye and YL profiles shortly after bounce (see also Lentz et al. 2012). The spike of YL 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 Ye. In the post-shock region, the minimal Ye 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 × 109 cm s−1 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.
Download figure:
Standard image High-resolution imageThe comoving-frame luminosities of ν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 after the initial rise. The dip is followed by a large pulse in that reaches 540 Bethe s−1 at 4 ms after bounce. The FWHM for is ∼7 ms. The peak 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 . 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
where ν is the neutrino energy and nν is the monochromatic neutrino number density. Nevertheless, we obtain after the initial phases, as expected. At 200 ms after bounce, we obtain , and . 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.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution image4.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 × 108 cm and −2.5 × 108 cm < z < 2.5 × 108 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, , 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−1, 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, , 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 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.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageSince 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, , 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 2563 cells. The finest level, which approximately covers the inner region of , 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
where N is the number of cores, TN is the wall clock time for the N-core run, N0 is the number of cores in a baseline run (i.e., 1536-core run), and T0 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 163. 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.
Download figure:
Standard image High-resolution image6. 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.
Footnotes
- 4
Unlike the simulation shown in Paper II, the simulation here used the 2010 CODATA recommended values of the fundamental physical constants. Therefore, the physical values of the states have changed slightly.