CRASH: A BLOCK-ADAPTIVE-MESH CODE FOR RADIATIVE SHOCK HYDRODYNAMICS—IMPLEMENTATION AND VERIFICATION

B. van der Holst; G. Tóth; I. V. Sokolov; K. G. Powell; J. P. Holloway; E. S. Myra; Q. Stout; M. L. Adams; J. E. Morel; S. Karni; B. Fryxell; R. P. Drake

doi:10.1088/0067-0049/194/2/23

1. INTRODUCTION

As photons travel through matter, the radiation field experiences changes due to net total emission, absorption, and scattering; see, for instance, Mihalas & Mihalas (1984), Pomraming (2005), and Drake (2006). At high enough energy density the radiation heats and accelerates the plasma. At a fundamental level, the radiation can be described by the time evolution of the spectral radiation intensity I_ν(r, t, n, ν), which is the radiation energy per unit area, per unit solid angle in the direction of photon propagation n, per unit interval of photon frequency ν, per unit time interval. Several methods have been developed to solve the radiation field in various degrees of physics fidelity.

In Monte Carlo radiative transfer methods, the radiation is statistically evaluated. Small photon packets are created with their energy and propagation direction statistically selected. The packets are propagated through matter using the radiation transfer equation (Nayakshin et al. 2009; Maselli et al. 2009; Baek et al. 2009). Characteristic methods use integration along rays of various lengths to solve for the angular structure of the radiation transport. A recent conservative, causal ray-tracing method was developed and combined with a short characteristic ray-tracing for the transfer calculations of ionizing radiation (Mellema et al. 2006). Solar surface magneto-convection simulations are increasingly realistic and use a three-dimensional (3D), non-gray, approximate local thermodynamic equilibrium (LTE), radiative transfer for the heating and cooling of plasma. These simulations are typically formulated for four frequency bins in the radiative transport equation (Vögler et al. 2005; Stein et al. 2007; Martínez-Sykora et al. 2009).

For some applications, simplifications to the radiation transfer can be made by calculating moments of the radiation intensity over the solid angle Ω. The spectral radiation energy and the spectral radiation energy flux are defined by the 0th and 1st moments as

$\begin{eqnarray} E_\nu ({\bf r},t,\nu) &=& \frac{1}{c}\int _{4\pi } I_\nu ({\bf r},t,{\bf n},\nu) d\Omega, \nonumber \\ {\bf F}_\nu ({\bf r},t,\nu) &=& \int _{4\pi } {\bf n} I_\nu ({\bf r},t,{\bf n},\nu) d\Omega. \end{eqnarray} \tag{ 1 }$

In addition, the spectral radiation pressure tensor P_ν is defined by the second moment

$\begin{equation} {\bf P}_\nu ({\bf r},t,\nu) = \frac{1}{c} \int _{4\pi } {\bf n}{\bf n} I_\nu ({\bf r},t,{\bf n},\nu)d\Omega. \end{equation} \tag{ 2 }$

A whole class of radiation transfer models is based on solving the corresponding radiation moment equations using a closure relation between the spectral pressure tensor and the spectral intensity (Mihalas & Mihalas 1984; Pomraming 2005; Drake 2006).

A radiation hydrodynamics code based on variable Eddington tensor (VET) methods (Stone et al. 1992) can still capture the angular structure of the radiation field by relating the spectral radiation pressure tensor to the spectral radiation intensity and the method is applicable for both the optically thin and thick regime. Optically thin versions of the VET method have been used in the context of cosmological reionization (Petkova & Springel 2009; Gneden et al. 2009). A more complete list of radiation hydrodynamic codes for the transfer of ionizing radiation in the context of astrophysical phenomena can be found in the radiative transfer comparison project (Iliev et al. 2009).

Further simplification assumes that the radiation pressure is isotropic and proportional to the radiation energy. This is the diffusion approximation. Several codes have been developed using this approximation. HYDRA (Marinak 2001) is an arbitrary Lagrange Eulerian code for two-dimensional (2D) and 3D radiation hydrodynamics. The radiation transfer model is based upon either flux-limited multi-group or implicit Monte Carlo radiation transport. The Eulerian code RAGE (Gittings et al. 2008) uses a cell-based adaptive mesh refinement (AMR) to achieve resolved radiative hydrodynamic flows. HYADES (Larsen & Lane 1994) solves the radiation hydrodynamic equations on a Lagrangian mesh, while CALE (Barton 1985) can use either a fixed Eulerian mesh, an embedded Lagrangian mesh, or a partially embedded, partially remapped mesh. Our newly developed radiation hydrodynamic solver uses an Eulerian grid together with a block-based AMR strategy.

We limit the discussion of the radiation hydrodynamics implementation in CRASH (Center for Radiative Shock Hydrodynamics) to plasmas in the absence of magnetic fields. Most of the description in this paper can, however, easily be extended to magnetohydrodynamic (MHD) plasmas as well. Indeed, since the CRASH code is essentially the MHD Block-Adaptive Tree Solarwind Roe Upwind Scheme (BATS-R-US) code (Powell et al. 1999; Tóth et al. 2011), extended with libraries containing radiation transport, equation-of-state (EOS), and opacity solvers, the implementation for the coupling between the radiation field and MHD plasmas is readily available. The CRASH code uses the recently developed Block-Adaptive Tree Library (BATL; Tóth et al. 2011). Here we will focus on the radiation implementation. Both the CRASH and BATS-R-US codes are publicly available as part of the Space Weather Modeling Framework (SWMF; Tóth et al. 2005) or can be used as stand-alone codes.

In the following, Section 2 introduces the radiation hydrodynamic equations for multi-material plasmas, in a form general enough to apply at high energy density. Section 3 describes the numerical algorithms to solve these equations. Section 4 presents verification tests for radiation and electron heat conduction on non-uniform meshes in one-dimensional (1D), 2D, and 3D slab geometry and in axially symmetric (rz) geometry. We also show a full system multi-material radiation hydrodynamic simulation on an adaptively refined mesh and demonstrate good scaling up to 1000 processors. The paper is summarized in Section 5.

2. EQUATIONS OF RADIATION HYDRODYNAMICS IN DENSE PLASMAS

The equations of radiation hydrodynamics describe the time evolution of both matter and radiation. For the applications that supported the work reported here, the code must be able to model matter as a high energy density plasma that is in LTE so that the population of all atomic and ion states can be obtained from statistical physics (see for instance Landau & Lifshitz 1980). We allow for multiple materials throughout the spatial domain of interest, but restrict the analysis to plasma flows that are far from relativistic. The materials can be heated to sufficiently high temperatures so that they can ionize and create free electrons, introducing the need for a time evolution equation for the electron energy density. The electrons transfer heat by thermal heat conduction and emit and absorb photon radiation. The radiation model discussed in this paper is non-equilibrium diffusion, in which the electron and radiation temperature can be different. We approximate the radiation transfer with a gray or multi-group flux-limited diffusion (FLD). This model is also of interest for application to a number of astrophysical problems.

In the following subsections, we will describe the radiative transfer equations for the evolution of the multi-group radiation energy densities (Section 2.1) in the FLD approximation (Section 2.5). The coupling of the radiation field to the two species hydrodynamic equations of electrons and ions are discussed in Section 2.2. In Section 2.3, the method for tracking the different materials is treated, while the lookup tables used for of the EOS and opacities are mentioned in Section 2.4.

2.1. Radiation Transport

In this section, we will build up the form of the radiation transport in the multi-group diffusion approximation that is used for the implementation in the CRASH code. The spectral pressure tensor, Equation (2), is often approximated in the form (Mihalas & Mihalas 1984)

$\begin{equation} {\bf P}_\nu ({\bf r},t,\nu) = E_\nu {\bf T}_\nu, \end{equation} \tag{ 3 }$

where

$\begin{equation} {\bf T}_\nu ({\bf r},t,\nu) = \frac{1}{2}(1-\chi _\nu){\bf I} + \frac{1}{2}(3\chi _\nu -1)\frac{{\bf F}_\nu {\bf F}_\nu }{|{\bf F}_\nu |^2}, \end{equation} \tag{ 4 }$

is the spectral Eddington tensor, χ_ν is the Eddington factor, and I is the identity matrix. The second term on the right-hand side is a dyad constructed from the direction of the spectral radiation flux. The pressure tensor can be used to arrive at a time evolution equation for the solid angle integrated spectral radiation energy (Buchler 1983)

$\begin{eqnarray} \frac{\partial E_\nu }{\partial t}+\nabla \cdot (E_\nu \mathbf {u}) -\nu \frac{\partial }{\partial \nu } ({\bf P}_\nu :\nabla {\bf u}) &=& {\rm diffusion + emission} \nonumber \\ && - {\rm absorption}, \end{eqnarray} \tag{ 5 }$

which contains the velocity u of the background plasma. Here the colon denotes the contraction of the two tensors P_ν and ∇u. The processes described by the symbolic terms on the right-hand side of Equation (5) will be described below.

Setting the Eddington factor χ_ν = 1/3 corresponds to the radiation diffusion model. In this case the radiation is assumed to be effectively isotropic and the spectral radiation pressure can be described by the scalar

$\begin{equation} p_\nu = \frac{1}{3}E_\nu =(\gamma _r - 1)E_\nu, \end{equation} \tag{ 6 }$

where we have introduced the adiabatic index of the radiation field, which in this case has the relativistic value γ_r = 4/3. The time evolution for the spectral energy density can then be simplified to

$\begin{eqnarray} \frac{\partial E_\nu }{\partial t}&+&\nabla \cdot (E_\nu \mathbf {u}) -(\gamma _r-1)(\mathbf {\nabla } \cdot \mathbf {u})\nu \frac{\partial E_\nu }{\partial \nu } \nonumber \\ &=& {\rm diffusion + emission - absorption}. \end{eqnarray} \tag{ 7 }$

The second and third terms on the left-hand side of Equation (7) express the change in the spectral energy density due to the advection and compression of the background plasma, which moves with the velocity u, as well as the frequency shift due to compression. In the free-streaming limit where the radiation hardly interacts with the matter, χ_ν approaches 1. In this paper, we will keep χ_ν = 1/3 and at the same time use an FLD for the free-streaming regime whenever needed (see Section 2.5).

The set of equations for the spectral energy density (7) still consists of an infinite number of equations, one for each frequency. A finite set of governing equations to describe the radiation transport in the multi-group diffusion approximation is obtained when we choose a set of frequency groups. Here we enumerate groups with the index, g = 1, ..., G. The interval of the photon frequencies, relating to the gth group is denoted as [ν_g−1/2, ν_g+1/2]. A discrete set of group energy densities, E_g, is introduced in terms of the integrals of the spectral energy density of the frequency group interval:

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

Now we can integrate Equation (7) to arrive at the desired set of the multi-group equations:

$\begin{eqnarray} \frac{\partial E_g}{\partial t} &+& \nabla \cdot (E_g \mathbf {u})+ (\gamma _r-1)E_g \mathbf {\nabla } \cdot \mathbf {u} - (\gamma _r-1)(\mathbf {\nabla } \cdot \mathbf {u})\nonumber\\ &\times&\int _{\nu _{g-1/2}}^{\nu _{g+1/2}} \frac{\partial (\nu E_\nu)}{\partial \nu }d\nu \nonumber \\ &=& \int _{\nu _{g-1/2}}^{\nu _{g+1/2}} {\rm (diffusion + emission - absorption)}\,d\nu.\quad \end{eqnarray} \tag{ 9 }$

The fourth term on the left-hand side is a frequency shift due to the plasma compression. This term is essentially a conservative advection along the frequency axis.

In the context of the multi-group radiation diffusion, a discussion about the stimulated emission is not less important than LTE. Excellent accounts on the stimulated emission exist in the literature, see for instance Zel'dovich & Raizer (2002). Here we merely summarize how the stimulated emission modifies the absorption opacity κ^a_ν obtained from, e.g., opacity tables. This is important when dealing with externally supplied opacity tables, since the CRASH code assumes that the absorption opacities are corrected. Integrating the total absorption and emission over all directions and summing up the two polarizations of the photons, the following expression can be derived for the emission and absorption

$\begin{equation} {\rm emission} - {\rm absorption} = c{\kappa _\nu ^a}^{\prime } (B_\nu - E_\nu), \end{equation} \tag{ 10 }$

where the effective absorption coefficient, κ^a_ν', is introduced to account for the correction due to stimulated emission:

$\begin{equation} {\kappa _\nu ^a}^{\prime } = \kappa _\nu ^a \left(1 - \exp \left[ -\frac{\varepsilon }{k_B T_e} \right] \right), \end{equation} \tag{ 11 }$

in which ε = hν is the photon energy, k_B is the Boltzmann constant, and T_e is the electron temperature. We also introduced the spectral energy density distribution of the blackbody radiation (the Planckian)

$\begin{equation} B_\nu = \frac{8\pi }{h^3 c^3} \frac{\varepsilon ^3}{\exp [ \varepsilon / (k_B T_e)] - 1}. \end{equation} \tag{ 12 }$

The total energy density in the Planck spectrum equals B = ∫^∞₀dνB_ν = aT⁴_e, where a = 8π⁵k⁴_B/(15h³c³) is the radiation constant.

We use the standard definition of the group Planck mean opacity κ_Pg and group Rosseland mean opacity κ_Rg (Mihalas & Mihalas 1984)

$\begin{eqnarray} \kappa _{Pg} &=& \frac{\int _{\nu _{g-1/2}}^{\nu _{g+1/2}}d\nu {\kappa _\nu ^a}^{\prime }B_\nu }{B_g}, \qquad \kappa _{Rg} = \frac{\frac{\partial B_g}{\partial T_e}}{\int _{\nu _{g-1/2}}^{\nu _{g+1/2}}d\nu \frac{1}{\kappa _\nu ^t} \frac{\partial B_\nu }{\partial T_e}}, \nonumber\\ B_g &=& \int _{\nu _{g-1/2}}^{\nu _{g+1/2}} d\nu B_\nu \end{eqnarray} \tag{ 13 }$

in which κ^t_ν is the spectral total opacity. The right-hand side of Equation (9) can now be written as (see for instance Mihalas & Mihalas 1984; Pomraming 2005)

$\begin{eqnarray} \frac{\partial E_g}{\partial t} &+& \nabla \cdot (E_g \mathbf {u}) + (\gamma _r - 1)E_g \mathbf {\nabla } \cdot \mathbf {u} - (\gamma _r - 1) \mathbf {\nabla } \cdot \mathbf {u}\nonumber\\ &\times& \int _{\nu _{g-1/2}}^{\nu _{g+1/2}} \frac{\partial (\nu E_\nu)}{\partial \nu }d\nu \nonumber \\ &=& \nabla \cdot (D_g\nabla E_g) + \sigma _g (B_g-E_g), \end{eqnarray} \tag{ 14 }$

where D_g = c/(3κ_Rg) is the radiation diffusion coefficient for radiation group g in the diffusion limit. The absorption and emission uses the coefficient σ_g = cκ_Pg. These group mean opacities are either supplied by lookup tables or by an opacity solver.

In a single group approximation (gray diffusion), the spectral energy density is integrated over all photon frequencies and the total radiation energy density is obtained by

$\begin{equation} E_r({\bf r},t) = \int _0^\infty {E_\nu d\nu }. \end{equation} \tag{ 15 }$

This amounts to summing up all groups E_r = ∑_gE_g. The gray radiation diffusion equation can be derived as (see for instance Mihalas & Mihalas 1984; Pomraming 2005; Drake 2006)

$\begin{equation} \frac{\partial E_r}{\partial t} + \nabla \cdot (E_r \mathbf {u}) + (\gamma _r - 1)E_r \mathbf {\nabla } \cdot \mathbf {u} = \nabla \cdot (D_r\nabla E_r) + \sigma _r (B-E_r), \end{equation} \tag{ 16 }$

where the diffusion coefficient D_r in the diffusion limit is now defined by the single group Rosseland mean opacity κ_R as D_r = c/(3κ_R), and the absorption coefficient σ_r is defined in terms of the single group Planck mean opacity κ_P as σ_r = cκ_P.

2.2. Hydrodynamics

In the CRASH code, a single fluid description is used so that all of the atomic and ionic species as well as the electrons move with the same bulk velocity u. The conservation of mass

$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \mathbf {u}) = 0 \end{equation} \tag{ 17 }$

provides the time evolution of the mass density ρ of all the materials in the simulation. The plasma velocity is obtained from the conservation of momentum

$\begin{equation} \frac{\partial \rho \mathbf {u} }{\partial t} + \nabla \cdot [ \rho \mathbf {u} \mathbf {u} + I(p+p_r)] = 0. \end{equation} \tag{ 18 }$

The total plasma pressure is the sum of the ion and electron pressures: p = p_i + p_e. The net force of the radiation on the plasma is given by the gradient of the total radiation pressure −∇p_r, where the total radiation pressure is obtained from the group radiation energies: p_r = (γ_r − 1)∑E_g.

In a high density plasma, the electrons are very strongly coupled to the ions by collisions. However, for higher temperatures, the electrons and ions get increasingly decoupled. At a shock front, where ions are preferentially heated by the shock wave, the electrons and ions are no longer in temperature equilibrium. Ion energy is transferred to the electrons by collisions, while the electrons in turn radiate energy. We therefore solve separate equations for the ion/atomic internal energy density E_i and the electron internal energy density E_e:

$\begin{equation} \frac{\partial E_i}{\partial t} + \nabla \cdot (E_i \mathbf {u})+p_i \mathbf {\nabla } \cdot \mathbf {u} = \sigma _{ie} (T_e - T_i), \end{equation} \tag{ 19 }$

$\begin{eqnarray} \frac{\partial E_e}{\partial t} + \nabla \cdot (E_e \mathbf {u})+p_e \mathbf {\nabla } \cdot \mathbf {u} &= &\nabla \cdot (C_e\nabla T_e) + \sigma _{ie} (T_i - T_e) \nonumber \\ &+& \sum _{g=1}^G \sigma _g (E_g-B_g). \end{eqnarray} \tag{ 20 }$

The coupling coefficient σ_ie = n_ak_B/τ_ie in the collisional energy exchange between the electrons and ions depends on the ion–electron relaxation time τ_ie(T_e, n_a, m) and the atomic number density n_a. Energy transfer depends also on the difference between ion temperature T_i and electron temperature T_e. In Equation (20), we have included electron thermal heat conduction with conductivity C_e(T_e, n_a, m). Since electrons are the species that are responsible for radiation emission and absorption, the energy exchange between the electrons and the radiation groups is added to Equation (20).

For the development of the numerical schemes in Section 3, we will use an equation for the conservation of the total energy density

$\begin{equation} e = \frac{\rho u^2}{2} + E_i + E_e + \sum _{g=1}^G E_g, \end{equation} \tag{ 21 }$

instead of the equation for ion internal energy (19). This is especially important in regions of the computational domain where hydrodynamic shocks can occur, so that we can recover the correct jump conditions. Conservation of total energy can be derived from Equations (14) and (17)–(20) as

$\begin{equation} \frac{\partial e}{\partial t} + \nabla \cdot [ (e+p+p_r) \mathbf {u} ] = \nabla \cdot (C_e\nabla T_e) + \sum _{g=1}^G \nabla \cdot (D_g\nabla E_g). \end{equation} \tag{ 22 }$

The frequency shift term in Equation (14), due to the plasma compression, does not show up in the conservation of total energy if we use energy conserving boundary conditions at the end points of the frequency domain, i.e., at ν = 0 and ν = ∞ in the analytical description or at the end points of the numerically truncated finite domain.

2.3. Level Sets and Material Identification

In many CRASH applications, we need a procedure to distinguish between different materials that may be present. We assume that materials do not mix, but differ from each other by their properties such as the EOS and opacities. If we use M different materials, then we can define for each material m = 1, ..., M the level set function d_m(r, t) (see for instance Kreiss & Olsson 2005; Olsson et al. 2007; Sussman & Pucket 2000) which is initially set to zero at the material interface, positive inside the material region, and negative outside the material region. Generally, we use a smooth and signed distance function in the initial state. At later times, the location of material m is obtained by means of a simple advection equation

$\begin{equation} \frac{\partial d_m}{\partial t} + \nabla \cdot (d_m \mathbf {u}) = d_m \mathbf {\nabla } \cdot \mathbf {u}. \end{equation} \tag{ 23 }$

For any given point in space and time, we can determine what the material is since analytically only one of the level set functions d_m can be positive at any given point. Numerical errors will create regions where this is not true in solutions to the discretized form of Equation (23). Whenever this happens, the material having the largest d_m is assigned as the material in that region. This is a simple approximation, and we may explore more sophisticated approaches in the future. The number of material levels M is configured at compile time.

2.4. Equation of State and Opacities

We have implemented EOS solvers and a code to calculate the frequency-averaged group opacities. This implementation will be reported elsewhere, but it is important to note that in the EOS and opacity solver, the temperature is assumed to be well below relativistic values: T ≪ 10⁵ eV. A non-relativistic speed of motion is also assumed, which simplifies the radiation transport equation and allows neglect of relativistic corrections for opacities. In this paper, we will assume that all necessary quantities are calculated and stored in lookup tables. Our EOS solver assumes that the corrections associated with ionization, excitation, and Coulomb interactions of partially ionized ion–electron plasma are all added to the energy of the electron gas and to the electron pressure. This is possible since these corrections are controlled by the electron temperature. The ion internal energy density, ion pressure, and ion specific heat in an isochoric process per unit of volume are simply

$\begin{equation} p_i = n_ak_BT_i, \;\; E_i = \frac{p_i}{\gamma -1}, \;\; C_{Vi} = \left(\frac{\partial E_i}{\partial T_i} \right)_\rho = \frac{n_ak_B}{\gamma -1}, \end{equation} \tag{ 24 }$

which are due to the contributions from ion translational motions, for which γ = 5/3.

The relations among electron internal energy density, pressure, density, and temperature are known as the EOS. To solve these relations is usually complex and time consuming. We therefore store these relations in invertible lookup tables. For each material m, our EOS tables have logarithmic lookup arguments (log T_e, log n_a). The list of quantities stored in these tables is indicated in Table 1. These lookup tables are populated with quantities that are needed for both single- and two-temperature simulations. For two-temperature plasma simulations, we need p_e, E_e, the electron specific heat C_Ve, and the electron-speed-of-sound gamma $\gamma _{S_e}$ . For convenience, we add the total pressure p = p_e + p_i, total internal energy density E = E_i + E_e, single-temperature specific heat C_V, and the single-temperature speed-of-sound gamma γ_S, which can be used in single-temperature simulations. We use high enough table resolutions so that it is sufficient to use a bilinear interpolation in the lookup arguments. If p_e or E_e (or p and E in single temperature mode) are known upon entry to the lookup instead of T_e, we do a binary search in the table to find the appropriate electron temperature. The latter works only if the necessary thermodynamic derivatives are sign definite, i.e., the table is invertible. Other thermodynamic quantities that are needed, but not stored in these lookup tables, can be derived. For example, the electron density can be obtained from the mean ionization $n_e = n_a \overline{Z}$ .

Table 1. Quantities Stored in the EOS Tables as a Function of log T_e[eV] and log n_a[m⁻³]

Quantity	Stored Quantity	Units
Total pressure p	p/n_a	eV
Total internal energy density E	E/n_a	eV
Electron pressure p_e	p_e/n_a	eV
Electron internal energy density E_e	E_e/n_a	eV
Specific heat C_V	C_V/(n_ak_B)
Electron specific heat C_Ve	C_Ve/(n_ak_B)
Speed of sound gamma	γ_S
Electron speed of sound gamma	$\gamma _{S_e}$
Inverse of ion–electron interaction time	1/τ_ie	s⁻¹
Electron conductivity	C_e	J m⁻¹ s⁻¹ K⁻¹
Mean ionization	$\overline{Z}$
Mean square ionization	$\overline{Z^2}$

Download table as: ASCII Typeset image

In addition, we have lookup tables for the averaged multi-group opacities. These tables are either constructed internally for a given frequency range, number of groups, and selected materials or externally supplied. For any material m, the logarithmic lookup arguments are (log ρ, log T_e). The stored quantities (see Table 2) are the specific Rosseland mean opacity κ_Rg/ρ and specific Planck mean opacity κ_Pg/ρ for all groups g = 1, ..., G used during a simulation. Planck opacities are assumed to be corrected for stimulated emission, as discussed in Section 2.1. The groups are always assumed to be logarithmically distributed in frequency space.

Table 2. Quantities Stored in the Opacity Tables as a Function of log ρ[kgm⁻³] and log T_e[eV]

Quantity	Symbol	Units
Specific group Rosseland mean opacities	κ_Rg/ρ	kg⁻¹ m²
Specific group Planck mean opacities	κ_Pg/ρ	kg⁻¹ m²

Download table as: ASCII Typeset image

2.5. Flux-limited Diffusion

Radiation diffusion, if uncorrected, can transport energy too efficiently in the optically thin free-streaming limit. In the diffusion limit, the radiation diffusion flux for each group follows Fick's law F_g = −D_g∇E_g, where the diffusion coefficient D_g depends on the Rosseland mean opacity κ_Rg for the group g via D_g = c/(3κ_Rg). However, this flux is potentially unbounded. In the optically thin free-streaming limit, the magnitude of the radiation flux can be at most cE_g in order to maintain causality. Various flux limiters exist in the literature (see, for instance, Minerbo 1978; Lund & Wilson 1980; and Levermore & Pomraming 1981) that ensure that the diffusion flux is limited by this free-streaming flux. We have implemented the so-called square-root flux limiter to obtain the correct propagation speed in the optically thin regime (Morel 2000). For this flux limiter, the diffusion coefficient is rewritten as

$\begin{equation} D_g = \frac{c}{\sqrt{(3\kappa _{Rg})^2 + \frac{|\nabla E_g|^2}{E_g^2}}}. \end{equation} \tag{ 25 }$

In the limit that the radiation length scale L_R = E_g/|∇E_g| is large, the diffusive limit is recovered. For a small radiation length scale, D_g = c|E_g|/|∇E_g| and the radiation propagates with the speed of light.

Similarly, we have implemented the option to limit the electron thermal heat flux (see Drake 2006 for more details on electron flux limiters). The classical Spitzer–Harm formula for the collisional electron conductivity is proportional to $T_e^{5/2}/\overline{Z^2}$ , where $\overline{Z^2}$ is the mean square ionization of the used material. The collisional model is only valid when the temperature scale length L_T = T_e/|∇T_e| is much larger than the collisional mean free path of the electrons λ_mfp. When the temperature scale length is only a few λ_mfp or smaller, this description breaks down. This may for instance happen in laser-irradiated plasmas. In such a case, one could determine the heat flux by solving the Fokker–Planck equation for the electrons, but this is computationally expensive. Instead, we use a simplified model to limit the electron heat flux. A free-streaming heat flux can be defined as the thermal energy density in the plasma transported at some characteristic thermal velocity: F_FS = n_ek_BT_ev_th, where $v_{th}=\sqrt{k_BT_e/m_e}$ . For practical applications, the maximum heat transport is usually only a fraction of this free-streaming flux: F = −(fF_FS/|∇T_e|)∇Te, where f is the so-called flux limiter. This heat flux model is the threshold model and is also used in other radiation hydrodynamics packages, such as HYADES (Larsen & Lane 1994). The flux-limited heat flux can now be defined as

$\begin{equation} {\bf F} = -\min \left(C_e,\frac{fF_{FS}}{|\nabla T_e|}\right)\nabla T_e. \end{equation} \tag{ 26 }$

The flux limiter f is an adjustable input parameter and can be tuned to let the simulated results better fit reality.

3. THE NUMERICAL METHOD

In this section, we describe the discretization of the set of multi-material, radiation hydrodynamic equations for the density (17), momentum (18), total energy (22), electron internal energies (20), radiation group energy (14), and material level-set functions (23). The equations are time integrated using an operator-split method to solve the equations in substeps. Formally, we may write this system as

$\begin{equation} \frac{\partial {\bf U}}{\partial t} = {\bf R}_{{\rm hydro}}({\bf U}) + {\bf R}_{{\rm frequency}}({\bf U}) + {\bf R}_{{\rm diffusion}}({\bf U}), \end{equation} \tag{ 27 }$

where U is the vector of state variables. We have split the right-hand sides of the equations into three parts and time advance the equations with an operator splitting method in the following order. (1) The right-hand side R_hydro describes the advection and pressure contributions (Section 3.1). This part is essentially the ideal hydrodynamic equations augmented with the advection and compression of the radiation energy, the electron internal energy, and the material level sets. (2) The right-hand side R_frequency is the advection of the radiation field in frequency space (Section 3.2). (3) The right-hand side R_diffusion takes care of the diffusion and energy exchange terms, which we will solve with an implicit scheme (Section 3.3). This choice of operator splitting is not unique. Instead of splitting the hydrodynamic advection operator and the extra advance operator for the frequency advection, one could attempt to discretize the frequency advection flux as an extra flux for the control volume of the four-dimensional (x, y, z, ν) space. However, since the CRASH code is built around the existing BATS-R-US code in 1D, 2D, and 3D, we opted for splitting the frequency advection from the hydro update. Boundary conditions are treated in Section 3.4.

3.1. Hydro Solve

The first step of the operator splitting is an update of the hydrodynamic equations, including the advection and compression of the radiation energy density, electron internal energy density, and the level sets. We have implemented two variants to solve the hydrodynamic equations: using (1) conservation of the total energy (Section 3.1.1) and (2) a non-conservative pressure formulation (Section 3.1.2). We can also combine the two discretizations in a hybrid manner within a multi-material simulation (Karni 1996).

3.1.1. Conservative

We have implemented several hydrodynamic shock-capturing schemes in the CRASH code: the HLLE scheme (Harten et al. 1983; Einfeldt et al. 1991), the Rusanov scheme (Yee 1989), and a Godunov scheme (Godunov 1959) with an exact Riemann solver. In this section, we will explain how we generalized the HLLE scheme for our system of equations that includes radiation, level sets, and an EOS. The other hydrodynamic schemes can be generalized in a similar fashion.

Typical hydrodynamic solvers in the literature assume constant γ. Our problem is to generalize the constant γ hydro solvers for the case of a spatially varying polytropic index, γ_e, which varies due to ionization, excitation, and Coulomb interactions. A method that is applicable to all the aforementioned, constant-γ, hydrodynamic shock-capturing schemes is one of splitting the electron internal energy E_e density into the sum of an ideal (translational) energy part p_e/(γ − 1) and an extra internal energy density E_X. Similarly, we can define an ideal total energy density

$\begin{equation} e_I = \frac{\rho u^2}{2} + \frac{p_i+p_e}{\gamma -1} + \sum _{g=1}^G E_g, \end{equation} \tag{ 28 }$

which is related to the total energy density by e = e_I + E_X. We time advance p_e with the ideal electron pressure equation and E_X by a conservative advection equation, and then apply a correction step as described below.

The time update with the operator R_hydro solves the following equations:

$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \mathbf {u}) = 0, \end{equation} \tag{ 29 }$

$\begin{equation} \frac{\partial \rho \mathbf {u} }{\partial t} + \nabla \cdot [ \rho \mathbf {u} \mathbf {u} + I(p+p_r)] = 0, \end{equation} \tag{ 30 }$

$\begin{equation} \frac{\partial e_I}{\partial t} + \nabla \cdot [ (e_I+p+p_r) \mathbf {u} ] = 0, \end{equation} \tag{ 31 }$

$\begin{equation} \frac{1}{\gamma -1}\frac{\partial p_e}{\partial t} + \frac{1}{\gamma -1}\nabla \cdot (p_e \mathbf {u})+p_e \mathbf {\nabla } \cdot \mathbf {u} = 0, \end{equation} \tag{ 32 }$

$\begin{equation} \frac{\partial E_X}{\partial t} + \nabla \cdot [ E_X \mathbf {u} ] =0, \end{equation} \tag{ 33 }$

$\begin{equation} \frac{\partial E_g}{\partial t} + \nabla \cdot (E_g \mathbf {u}) + (\gamma _r - 1)E_g \mathbf {\nabla } \cdot \mathbf {u} = 0, \end{equation} \tag{ 34 }$

$\begin{equation} \frac{\partial d_m}{\partial t} + \nabla \cdot (d_m \mathbf {u}) - d_m \mathbf {\nabla } \cdot \mathbf {u} =0, \end{equation} \tag{ 35 }$

where the frequency advection, diffusion, and energy exchange terms are omitted in this first operator step. After each time advance from time tⁿ to time tⁿ⁺¹, we have to correct e, e_I, p_e, and E_X. We denote the uncorrected variables with a superscript *, then we recover at time level n + 1 the true electron internal energy Eⁿ⁺¹_e and the true total energy density eⁿ⁺¹ by

$\begin{equation} E_e^{n+1} = \frac{p_e^*}{\gamma -1} + E_X^*, \end{equation} \tag{ 36 }$

$\begin{equation} e^{n+1} = e_I^* + E_X^*. \end{equation} \tag{ 37 }$

Since both e_I and E_X follow a conservation law, the total energy density e is also conserved. The true electron pressure is recovered from the updated electron internal energy and mass density by means of the EOS:

$\begin{equation} p_e^{n+1} = p_{\rm EOS}\big(\rho ^{n+1}, E_e^{n+1},m\big), \end{equation} \tag{ 38 }$

where the function p_EOS can be either a calculated EOS or an EOS lookup table for material m, determined by the level set functions dⁿ⁺¹_m (Section 2.3). The extra internal energy E_X is reset as the difference between the true electron internal energy and the ideal electron internal energy for γ = 5/3:

$\begin{equation} E_X^{n+1} = E_e^{n+1} - \frac{p_e^{n+1}}{\gamma -1}. \end{equation} \tag{ 39 }$

This is positive because the EOS state p_EOS satisfies E_e − p_e/(γ − 1) ⩾ 0 at all times. The ideal part of the total energy density at time level n + 1 can now be updated as

$\begin{equation} e_I^{n+1} = e^{n+1} - E_X^{n+1}. \end{equation} \tag{ 40 }$

We have now recovered eⁿ⁺¹, e_Iⁿ⁺¹, pⁿ⁺¹_e, and Eⁿ⁺¹_X at time tⁿ⁺¹.

We time advance the hydrodynamic equations to the time level * with a shock-capturing scheme with a constant γ = 5/3. For an ideal EOS, the speed of sound of Equations (29)–(34) can be derived as

$\begin{equation} c_s = \sqrt{\frac{\gamma (p_i+p_e) + \gamma _r p_r}{\rho }}, \end{equation} \tag{ 41 }$

which includes modifications due to the presence of the total radiation pressure. This speed of sound will be used in the hydro scheme below. Since the CRASH EOS solver always satisfies E_X ⩾ 0 and γ_e ⩽ 5/3, the speed of sound for the ideal EOS is always an upper bound for the true speed of sound.

We use shock-capturing schemes to advance Equations (29)– (35). In the following, we denote the (near) conservative variables by U = (ρ, ρu, e_I, p_e, E_X, E_g, d_m) and let U be grid cell averages in the standard finite volume sense. If we assume for the moment a 1D grid with spacing Δx, cell center index i, and cell face between cell i and i + 1 identified by half indices i + 1/2, then we can write the two-stage Runge–Kutta hydro update as

$\begin{equation} U_i^{n+1/2} = U_i^n - \frac{\Delta t}{2\Delta x}\left(f_{i+1/2}^n - f_{i-1/2}^n \right), \end{equation} \tag{ 42 }$

$\begin{equation} U_i^{n+1} = U_i^n - \frac{\Delta t}{\Delta x}\left(f_{i+1/2}^{n+1/2} - f_{i-1/2}^{n+1/2} \right), \end{equation} \tag{ 43 }$

where f is the numerical flux. In particular, the HLLE flux f equals the physical flux F(U^R_i+1/2) when c⁺_s = u_i + c_s ⩽ 0, F(U^L_i+1/2) when c⁻_s = u_i − c_s ⩾ 0, and in all other cases it uses the weighted flux

$\begin{equation} f_{i+1/2} = \frac{c_s^+F\big(U_{i+1/2}^L\big) - c_s^-F\big(U_{i+1/2}^R\big) + c_s^+c_s^-\big(U_{i+1/2}^R - U_{i+1/2}^L\big)}{c_s^+ - c_s^-}. \end{equation} \tag{ 44 }$

Here the left and right cell face states are

$\begin{equation} U_{i+1/2}^L = U_i + \frac{1}{2}\bar{\Delta }^LU_i, \end{equation} \tag{ 45 }$

$\begin{equation} U_{i+1/2}^R = U_{i+1} - \frac{1}{2}\bar{\Delta }^RU_{i+1}. \end{equation} \tag{ 46 }$

We use the generalized Koren limiter and define the limited slopes as

$\begin{eqnarray} &&\bar{\Delta }^LU_i = \nonumber\\ &&{\rm minmod}\left[ \beta (U_{i+1}-U_i), \beta (U_i - U_{i-1}), \frac{2U_{i+1}-U_i-U_{i-1}}{3} \right],\nonumber\\ \end{eqnarray} \tag{ 47 }$

$\begin{eqnarray} &&\bar{\Delta }^RU_i =\nonumber\\ &&{\rm minmod}\left[ \beta (U_{i+1}-U_i), \beta (U_i - U_{i-1}), \frac{U_{i+1}-U_i-2U_{i-1}}{3} \right],\nonumber\\ \end{eqnarray} \tag{ 48 }$

for the extrapolations from the left and right. This reconstruction can be third order in smooth regions away from extrema (Koren 1993; Tóth et al. 2008). The parameter β can be changed between 1 and 2, but in simulations with AMR we have best experience with β = 3/2. We generally apply the slope limiters on the primitive variables (ρ, u, p_i, p_e, E_X/ρ, E_g, d_m), instead of the conservative variables. We apply the slope limiter on E_X/ρ, instead of E_X, since E_X/ρ is smoother at shocks and across material interfaces. A multi-dimensional update is obtained by adding the fluxes for each direction in a dimensionally unsplit manner.

After each stage of the two step Runge–Kutta, we correct for the EOS effects via the update procedure outlined in Equations (36)–(40).

3.1.2. Non-conservative Pressure Equations

In regions away from shocks, it is sometimes more important to preserve pressure balance than to have a shock-capturing scheme that recovers the correct jump conditions. This is especially important at material interfaces. We therefore have implemented the option to solve the hydro part of the pressure equations

$\begin{eqnarray} \frac{\partial p_i}{\partial t} &+& \nabla \cdot (p_i \mathbf {u}) + (\gamma -1)p_i \mathbf {\nabla } \cdot \mathbf {u} = 0, \end{eqnarray} \tag{ 49 }$

$\begin{eqnarray} &&\frac{\partial p_e}{\partial t} + \nabla \cdot (p_e \mathbf {u}) + (\gamma _{Se} -1)p_e \mathbf {\nabla } \cdot \mathbf {u} = 0, \end{eqnarray} \tag{ 50 }$

instead of the equations for the total energy (31) and the electron internal energy (32). As long as the speed-of-sound gamma for the electrons

$\begin{equation} \gamma _{S_e} = \frac{\rho }{p_e} \left(\frac{\partial p_e}{\partial \rho }\right)_{S_e} \end{equation} \tag{ 51 }$

is smaller than γ = 5/3, the numerical scheme is stable. Contrary to the energy conserving scheme, the pressure-based scheme can directly include the EOS and we no longer need the time evolution of the extra internal energy density (33). The EOS contribution in the electron-pressure Equation (50) is implemented as a source term −(γ_Se − γ)p_e∇ · u added to the ideal electron pressure equation.

To facilitate using both the shock-capturing properties and the pressure balance at the material interfaces during CRASH simulations, we have several criteria to switch between them automatically. One of the criteria, for instance, uses a detection of steep pressure gradients as a shock identification. The user can select the magnitude of the pressure gradient above which the scheme switches to the conservative energy equations.

3.2. Frequency Advection

The set of multi-group Equation (14) contains an integral over the group photon frequencies. Performing this integration, the frequency advection update by the R_frequency operator can be written as

$\begin{equation} \frac{{\partial }E_g}{\partial t} -(\gamma _r-1)(\nabla \cdot {\bf u})\left[ \nu _{g+1/2} E_\nu (\nu _{g+1/2}){-} \nu _{g-1/2} E_\nu (\nu _{g-1/2})\right] = 0. \end{equation} \tag{ 52 }$

These equations, however, still depend on the, as yet, unassigned photon group frequencies ν_g and the spectral radiation energy density E_ν. We will now restrict the analysis to a frequency grid that is uniformly spaced in the frequency logarithm, i.e.,

$\begin{equation} \ln (\nu _{g+1/2}) - \ln (\nu _{g-1/2}) = \Delta (\ln \nu) = {\rm constant}. \end{equation} \tag{ 53 }$

For a large enough number of frequency groups G, the group energy E_g can then be approximated as the product of the photon frequency, spectral radiation energy E_ν, and the logarithmic group spacing Δ(ln ν):

$\begin{equation} E_g = \int _{\nu _{g-1/2}}^{\nu _{g+1/2}}E_\nu d\nu = \int _{\ln \nu _{g-1/2}}^{\ln \nu _{g+1/2}}E_\nu \nu d(\ln \nu) \approx E_\nu \nu \Delta (\ln \nu). \end{equation} \tag{ 54 }$

Using this approximation in Equation (52), we obtain our final form of the frequency advection

$\begin{equation} \frac{\partial E_g}{\partial t} + u_\nu \frac{E_{g+1/2}- E_{g-1/2}}{\Delta (\ln \nu)} = 0, \end{equation} \tag{ 55 }$

where u_ν = −(γ_r − 1)∇ · u is the frequency advection speed. The values E_g±1/2 are to be interpolated from the mesh-centered values E_g toward the group boundaries.

The frequency advection is a conservative linear advection in the log-frequency coordinate, for which the physical flux is defined as F_g−1/2 = u_νE_g−1/2. For the boundary conditions in the frequency domain we assume zero radiation flux so that no radiation can leak at the edges of the frequency domain. Equation (55) can be discretized with the one-stage second-order upwind scheme

$\begin{equation} E_g^{n+1} = E_g^* - \Delta t \frac{f_{g+1/2}-f_{g-1/2}}{\Delta (\ln \nu)}, \end{equation} \tag{ 56 }$

where time level * is now the state after the hydro update and the numerical flux is

$\begin{eqnarray} f_{g-1/2} &=& u_\nu \left[E_g-\frac{1-{}C}{2}\bar{\Delta } (E_{g+1}-E_{g},E_{g}-E_{g-1})\right]\!, \nonumber\\ u_\nu &\le& 0,\nonumber\\ f_{g-1/2} &=& u_\nu \left[E_{g-1}+\frac{1-{}C}{2}\bar{\Delta } (E_g-E_{g-1},E_{g-1}-E_{g-2})\right]\!, \nonumber\\ u_\nu &\ge& 0, \end{eqnarray} \tag{ 57 }$

and we use the superbee limiter (Roe 1986) for the limited slope $\bar{\Delta }$ . The Courant–Friedrichs–Lewy (CFL) number C = |u_ν|Δt/Δ(ln ν) depends on the hydrodynamic time step Δt. If C is larger than 1, the frequency advection is sub-cycled with the number of steps equal to the smallest integer value larger than C.

3.3. Implicit Diffusion and Energy Exchange

The stiff parts of the radio hydrodynamic equations are solved implicitly in an operator-split fashion. These stiff parts are the radiation energy diffusion, electron heat conduction, and the energy exchange between the electrons and each energy group g and between the electrons and ions. In this section, we will describe two implicit schemes that are implemented: (1) solving all radiation groups, together with electron and ion temperatures in a coupled manner (Section 3.3.1) and (2) solving each radiation group energy and the electron temperature independently (Section 3.3.2). Our strategy for resolution changes is described in Appendix A, while the modifications for the rz-geometry are explained in Appendix B.

3.3.1. Coupled Implicit Scheme

Discretizing the diffusion and energy exchange terms of Equations (19)–(20), and (14) implicitly in time leads to

$\begin{eqnarray} &&\frac{E_i^{n+1}-E_i^*}{\Delta t} = \sigma _{ie}^*\big(T_e^{n+1} - T_i^{n+1}\big), \end{eqnarray} \tag{ 58 }$

$\begin{eqnarray} \frac{E_e^{n+1}-E_e^*}{\Delta t} &=& \sigma _{ie}^*\big(T_i^{n+1} - T_e^{n+1}\big) + \nabla \cdot C_e^*\nabla T_e^{n+1} \nonumber \\ & +& \sum _{g=1}^G \sigma _g^*\big(E_g^{n+1} - B_g^{n+1}\big), \end{eqnarray} \tag{ 59 }$

$\begin{eqnarray} &&\frac{E_g^{n+1}-E_g^*}{\Delta t} = \sigma _g^*\big(B_g^{n+1} - E_g^{n+1}\big) + \nabla \cdot D_g^*\nabla E_g^{n+1}, \end{eqnarray} \tag{ 60 }$

where time level * now corresponds to the state after the hydro update and the frequency advection. The coupling coefficients σ*_ie and σ*_g and the diffusion coefficients C*_e and D*_g are taken at time level * (frozen coefficients). One can either (1) solve the coupled system of G + 2 Equations (58)–(60) implicitly or (2) solve Equation (58) for the ion internal energy Eⁿ⁺¹_i, substitute the solution back into Equation (59), and solve the resulting reduced set of G + 1 Equations (59)–(60) implicitly. Here we describe the second scheme, because it is more efficient, especially for a small number of groups, e.g., for gray radiation diffusion. Note that if we had included ion heat conduction in Equation (60), then we would have to solve the entire coupled system of equations.

First, we introduce the ion Planck function B_i = aT⁴_i as a new variable similar to the electron Planck function B = aT⁴_e, and replace E_i and E_e with these variables using the chain rule

$\begin{equation} \frac{\partial E_i}{\partial t} = \frac{\partial E_i}{\partial T_i} \frac{\partial T_i}{\partial B_i} \frac{\partial B_i}{\partial t} = \frac{C_{Vi}}{4aT_i^3}\frac{\partial B_i}{\partial t}, \qquad \frac{\partial E_e}{\partial t} = \frac{C_{Ve}}{4aT_e^3}\frac{\partial B}{\partial t}, \end{equation} \tag{ 61 }$

where C_Vi and C_Ve are the specific heats of the ions and electrons, respectively. Now Equation (58) can be replaced with

$\begin{equation} B_i^{n+1} = B_i^* + \Delta t \sigma _{ie}^{\prime } \big(B^{n+1} - B_i^{n+1}\big), \end{equation} \tag{ 62 }$

where

$\begin{equation} \sigma _{ie}^{\prime } = \sigma _{ie}^*\frac{4aT_i^3}{C_{Vi}} \frac{1}{a(T_e+T_i)\big(T_e^2+T_i^2\big)} \end{equation} \tag{ 63 }$

is again taken at time level *. The numerator comes from (T⁴_e − T⁴_i)/(T_e − T_i). Equation (62) can be solved for Bⁿ⁺¹_i. This result can be substituted into the electron internal energy Equation (59) to obtain

$\begin{eqnarray} \frac{C_{Ve}^{\prime }}{\Delta t} (B^{n+1} - B^*) &=& \sigma _{ie}^{\prime \prime }(B_i^* - B^{n+1}) + \nabla \cdot C_e^{\prime }\nabla B^{n+1} \nonumber\\ & +& \sum _{g=1}^G \sigma _g^*\big(E_g^{n+1} - w_g^* B^{n+1}\big), \end{eqnarray} \tag{ 64 }$

where we have introduced new coefficients at time level *:

$\begin{equation} \sigma _{ie}^{\prime \prime } = \frac{C_{Vi}}{4aT_i^3} \frac{\sigma _{ie}^{\prime }}{1+\Delta t \sigma _{ie}^{\prime }}, \quad C_e^{\prime } = \frac{C_e^*}{4aT_e^3}, \quad C_{Ve}^{\prime } = \frac{C_{Ve}^*}{4aT_e^3}. \end{equation} \tag{ 65 }$

The Planck weights w*_g = B*_g/B* satisfy ∑_gw_g = 1. It is convenient to introduce the changes ΔB = Bⁿ⁺¹ − B* and ΔE_g = Eⁿ⁺¹_g − E*_g to arrive at

$\begin{eqnarray} &&\left[ \frac{C_{Ve}^{\prime }}{\Delta t} + \sigma _{ie}^{\prime \prime } - \nabla \cdot C_e^{\prime }\nabla \right] \Delta B - \sum _{g=1}^G \sigma _g^* \big(\Delta E_g - w_g^* \Delta B\big) \qquad\quad \nonumber \\ &&\quad=\sigma _{ie}^{\prime \prime }\big(B_i^* - B^*\big) + \nabla \cdot C_e^{\prime }\nabla B^* + \sum _{g=1}^G\sigma _g^*\big(E_g^* - w_g^* B^*\big), \end{eqnarray} \tag{ 66 }$

$\begin{eqnarray} &&\left[ \frac{1}{\Delta t} - \nabla \cdot D_g^*\nabla \right] \Delta E_g - \sigma _g^*\big(w_g^*\Delta B - \Delta E_g\big) \nonumber\\ &&\quad= \sigma _g^*\big(w_g^*B^* - E_g^*\big) + \nabla \cdot D_g^* \nabla E_g^*. \end{eqnarray} \tag{ 67 }$

This is a coupled system of G + 1 linearized equations for the changes ΔB and ΔE_g. The right-hand sides are all at time level *.

A discrete set of equations is obtained by applying the standard finite volume method to Equations (66) and (67) and partitioning the domain in a set of control volumes V_i, enumerated by a single index i = 1, ..., I. As an example, the fluxes F_gij associated with the radiation diffusion operator may be obtained by approximating the gradient of the group energy density with a simple central difference in the uniform part of the mesh:

$\begin{equation} {-} \int _{V_i} \nabla \cdot (D_g\nabla E_g)dV = \sum _{j} F_{{\rm gij}} = \sum _{j} S_{ij}D_{{\rm gij}}\frac{E_{gi} - E_{gj}}{|{\bf x}_i - {\bf x}_j|}, \end{equation} \tag{ 68 }$

where the index j enumerates the control volumes which have a common face with the control volume i, the face area being S_ij, and the distance between cell centers is |x_i − x_j|. Note that we assumed here an orthogonal mesh. Generalization to curvilinear grids can be done as shown in Tóth et al. (2008). The diffusion coefficients at a face are obtained by simple averaging of the cell centered diffusion coefficient: D_gij = (D_gi + D_gj)/2. The discretization of the diffusion operator at resolution changes is described in Appendix A.

The linear system (66)–(67) can be written in a more compact form as the linearized implicit backward Euler scheme

$\begin{equation} \left({\bf I} - \Delta t \frac{\partial {\bf R}}{\partial {\bf U}} \right) \Delta {\bf U} = \Delta t {\bf R}({\bf U}^*), \end{equation} \tag{ 69 }$

where U are the I × (G + 1) state variables B and E_g for all I control volumes and ΔU = Uⁿ⁺¹ − U*. R is defined by the spatially discretized version of the right-hand side of Equations (66) and (67). The matrix A = I − Δt∂R/∂U is a I × I block matrix consisting of (G + 1) × (G + 1) sub-matrices. This matrix A is in general non-symmetric due to the Planck weight w*_g in the energy exchange between the radiation and electrons. To solve this system we use Krylov sub-space type iterative solvers, like GMRES (Saad & Schultz 1986) or Bi-CGSTAB (van der Vorst 1992). To accelerate the convergence of the iterative scheme, we use a preconditioner. In the current implementation of CRASH, we use the Block Incomplete Lower-Upper decomposition (BILU) preconditioner, which is applied for each adaptive mesh refinement block independently. For gray radiation diffusion the Planck weight is one, and the matrix A can be proven to be symmetric positive definite (SPD) for commonly used boundary conditions (see, for example, Edwards 1996). In that case, we can use a preconditioned conjugate gradient (PCG) scheme (see, for instance, Eisenstat 1981).

For some verification tests, we can attempt to go second order in time under the assumption of temporally constant coefficients using the Crank–Nicolson scheme

$\begin{equation} \frac{{\bf U}^{n+1} - {\bf U}^*}{\Delta t} = (1 - \alpha) {\bf R}({\bf U}^*) + \alpha {\bf R}({\bf U}^{n+1}), \end{equation} \tag{ 70 }$

with α = 1/2. The implicit residual can again be linearized R(Uⁿ⁺¹) = R(U*) + (∂R/∂U)*ΔU to obtain the linear system of equations

$\begin{equation} \left({\bf I} - \alpha \Delta t \frac{\partial {\bf R}}{\partial {\bf U}} \right) \Delta {\bf U} = \Delta t {\bf R}({\bf U}^*). \end{equation} \tag{ 71 }$

We use the same iterative solvers as for the backward Euler scheme.

Finally, we show how we use the solution ΔB and ΔE_g for g = 1, ..., G from the non-conservative Equations (66) and (67) to advance the solution of the original Equations (58)–(60) and still conserve the total energy. One needs to express the fluxes and energies on the right-hand side in the latter equations in terms of Bⁿ⁺¹ and Eⁿ⁺¹_g while still keeping the coefficients frozen. After some algebra we obtain

$\begin{eqnarray} E_i^{n+1} &=& E_i^* + \Delta t \sigma _{ie}^{\prime \prime } \big(B^{n+1} - B_i^*\big), \end{eqnarray} \tag{ 72 }$

$\begin{eqnarray} E_e^{n+1} &=& E_e^* + C_{Ve}^{\prime }(B^{n+1} - B^*), \end{eqnarray} \tag{ 73 }$

$\begin{eqnarray} E_g^{n+1} &=& E_g^* + \Delta E_g. \end{eqnarray} \tag{ 74 }$

This update conserves the total energy to round-off error. Note that at this final stage, taking too large time step may result in negative ion internal energy Eⁿ⁺¹_i if Bⁿ⁺¹ ≪ B*_i and negative electron internal energy Eⁿ⁺¹_e if Bⁿ⁺¹ ≪ B*. If this happens, the advance might be redone with a smaller time step, to limit the drop in B, or by some other time step control scheme. A generalization of the conservative update to the Crank–Nicolson scheme is also implemented for verification tests with time constant coefficients.

For completeness, we mention that in the absence of radiation we solve during the implicit step for the temperatures T_e and T_i instead of the radiation-energy-like variables aT⁴_e and aT⁴_i. In that case the corresponding matrix A is always SPD. In principle, the formulation in temperatures can be generalized to radiation as well. In Landau & Lifshitz (1980), a spectral temperature T_ν(E_ν, ν) is defined, such that the spectral energy density is locally equal to the spectral Planckian energy density at the temperature T_ν: E_ν = B_ν(T_ν, ν). This relationship is a one-to-one map. A group temperature, T_g, can also be introduced as the discrete analog such that the group energy density can be obtained by

$\begin{equation} E_g(T_g) = \int _{\nu _{g-1/2}}^{\nu _{g+1/2}} B_\nu (T_g,\nu)d\nu. \end{equation} \tag{ 75 }$

Equation (60) can be recast as equation for the group temperature T_g. This introduces the group specific heat of the radiation C_g = dE_g/dT_g. The set of Equations (58)–(60) reformulated as an implicit backward Euler scheme for the temperatures T_i, T_e, and T_g can in a similar way as in Edwards (1996) be proven to be SPD. While this scheme has the advantage of being SPD, the conservative update of the group energy density Eⁿ⁺¹_g = E*_g + C*_gΔT_g might result in negative energy density Eⁿ⁺¹_g for time steps that are too large.

3.3.2. Decoupled Implicit Scheme

The coupled implicit scheme of Section 3.3.1 requires solution of a large system of equations (G + 1 variables per mesh cell). The preconditioning of such a system can be computationally expensive and requires overall a lot of memory. We therefore also implemented a decoupled implicit scheme that solves each equation independently.

For some applications, the electron temperature does not change much in exchanging energy with the radiation. This is typically so if the electrons have a much larger energy density than the radiation, so that T_e changes little due to interaction with the radiation in a single time step. In that case, we first solve for the electron and ion temperatures without the contributions from the radiation-electron energy exchange. Let again time level * indicate the state after the hydro update and frequency advection, and again freeze C*_e, D*_g, σ*_ie, σ*_g at time level *. Discretization in time now leads to

$\begin{eqnarray} \frac{E_i^{n+1}-E_i^*}{\Delta t} &=& \sigma _{ie}^*\big(T_e^{**} - T_i^{n+1}\big), \end{eqnarray} \tag{ 76 }$

$\begin{eqnarray} \frac{E_e^{**}-E_e^*}{\Delta t} &=& \sigma _{ie}^*\big(T_i^{n+1} - T_e^{**}\big) + \nabla \cdot C_e^*\nabla T_e^{**}, \end{eqnarray} \tag{ 77 }$

where the time level ** of E_e indicates that we still have to perform an extra update to time level n + 1 with the radiation-electron energy exchange. Each radiation group energy density is solved independently using time level * for the electron temperature in B*_g:

$\begin{equation} \frac{E_g^{n+1}-E_g^*}{\Delta t} = \sigma _g^*\big(B_g^* - E_g^{n+1}\big) + \nabla \cdot D_g^*\nabla E_g^{n+1}, \end{equation} \tag{ 78 }$

where we have exploited the assumption that B*_g is not stiff.

Equations (76)–(78) can be recast in equations for the G + 1 independent changes ΔB = B** − B* and ΔE_g = Eⁿ⁺¹_g − E*_g:

$\begin{eqnarray} \left[ \frac{C_{Ve}^{\prime }}{\Delta t} + \sigma _{ie}^{\prime \prime } - \nabla \cdot C_e^{\prime }\nabla \right] \Delta B &=& \sigma _{ie}^{\prime \prime }\big(B_i^* - B^*\big) + \nabla \cdot C_e^{\prime }\nabla B^*,\nonumber\\ \end{eqnarray} \tag{ 79 }$

$\begin{eqnarray} \left[ \frac{1}{\Delta t} +\sigma _g^* - \nabla \cdot D_g^*\nabla \right] \Delta E_g &=& \sigma _g^*\big(w_g^*B^* - E_g^*\big) + \nabla \cdot D_g^* \nabla E_g^*.\nonumber\\ \end{eqnarray} \tag{ 80 }$

where we have used the definitions (61), (63), and (65) of the coefficients, frozen at time level *. Each equation for the changes is in the form of the linearized implicit backward Euler scheme (69) and can be solved independently with iterative solvers like GMRES and Bi-CGSTAB using a BILU preconditioner. As long as the boundary conditions are such that the matrices are symmetric and positive definite, a PCG method might also be used.

In a manner similar to the coupled implicit scheme, a conservative update for the energy densities can be derived as

$\begin{equation} E_g^{n+1} = E_g^* + \Delta E_g, \end{equation} \tag{ 81 }$

$\begin{equation} E_i^{n+1} = E_i^* + \Delta t \sigma _{ie}^{\prime \prime } \big(B^{**} - B_i^*\big), \end{equation} \tag{ 82 }$

$\begin{equation} E_e^{n+1} = E_e^* + C_{Ve}^{\prime }(B^{**} - B^*) + \Delta t \sum _{g=1}^G \sigma _g^* \big(E_g^{n+1}-w_g^*B^*\big), \end{equation} \tag{ 83 }$

which preserves the total energy to round-off errors. The main difference between the conservative update in the coupled and decoupled schemes is that here the energy exchange between the radiation and electrons is added afterward as the last term in Equation (83).

This scheme requires less computational time for preconditioning and for the Krylov solver than the coupled implicit algorithm. However, it generally needs more message passing in parallel computations. It is therefore not always guaranteed that the decoupled scheme is faster. The memory usage is always smaller.

3.4. Boundary Conditions

The CRASH code allows for any user specified type of boundary conditions. Several commonly used boundary conditions are readily available in the main code for convenience, e.g., fixed, extrapolation with zero gradient, periodic, and reflective.

For the radiation field, we have implemented a zero or fixed incoming flux boundary condition that is used instead of the extrapolation with zero gradient. This type of boundary condition is useful if there are no sources of radiation outside the computational domain and we assume that outflowing radiation does not return back into the computational domain (zero albedo). Note that simple extrapolation with zero gradient can make the radiation diffusion problem ill-posed. The boundary condition is derived as follows. Radiation diffusion approximation corresponds to a linear-in-angle intensity distribution

$\begin{equation} I_g = \frac{c}{4\pi }E_g + \frac{3}{4\pi } {\bf F}_g\cdot {\bf n}, \end{equation} \tag{ 84 }$

so we can calculate the radiation flux through a boundary surface. If we define the outward pointing normal vector of the boundary as n_b, the net flux of radiation energy inward through this boundary is

$\begin{equation} F^{\rm in}_g = {-}\int _{{\bf n}\cdot {\bf n}_b<0} {\bf n}_b\cdot {\bf n} I_gd\Omega = \frac{cE_g}{4} -\frac{1}{2}{\bf n}_b\cdot {\bf F}_g, \end{equation} \tag{ 85 }$

where the closure (84) is used. In the radiation diffusion model, the flux is written as F_g = −D_g∇E_g, where the diffusion coefficient D_g is a non-linear function of E_g and ∇E_g in an FLD model. The boundary condition satisfies

$\begin{equation} E_g + \frac{2D_g}{c} {\bf n}_b\cdot \nabla E_g = \frac{4}{c}F_g^{\rm in}. \end{equation} \tag{ 86 }$

For the left boundary in the x-direction, for instance, this can be discretized as

$\begin{equation} \frac{E_{g0}+E_{g1}}{2} - \frac{2D_g}{c} \frac{E_{g1} - E_{g0}}{\Delta x} = \frac{4}{c}F_g^{\rm in}, \end{equation} \tag{ 87 }$

where the index 1 corresponds to the last physical cell and 0 to the ghost cell. This equation can be solved for the ghost cell value. For zero incoming radiation flux boundary conditions we set Fⁱⁿ_g = 0.

4. CODE VERIFICATION

To test the CRASH as well as the BATS-R-US and SWMF codes, we have implemented numerous tests. These tests are subdivided in two categories: functionality tests and verification tests. Both test suites are performed automatically and return pass or fail messages depending on whether or not certain predefined tolerance criteria are met. This automated testing process provides software quality confidence especially when used in combination with a software version control system like Concurrent Versions System to recover previous correctly performing code.

The functionality tests are performed nightly on several computer platforms with different compilers and numbers of processors. They consist of unit tests and full system tests. Unit tests are designed to test a particular unit, for example, a linear equation solver. The full system tests, on the other hand, exercise the code in the way end-users will use it for their research applications. We always try to cover as much code as possible with these tests so that we can discover bugs and other unwanted side effects early on.

To test the correctness of the implemented algorithms we have also constructed a suite of verification tests. This suite is executed daily on a dedicated parallel computer and runs specific simulations to quantify against analytic and semi-analytic solutions, whenever possible. The CRASH test repository currently covers a wide range of tests for hydrodynamics, multi-material advection methods, gray and multi-group radiation diffusion, and heat conduction, to mention a few. These are performed to test for grid and/or time convergence, as deemed necessary. We also simulate full system laboratory experiment configurations in various geometries, dimensionalities, and physics fidelities. The results are either validated against laboratory experiments or simply used to check that the code keeps performing these simulations as expected. Once a week, we also perform a parallel scalability test on a large parallel computer to verify that the code does not degrade in performance during further development of the software.

In the following sub-sections, we highlight some specific verification tests related to the implicit radiation (Section 4.2) and heat conduction (Section 4.3) solver. These tests cover both Cartesian and rz-geometry, and some of them also involve the hydrodynamic solver. We demonstrate a 3D full system test in Section 4.4 and describe the parallel scalability in Section 4.5.

4.1. Error Assessment

For assessment of the accuracy of the solutions in the test suites, an appropriate definition of the numerical errors has to be defined. We use two types of errors to quantify the verification analysis. The relative L1 error is defined as

$\begin{equation} E_{{\rm L1}} = \sum _{\alpha =1}^N \frac{\sum _{i=1}^I |{\bf U}_{\alpha i} -{\bf V}_{\alpha i}|}{\sum _{i=1}^I|{\bf V}_{\alpha i}|}, \end{equation} \tag{ 88 }$

where α = 1, ..., N indexes the state variables of numerical solution vector U and the reference solution V, and i = 1, ..., I indexes the grid cells of the entire computational domain. For test problems with smooth solutions, we will also use the relative maximum error defined by

$\begin{equation} E_{{{\rm L}\infty }} = \sum _{\alpha =1}^N \frac{\max _{i=1,\ldots,I} |{\bf U}_{\alpha i}-{\bf V}_{\alpha i}|}{\max _{i=1,\ldots,I}| {\bf V}_{\alpha i}|}. \end{equation} \tag{ 89 }$

Quite often, the reference solution is defined on a grid with higher resolution than that of the numerical solution. In that case, we first coarsen the reference solution to the resolution of the numerical solution.

4.2. Radiation Tests

4.2.1. Su–Olson Test

Su & Olson (1996) developed a 1D Marshak wave test, to check the accuracy of the scheme and the correctness of the implementation of the time-dependent non-equilibrium gray radiation diffusion model. In this test, radiation propagates through a cold medium that is initially absent of radiation. The equations are linearized by the choice of the specific heat of the material C_V = 4aT³ as well as by setting the Rosseland and Planck opacities to the same uniform and time-independent constant κ_R = κ_P = κ. The cold medium is defined on a half-space of the slab geometry 0 ⩽ x < ∞. At the boundary on the left, a radiative source is specified, creating an incident radiation flux of Fⁱⁿ = aT⁴_in, where T_in = 1 keV. As time progresses, the radiation diffuses through the initially cold medium and by energy exchange between radiation and matter, the material temperature rises. In Su & Olson (1996), a semi-analytical solution is derived for the time evolution of the radiation energy and material temperature. We use this solution for our verification test.

For convenience, we locate the right boundary at a finite distance x = 5 cm and impose a zero incoming radiation flux on that boundary. We decompose the computational domain into 6 grid blocks at the base level with 10 cells per block. Between x = 5/6 cm and x = 5/3 cm, the domain is refined by one level of AMR. During the time evolution, radiation diffuses to the right through the resolution changes. The system is time-evolved with the implicit radiation diffusion solver by using a PCG method until the final time 0.02 ns. The solver steps through a series of fixed time steps of 5 × 10⁻⁴ ns and we use a Crank–Nicolson approach to achieve second-order accurate time integration. Note that this is possible because coefficients of the matrix to be solved are not time-dependent. The computed radiation and material temperatures at the final time are shown in Figure 1 and agree well with the semi-analytical solution.

**Figure 1.** Material (T_mat) and radiation (T_rad) temperature solution of the Su & Olson (1996) non-equilibrium Marshak radiation diffusion problem obtained with the CRASH code on a non-uniform grid. The reference temperatures of the analytical method of Su & Olson (1996) are shown as lines.
Download figure:
Standard image High-resolution image

Figure 2 shows the relative L1 error of the radiation and material temperatures versus increasing grid resolution of the base level grid. We did not use the semi-analytical solution as the reference, since it is difficult to get an accurate enough solution with the quadrature method as mentioned by Su & Olson (1996). Instead, we use a very high resolution (1920 cells) numerical reference solution obtained with the CRASH code. Four different base level resolutions with 60, 120, 240, and 480 cells are used to demonstrate the second-order convergence. The time step is proportional to the cell size Δx.

4.2.2. Lowrie's Non-equilibrium Radiation Hydrodynamic Solutions

Lowrie & Edwards (2008) designed several shock tube problems for the non-equilibrium gray radiation diffusion coupled to the hydrodynamic equations that can be used for code verification. These solutions are planar radiative shock waves where the material and radiation temperatures are out of equilibrium near the shock, but are assumed to be in radiative equilibrium far from the shock. Depending on the Mach number of the pre-shock state, a wide range of shock behavior can occur. For the CRASH test suite, we selected a few of the semi-analytic solutions from Lowrie & Edwards (2008). In this section we will describe the Mach 1.05 flow with uniform opacities as an example. Here the shock is smoothed by energy exchange with the diffusive radiation. Another more challenging Mach 5 problem with non-uniform opacities will be described in Section 4.3.3.

The Mach 1.05 test is performed on a 2D non-uniform grid. The initial condition is taken to be the same as the original steady state reference solution. Since the system of equations is Galilean invariant, we can add an additional velocity −1.05 so that the velocity on the left boundary is zero while the smoothed shock will now move to the left. This new initial condition as well as the velocity vector are rotated by tan⁻¹(1/2) ≈ 26 fdg 56. This means that there is a translational symmetry in the (−1, 2) direction of the xy-plane as shown in Figure 3. The computational domain is −0.12 < x < 0.12 by −0.02 < y < 0.02 decomposed in 3 × 3 grid blocks of 24 × 4 cells each. We apply one level of refinement inside the region −0.04 < x < 0.04 by −0.02/3 < y < 0.02/3. The initial smoothed shock starts at the right boundary of the refined grid and we time-evolve the solution until it reaches the resolution change on the left as shown in Figure 3. For the boundary conditions in the x-direction, we use zero radiation influx conditions for the radiation field, while a zero gradient is applied to the remaining state variables. On the y boundaries, we apply a sheared zero gradient in the (−1, 2) direction for all variables.

**Figure 3.** Rotated shock-tube test on a 2D AMR grid based on the Mach number 1.05 non-equilibrium gray radiation hydrodynamic test in Lowrie & Edwards (2008). Shown is the radiation temperature in color contour at the initial (top panel) and final (bottom panel) times. Black crosses indicate cell centers.
Download figure:
Standard image High-resolution image

The hydrodynamic equations are time-evolved with the HLLE scheme with a CFL number 0.8. We use the generalized Koren limiter with β = 3/2 for the slope reconstruction. For the implicit radiation diffusion solver, we use the GMRES iterative solver in combination with a BILU preconditioner. The specific heat is time-dependent since it depends on the density, therefore the implicit scheme is only first-order accurate in time. To enable second-order grid convergence for this smooth test problem, we compensate for this by reducing the CFL number proportional to the grid cell size, in other words Δt ∝ Δx², so that second-order accuracy with respect to Δx can be achieved. We increase the spatial resolution by each time doubling the number of grid blocks at the base level in both the x- and y-directions.

The convergence of the numerically obtained material and radiation temperatures along the y = 0 cut at the final time t = 0.07 is shown in Figure 4. The solid, dotted, and dashed lines correspond to the solutions with the 3 × 3, 6 × 6, and 12 × 12 base level grid blocks, respectively. The advected semi-analytical reference solution is shown as a blue line for comparison.

To assess the order of accuracy, the grid convergence is shown in Figure 5 for the three resolutions. The relative L1 error is calculated using the density, velocity components, and both the material and radiation temperatures. We obtain second-order convergence for both the conservative and the non-conservative (using the pressure equation instead of the total energy) hydrodynamic schemes. The latter scheme can be used because in the Mach 1.05 test the hydro shock is smoothed out by the interaction with the radiation.

4.2.3. Double Light Front

As a test for the multi-group radiation diffusion model we developed a double light front test problem. This test is used to verify the implementation of both the group diffusion and flux limiters. At the light front, the discontinuity in the radiation field switches on the flux limiter. This limiter is used to correct the radiation propagation speed in the optically thin free-streaming regime. With the light front test we can then check that we obtain the speed of light propagation of the front and that the front maintains as much as possible the initial discontinuity.

This test is constructed as follows. We use a 1D computational domain of the size of 1 m in the x-direction. On this domain, we initialize the two radiation group energy densities E_g (g = 1, 2) with a very small, positive number to avoid division by zero in the flux limiter. Also the Rosseland mean opacities are set to a small number corresponding to strong radiation diffusion, while the Planck mean opacities are set to zero corresponding to an optically thin medium. The radiation energy density of the first group enters from the left boundary by applying a fixed boundary condition with value one in arbitrary units. On the right boundary this group is extrapolated with zero gradient. Note that these are the proper boundary conditions in the free-streaming limit and not the diffusive flux boundary conditions described in Section 3.4. The second radiation group enters from the right boundary with density 1, and it is extrapolated with zero gradient at the left boundary. We time-evolve both groups for 0.5 m/c seconds. The analytic solution is then two discontinuities that have reached x = 0.5 m, since both fronts propagate with the speed of light c.

The computational domain is non-uniform. In the coarsest resolution, there are 10 grid blocks of 4 cells each at the base level. Inside the regions 0.1 < x < 0.2 and 0.8 < x < 0.9, we use one level of refinement. The total time evolution is divided into 400 fixed time steps. We use GMRES for the radiation diffusion solver in combination with a BILU preconditioner. For the grid convergence, we reduce the fixed time step quadratically with grid resolution. This time step reduction mimics second-order discretization in time. In Figure 6, the two group energy densities are shown for the base level grid resolutions 40, 80, 160, and 320. Clearly, with increasing number of cells, the solution converges toward the reference discontinuous fronts at x = 0.5.

**Figure 6.** Solutions for the 1D double light front test for four different non-uniform grid resolutions. The radiation energy for group 1 (left panel) enters from the left boundary, for group 2 (right panel) it enters from the right boundary. The symbols for base resolution 80 show one level of grid refinement for 0.1 < x < 0.2 and 0.8 < x < 0.9.
Download figure:
Standard image High-resolution image

In Figure 7, the grid convergence is shown for the four resolutions. The relative L1 error is calculated using both radiation group energy densities and compared to the analytical reference solution with the discontinuities at x = 0.5. In Gittings et al. (2008), it was stated that for a second-order difference scheme the convergence rate for a contact discontinuity is 2/3. Indeed, we find this type of convergence rate, due to numerical diffusion of the discontinuities, for the light front test. We have also performed the tests in the y- and z-directions to further verify the implementation.

4.2.4. Relaxation of Radiation Energy Test

This test is designed to check the relaxation rate between material and radiation. The energy exchange between the material and radiation groups can be written as

$\begin{equation} C_V \frac{\partial T}{\partial t} = \sum _{g=1}^G \sigma _g(E_g - B_g), \end{equation} \tag{ 90 }$

$\begin{equation} \frac{\partial E_g}{\partial t} = \sigma _g(B_g - E_g). \end{equation} \tag{ 91 }$

For a single radiation group, an analytic expression can be found to describe the relaxation in time. However, for an arbitrary number of groups, a time-dependent analytic solution is less obvious, except for some rather artificial cases. Here we assume an extremely large value of the specific heat C_V to make analysis more tractable. In this case, the material temperature is time-independent, so that B_g is likewise time independent. The solution is then $E_g = B_g(1-e^{-\sigma _g t})$ assuming E_g(t = 0) = 0 initially. At time t = 1/σ_g, the group radiation energy density is E_g = B_g(1 − 1/e). Note that this test only needs one computational mesh cell in the spatial domain. We set T = 1 keV and the resulting Planckian spectrum, defined by B_g, is depicted by the dotted line in the left panel of Figure 8. We use 80 groups logarithmically distributed over the photon energy domain in the range of 0.1 eV to 20 keV. The computed E_g at time t = 1/σ_g are shown as + points. For the simulation we used the GMRES iterative solver and the Crank–Nicolson scheme. To assess the error, we repeated the test with time steps of 1/20, 1/40, and 1/80 of the simulation time. The second-order convergence rate is demonstrated in the right panel of Figure 8.

**Figure 8.** Relaxation of radiation energy test for 80 groups. Left panel is for the time-independent spectrum *B_g* (dotted line) and the group radiation energy solution *E_g* at time 1/σ_g (+ points) vs. the photon energies after 80 time steps. The analytical reference solution is shown as a solid line. Right panel shows the relative maximum error for 20, 40, and 80 time steps demonstrating second-order convergence rate.
Download figure:
Standard image High-resolution image

4.3. Heat Conduction Tests

4.3.1. Uniform Heat Conduction in rz-geometry

This test is designed to verify the implicit heat conduction solver in rz-geometry. It tests the time evolution of the temperature profile using uniform and time-independent heat conductivity. In rz-geometry, the equation of the electron temperature for purely heat conductive plasma follows

$\begin{equation} C_{Ve}\frac{\partial T_e}{\partial t} = \frac{1}{r} \frac{\partial }{\partial r}\left(rC_e \frac{\partial T_e}{\partial r} \right) + \frac{\partial }{\partial z}\left(C_e \frac{\partial T_e}{\partial z} \right)\!. \end{equation} \tag{ 92 }$

We set the electron specific heat C_Ve = 1 and assume electron conductivity C_e to be constant. In this case, a solution can be written as a product of a Gaussian profile in the z-direction and an elevated Bessel function J₀ in the r-direction (Arfken 1985):

$\begin{equation} T_e = T_{\rm min} + T_0 \frac{1}{\sqrt{4\pi C_e t}} e^{-\frac{z^2}{4C_e t}} J_0(b r)e^{-b^2 C_e t}, \end{equation} \tag{ 93 }$

where b ≈ 3.8317 is the first root of the derivative of J₀(r). We select the following values for the input parameters: T_min = 3, T₀ = 10, and C_e = 0.1 in dimensionless units.

The computational domain is −3 < z < 3 and 0 < r < 1 discretized with 3 × 3 grid blocks of 30 × 30 cells each. In the region −1 < z < 1 and 1/3 < r < 2/3, we apply one level of mesh refinement. We impose a symmetry condition for the electron temperature on the axis. On all other boundaries the electron temperature is fixed to the time-dependent reference solution. We time-evolve this heat conduction problem with a PCG method from time t = 1 to the final time at t = 1.5. The Crank–Nicolson approach is used to achieve second-order accurate time integration.

The initial and final solutions for the electron temperature are shown in Figure 9 in color contour in the rz-plane. The heat conduction has diffused the temperature in time to a more uniform state. The black line indicates the region in which the mesh refinement was applied. The relative maximum error of the numerically obtained electron temperature versus the analytical solution is shown in Figure 10. Here we used the non-uniform grid with base resolutions of 90², 180², 360², and 720² cells and set the time step proportional to the cell size to demonstrate a second-order convergence rate.

**Figure 10.** Relative maximum error for the uniform heat conduction test on a non-uniform grid in the rz-geometry.
Download figure:
Standard image High-resolution image

4.3.2. Reinicke–Meyer-ter Vehn Test

The Reinicke & Meyer-ter-Vehn (1991) problem tests both the hydrodynamic and heat conduction implementations. This test generalizes the well-known Sedov–Taylor strong point explosion in single-temperature hydrodynamics by including the heat conduction. The heat conductivity is parameterized as a non-linear function of the density and material temperature: C_e = ρ^aT^b. We select the spherically symmetric self-similar solution of Reinicke & Meyer-ter-Vehn (1991) with coefficients a = −2 and b = 13/2 and the adiabatic index is γ = 5/4. This solution produces, similar to the Sedov–Taylor blast-wave, an expanding shock front through an ambient medium. However, at very high temperatures, thermal heat conduction dominates the fluid flows, so that a thermal front precedes the shock front. With the selected parameters, the heat front is always at twice the distance from the origin of the explosion as is the shock front.

We perform the test in rz-geometry. The computational domain is divided in 200 × 200 cells. The boundary conditions along the r and z axes are reflective. The two other boundaries, away from explosion, are prescribed by the self-similar solution. The time evolution is numerically performed as follows. For the hydrodynamics, we use the HLLE scheme with the CFL number set to 0.8. Since this test is performed on a uniform grid without AMR, we can use the generalized Koren limiter with β = 2. This is the same as the original Koren (1993) slope limiter. The heat conduction is solved implicitly with the PCG method. The test is initialized with the spherical self-similar solution with the shock front located at the spherical radius 0.225 and the heat front is at 0.45. The simulation is stopped once the shock front has reached 0.45 and the heat front is at 0.9.

A 1D slice along the r-axis of the solution at the final time is shown in Figure 11. We normalize the output similar to Reinicke & Meyer-ter-Vehn (1991). The temperature is normalized by the central temperature, while the density and radial velocity are normalized by their values of the post-shock state at the shock. The numerical solution obtained by the CRASH code is shown as + symbols and is close to the self-similar reference solution, shown as solid lines. Note that the temperature is smooth due to the heat conduction, except for the discontinuous derivative at the heat front. The wiggle at r = 0.3 in the density and radial velocity is due to the diffusion of the analytical shock discontinuity in the initial condition during the first few time steps. In the left panel of Figure 12, the spherical expansion of temperature at the final time is shown. Clearly, the Cartesian grid with the rz-geometry does not significantly distort the spherical symmetry of the solution. The spatial distribution of the error in the temperature is shown in the right panel. The errors are largest at the discontinuities of the shock and heat fronts as expected.

**Figure 12.** Temperature (left panel) and temperature error compared to the reference solution (right panel) for the Reinicke–Meyer-ter Vehn test in rz-geometry.
Download figure:
Standard image High-resolution image

A grid convergence study is performed with resolutions of 200², 400², and 800² cells. The relative L1 error in Figure 13 is calculated using the density, velocity components, and the material temperature. The convergence rate is first-order due to the shock and heat front.

4.3.3. Heat Conduction Version of Lowrie's Test

Any of the verification tests for non-equilibrium gray diffusion coupled to the single temperature hydrodynamics can be reworked as a test for the hydrodynamic equations for the ions coupled to the electron pressure equation with electron heat conduction and energy exchange between the electrons and ions. As an example, we transform one of the non-equilibrium gray diffusion tests of Lowrie & Edwards (2008) to verify the heat conduction implementation.

The electron energy density Equation (20) without the radiation interaction can be written as

$\begin{equation} \frac{\partial E_e}{\partial t} + \nabla \cdot [ E_e \mathbf {u} ] + p_e \mathbf {\nabla } \cdot \mathbf {u} = \nabla \cdot [ C_e\nabla T_e ] + \sigma _{ie} (T_i - T_e), \end{equation} \tag{ 94 }$

where the heat conduction and energy exchange terms on the right-hand side depend on the gradients and differences of the temperatures. The equation for the gray radiation energy density (16) on the other hand depends on the gradients and differences of energy densities. By defining the radiation temperature T_r by E_r = aT⁴_r and using the definition of the Planckian B = aT⁴_e, we can rewrite the energy density equation for the radiation as

$\begin{eqnarray} \frac{\partial E_r}{\partial t} + \nabla \cdot [ E_r \mathbf {u} ] + \frac{1}{3} E_r \mathbf {\nabla } \cdot \mathbf {u} &=& \nabla \cdot [ D_r \nabla E_r ]\! + c\kappa _P (aT^4 - E_r\!) \nonumber \\ &=& \nabla \cdot [ \overline{D}_r \nabla T_r ] + c\overline{\kappa }_P (T-T_r),\nonumber \\ \end{eqnarray} \tag{ 95 }$

where $\overline{D}_r = D_r 4aT_r^3$ and $c\overline{\kappa }_P = c\kappa _Pa(T^2+T_r^2)(T+T_r)$ are the new coefficients that appear due to this transformation. Equations (94) and (95) are now of the same form. To translate a gray diffusion test to a heat conduction test, we reinterpret $\overline{D}_r$ as the heat conductivity C_e and $c\overline{\kappa }_P$ as the relaxation coefficient σ_ie in the ion–electron energy exchange. In addition, the material temperature T and radiation temperature T_r have to be reinterpreted as the ion temperature T_i and electron temperature T_e, respectively. Note that we also have to relate the electron pressure and internal energy by p_e = E_e/3 similar to the radiation field corresponding to γ_e = 4/3, and let the electron internal energy and electron specific heat depend on the electron temperature as E_e = aT⁴_e and C_Ve = 4aT³_e, respectively.

As an example, we transform the Mach 5 non-equilibrium gray diffusion shock tube problem of Lowrie & Edwards (2008). It uses non-uniform opacities that depend on the density and temperature defined by D_r = 0.0175(γT)^7/2/ρ and cκ_P = 10⁶/D_r. The above described procedure is used to translate this problem to an electron heat conduction test with energy exchange between the electron and ions. The heat conductivity for this test is C_e = 4aT³_e0.0175(γT_i)^7/2/ρ and the relaxation coefficient between the electron and ions is σ_ie = a(T²_i + T²_e)(T_i + T_e)4aT³_e10⁶/C_e.

We perform this Mach 5 heat conduction test on a 2D non-uniform grid. For the initial condition, the 1D semi-analytical steady state reference solution of Lowrie & Edwards (2008) is used. There is a Mach 5 pre-shock flow on the left side of the tube resulting in an embedded hydro shock as well as a steep thermal front (a look at Figure 14 will help to understand this shock tube problem). We add an additional velocity of Mach −5 so that the pre-shock velocity is zero and the shock is no longer steady, but instead will move to the left with a velocity −5 (in units in which the pre-shock speed of sound is 1). The problem is rotated counterclockwise on the grid by tan⁻¹(1/2). The translational symmetry is now in the (−1, 2) direction in the xy-plane similar to the Mach 1.05 shock tube problem described in Section 4.2.2. The computational domain is −0.0384 < x < 0.0384 by −0.0048 < y < 0.0048. Inside the area −0.0128 < x < 0.0128 and −0.0016 < y < 0.0016, we apply one level of refinement. This refinement is set up such that both the heat front as well as the shock front propagate through the resolution change on the left (from fine to coarse) and right (from coarse to fine), respectively. For the boundary conditions in the x-direction, we fix the state on the right side with the semi-analytical solution, but for the left side we use zero gradient. On the y boundaries, we apply a sheared zero gradient in the (−1, 2) direction.

For the evolution until the final time t = 0.0025, we use the HLLE scheme together with the generalized Koren limiter with β = 3/2 to solve the hydrodynamic equations. The CFL number is set to 0.8. The heat conduction and energy exchange between electrons and ions are solved implicitly with the backward Euler scheme using the GMRES iterative solver in combination with a BILU preconditioner.

In Figure 14, the electron (right panel) and ion (left panel) temperatures are shown at the final time along the x-axis. The semi-analytical reference solution is shown as a blue line, while the numerical solution is shown with + symbols for a simulation with 192 × 24 cells at the base level in the x- and y-directions. The hydro shock is located near x ≈ 0.0085 and shows up in the ion temperature as a jump in the temperature, followed directly behind the shock by a strong relaxation due to the energy exchange between the ions and electrons. The electron temperature stays smooth due to strong heat conduction. The heat front is seen with a steep foot at x ≈ −0.022. This front corresponds to the radiative precursor in the non-equilibirum gray diffusion tests of Lowrie & Edwards (2008). We repeat the test with four different resolutions at the base level: 192 × 24, 384 × 48, 768 × 96, and 1536 × 192 cells in the x- and y-directions. The insets in both panels of Figure 14 show the four resolutions as solid, dotted, dashed, dash-dotted lines, respectively. In the left panel, the zoom-in shows the convergence of the ion temperature toward the embedded hydro shock and the temperature relaxation. In the right panel, the blow-up shows the convergence toward the reference precursor front. Note that no spurious oscillations appear near the shock or near the precursor.

Due to the discontinuity in both the shock and heat precursor, the convergence rate can be at most first-order. Indeed, in Figure 15 the relative L1 error shows first-order accuracy. The error is calculated using all the density, velocity components, and both temperatures. Note that the spike in the ion temperature is spatially so small that a huge number of grid cells are needed to get a fully resolved shock and relaxation state.

4.4. Full System Tests

The CRASH test repository contains a range of full system configurations to be used for validation with future laboratory experiments. In Figure 16, we show the configuration of a 3D elliptical nozzle through which a fast shock of the order of 150 km s⁻¹ will be launched, which is still significantly slower than the speed of light. The shock wave is produced by a 1.1 ns laser pulse from the left with 4 kJ of energy irradiating a 20 μm thick beryllium disk, initially located at x = 0. A layer of gold is glued to the plastic tube to protect the outside of the tube from the laser-driven shock. The plastic (polyimide) tube is circular for x < 500 μm with a radius of 600 μm. Beyond x = 750 μm the tube is made elliptical by flattening the tube in the y-direction by a factor two.

**Figure 16.** Geometry of the 3D elliptic nozzle experiment after 1.1 ns, consisting of 5 materials: beryllium (blue), xenon (black), polyimide (green), gold (yellow), and acrylic (red) in both panels. The radius of the inside of the polyimide tube is 600 μm in the y = 0 plane (left panel). In the z = 0 plane (right panel), the radius of the inner tube is 600 μm for x < 500 μm, but shrinks to 300 μm beyond x = 750 μm. The lines represent the mesh refinement at material interfaces and shock fronts.
Download figure:
Standard image High-resolution image

A laser energy deposition library is currently under construction and the implemented will be reported elsewhere. For this paper, we will instead perform the first part of the simulation with the 2D, Lagrangian, radiation hydrodynamics code HYADES (Larsen & Lane 1994) to time advance the heating due to the irradiation by the laser beams and the response of the plasma until 1.1 ns. This laser pulse first shocks and then accelerates the beryllium to the right. After 1.1 ns, the output of HYADES is used as an initial condition of the CRASH code.

This simulation is performed for a two-temperature, electron and ion, plasma. For the radiation, we use the FLD approximation with 30 groups. The photon energy is in the range of 0.1 eV to 20 keV, logarithmically distributed over the groups. Due to the symmetry in the problem we only simulate one quadrant (y>0 and z>0), with reflective boundary conditions at y = 0 and z = 0. At all other boundaries we use an extrapolation with zero gradient for the plasma and a zero incoming flux boundary for the radiation. The domain size is [ − 150, 3900] × [0, 900] × [0, 900] microns for the x, y, z coordinates. The base level grid consists of 120 × 20 × 20 blocks of 4 × 4 × 4 mesh cells. One level of dynamic mesh refinement is used at material interfaces and the shock front. Overall, the effective resolution is 960 × 160 × 160 cells and there are approximately 4.5 million finite volume cells. The hydrodynamic equations are solved with the HLLE scheme with a CFL number 0.8 together with the generalized Koren limiter with β = 3/2. The diffusion and energy exchange of the radiation groups as well as the heat conduction are solved with the decoupled implicit scheme using a Bi-CGSTAB iterative solver. The simulation from 1.1 ns to 13 ns physical time took 1 hr and 55 minutes on 480 cores of the FLUX supercomputer at the University of Michigan.

In Figure 17, we show the shock structure at 13 ns. The accelerated beryllium compresses the xenon directly to the right of the interface, which is seen as a high density plasma near x = 1700 μm in the top right panel of Figure 17. This drives a primary shock and the velocity jump at x ≈ 1700 μm is seen in the middle left panel. Behind the shock front, the ions are heated as depicted from the middle right panel, followed directly behind the shock by cooling due to the energy exchange between ions and electrons. Early on, electron heating produces ionization and emission of radiation, and the radiation in turn heats and ionizes the material ahead of the primary shock. The radiation temperature, as measured by the total radiation energy density, is shown in the bottom left panel. The photons will interact again with the matter, sometimes after traveling some distance. This is the source of the wall shock seen ahead of the primary shock (Doss et al. 2009, 2011): photons traveling ahead of the shock interact with the plastic wall, heat it, and in turn drive a shock off the wall into the xenon. The ablation of plastic is depicted in the top left panel as a radially inward moving polyimide (in green color) near and even ahead of the primary shock. The compressed xenon due to plastic ablation is seen in the top right panel as a faint density feature that is ahead of the primary shock front, between x = 1700 μm and x = 2000 μm. Interaction between the photons and matter is also seen by the radiative precursor to the right of the radiative shock elevating the electron temperature ahead of the shock in the bottom right panel. This is due to the strong coupling between the electrons and radiation field. The reader is referred to Drake et al. (2011) for more details on radiative effects in radiative shock tubes.

4.5. Parallel Performance

We present parallel scaling studies on the Pleiades supercomputer at NASA Ames. This computer is an SGI ICE cluster connected with InfiniBand. Figure 18 shows strong scaling for a problem size that is independent of the number of processors. This 3D simulation is a circular tube version of the full system test described in Section 4.4. It uses five materials, 30 radiation groups, and separate electron and ion temperatures. The grid contains 80 × 8 × 8 blocks of 4 × 4 × 4 cells each at the base level and in addition two time-dependent refinement levels. There are overall approximately 2.6 million cells in this problem. We use lookup tables for the EOS and opacities, so that the computational time for that is negligible. For the hydrodynamic equations, we use the HLLE scheme together with the generalized Koren limiter with β = 3/2. The radiation diffusion, electron heat conduction, and energy exchange terms are solved implicitly with the decoupled scheme, using the Bi-CGSTAB iterative solver. This simulation is performed for 20 time steps for the number of cores varying from 128 to 2048, but excludes file I/O to measure the performance of the implicit solver. Up to 1024 cores, we get good scaling. However, for more cores we observe saturation in the performance.

**Figure 18.** Strong scaling of the CRASH code, running a 3D CRASH application with five material level sets, electron and ion temperature, 30 radiation groups, and two levels of time-dependent mesh refinement.
Download figure:
Standard image High-resolution image

5. SUMMARY

We have extended the BATS-R-US code (Powell et al. 1999; Tóth et al. 2011) with a new radiation transfer and heat conduction library. This new combination, together with the EOS and multi-group opacity solver, is called the CRASH code. This code uses the recently developed parallel BATL (Tóth et al. 2011) to enable highly resolved radiation hydrodynamic solutions. The implemented radiation hydrodynamic schemes solve for gray or multi-group radiation diffusion models in the FLD approximation.

In high-energy-density plasmas, the electrons are most of the time strongly coupled to the ions by collisions. An important exception is at hydrodynamic shocks, where the ions are heated by the shock wave and the electrons and ions are out of temperature equilibrium. Since radiative shocks are the main application for CRASH, we have implemented a separate electron pressure equation with the electron thermal heat conduction. For the electron heat conduction, we have added the option of a flux limiter to limit the thermal flux with the free-streaming heat flux.

The multi-material radiation hydrodynamic equations are solved with an operator-split method that consists of three substeps: (1) solving the hydrodynamic equations with standard finite volume shock-capturing schemes, (2) the linear advection of the radiation in frequency-logarithm space, and (3) the implicit solution of the radiation, heat conduction, and energy exchanges. For the implicit solver, standard Krylov solvers are used together with a BILU preconditioner. This preconditioner scales well up to 1000 processors on present-day parallel computers. For future work, we may explore for the implicit multi-group diffusion a multi-level preconditioner to better scale the radiation solver beyond 1000 processors.

We have presented a suite of verification tests that benchmark the performance. These tests verify the correctness and accuracy of our implementation of the gray and multi-group radiation diffusion algorithm and the heat conduction in 1D, 2D, and 3D slab and 2D rz geometry. To demonstrate the full capability of the implementation, we have presented a 3D multi-material simulation of a radiative shock wave propagating through an elliptical nozzle. This configuration will be used in future validation studies.

We caution the reader that another CRASH code has already been discussed in the literature (Ciardi et al. 2001; Maselli et al. 2003, 2009; Pierleoni et al. 2009). This code is used for cosmological reionization around the first stars using Monte Carlo radiative transfer methods. Our code is different and uses multi-group radiation diffusion models in the FLD approximation. The applications are for astrophysics and laboratory astrophysics.

The authors express their gratitude to M. J. Grosskopf and E. Rutter for providing HYADES output. Discussions with B. Torralva and C. C. Chou are gratefully acknowledged. This work was funded by the Predictive Sciences Academic Alliances Program in DOE/NNSA-ASC via grant DEFC52-08NA28616 and by the University of Michigan. The simulations were performed on the NASA Advanced Supercomputing system Pleiades.

APPENDIX A: DISCRETIZATION OF THE DIFFUSION OPERATOR AT RESOLUTION CHANGES

In Sections 3.3.1 and 3.3.2, the diffusion operator is discretized on a uniform mesh with a standard finite volume method in combination with a central difference approximation for the gradient in the flux calculation as in Equation (68). The diffusion coefficient needed at the face is obtained by simple arithmetic averaging of the left and right cell center diffusion coefficients. The generalization to resolution changes as in Figure 19 is less straightforward. In the following, we denote the fine cell centers by a and b, the coarse cell center by c. The flux densities at resolution changes in the direction orthogonal to the interface are denoted by F₁ and F₂ at the fine faces, and F₃ at the coarse face.

**Figure 19.** Cell and face centers at the adaptive interface in 2D.
Download figure:
Standard image High-resolution image

In Edwards (1996), a strategy was developed to discretize the diffusion operator on an adaptive mesh in the context of reservoir simulations. The main ingredients of the method are (1) require the continuity of the flux at the resolution change in the strong sense, i.e., F₁ = F₂ = F₃ and (2) discretize the gradient in the diffusion flux by a one-sided difference. An expression was found for the diffusion flux F = −D∇E in which the diffusion coefficient is replaced by a weighted harmonic average of the cell centered values D_a, D_b, D_c. In Gittings et al. (2008), it was argued that this discretization does not properly propagate the self-similar Marshak waves of the radiation diffusion model, unless the cell centered diffusion coefficients are calculated on a common face temperature.

In the code discussed in this paper, we follow a different approach that replaces the harmonic average of the diffusion coefficient in Edwards (1996) with an arithmetic average and for the flux densities normal to the resolution change interface we obtain

$\begin{equation} F_1 = F_2 = F_3 = -\frac{2D}{3\Delta x}\left[ E_c - (E_a + E_b)/2 \right], \end{equation} \tag{ A1 }$

where Δx is the fine cell size and the diffusion coefficient at the face is averaged as

$\begin{equation} D_3 = [ D_c + (D_a+D_b) ] /3. \end{equation} \tag{ A2 }$

We have demonstrated with verification tests including those discussed above that this change produces properly propagating radiative precursor and shock fronts. Generalizations to 1D and 3D are straightforward.

APPENDIX B: rz-GEOMETRY

Incorporating rz-geometry in a finite volume formulation is as follows: the radial cell face area and the cell volume must be made proportional to the distance r from the symmetry axis. In addition, the r component of the momentum Equation (18) is modified as

$\begin{equation} \frac{\partial \rho u_r}{\partial t} + \nabla \cdot [ \rho \mathbf {u} u_r + {\bf \hat{r}}(p+p_r)] = \frac{p+p_r}{r}, \end{equation} \tag{ B1 }$

where ${\bf \hat{r}}$ is the unit vector in the r-direction and $u_r= \mathbf {u} \cdot {\bf \hat{r}}$ . This correction reflects that the pressure term is a gradient, not a divergence.

CRASH: A BLOCK-ADAPTIVE-MESH CODE FOR RADIATIVE SHOCK HYDRODYNAMICS—IMPLEMENTATION AND VERIFICATION

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. EQUATIONS OF RADIATION HYDRODYNAMICS IN DENSE PLASMAS

2.1. Radiation Transport

2.2. Hydrodynamics

2.3. Level Sets and Material Identification

2.4. Equation of State and Opacities

2.5. Flux-limited Diffusion

3. THE NUMERICAL METHOD

3.1. Hydro Solve

3.1.1. Conservative

3.1.2. Non-conservative Pressure Equations

3.2. Frequency Advection

3.3. Implicit Diffusion and Energy Exchange

3.3.1. Coupled Implicit Scheme

3.3.2. Decoupled Implicit Scheme

3.4. Boundary Conditions

4. CODE VERIFICATION

4.1. Error Assessment

4.2. Radiation Tests

4.2.1. Su–Olson Test

4.2.2. Lowrie's Non-equilibrium Radiation Hydrodynamic Solutions

4.2.3. Double Light Front

4.2.4. Relaxation of Radiation Energy Test

4.3. Heat Conduction Tests

4.3.1. Uniform Heat Conduction in rz-geometry

4.3.2. Reinicke–Meyer-ter Vehn Test

4.3.3. Heat Conduction Version of Lowrie's Test

4.4. Full System Tests

4.5. Parallel Performance

5. SUMMARY

APPENDIX A: DISCRETIZATION OF THE DIFFUSION OPERATOR AT RESOLUTION CHANGES

APPENDIX B: rz-GEOMETRY

CRASH: A BLOCK-ADAPTIVE-MESH CODE FOR RADIATIVE SHOCK HYDRODYNAMICS—IMPLEMENTATION AND VERIFICATION

Article metrics

Permissions

Share this article

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. EQUATIONS OF RADIATION HYDRODYNAMICS IN DENSE PLASMAS

2.1. Radiation Transport

2.2. Hydrodynamics

2.3. Level Sets and Material Identification

2.4. Equation of State and Opacities

2.5. Flux-limited Diffusion

3. THE NUMERICAL METHOD

3.1. Hydro Solve

3.1.1. Conservative

3.1.2. Non-conservative Pressure Equations

3.2. Frequency Advection

3.3. Implicit Diffusion and Energy Exchange

3.3.1. Coupled Implicit Scheme

3.3.2. Decoupled Implicit Scheme

3.4. Boundary Conditions

4. CODE VERIFICATION

4.1. Error Assessment

4.2. Radiation Tests

4.2.1. Su–Olson Test

4.2.2. Lowrie's Non-equilibrium Radiation Hydrodynamic Solutions

4.2.3. Double Light Front

4.2.4. Relaxation of Radiation Energy Test

4.3. Heat Conduction Tests

4.3.1. Uniform Heat Conduction in rz-geometry

4.3.2. Reinicke–Meyer-ter Vehn Test

4.3.3. Heat Conduction Version of Lowrie's Test

4.4. Full System Tests

4.5. Parallel Performance

5. SUMMARY

APPENDIX A: DISCRETIZATION OF THE DIFFUSION OPERATOR AT RESOLUTION CHANGES

APPENDIX B: rz-GEOMETRY