A&A 464, 429-435 (2007)
DOI: 10.1051/0004-6361:20065486

HERACLES: a three-dimensional radiation hydrodynamics code

M. González^1,2 - E. Audit^1,2 - P. Huynh³

1 - CEA-Saclay, DSM/DAPNIA/Service d'Astrophysique, 91191 Gif-sur-Yvette Cedex, France
2 - AIM - Unité Mixte de Recherche CEA, CNRS, Université Paris VII, UMR No. 7158, France
3 - CEA-Saclay, DSI/DIR, 91191 Gif-sur-Yvette Cedex, France

Received 24 April 2006 / Accepted 8 November 2006

Abstract
Aims. We present a new three-dimensional radiation hydrodynamics code called HERACLES that uses an original moment method to solve the radiative transfer.
Methods. The radiation transfer is modelled using a two-moment model and a closure relation that allows large angular anisotropies in the radiation field to be preserved and reproduced. The radiative equations thus obtained are solved by a second-order Godunov-type method and integrated implicitly by using iterative solvers. HERACLES has been parallelized with the MPI library and implemented in Cartesian, cylindrical, and spherical coordinates. To characterize the accuracy of HERACLES and to compare it with other codes, we performed a series of tests including purely radiative tests and radiation-hydrodynamics ones.
Results. The results show that the physical model used in HERACLES for the transfer is fairly accurate in both the diffusion and transport limit, but also for semi-transparent regions.
Conclusions. This makes HERACLES very well-suited to studying many astrophysical problems such as radiative shocks, molecular jets of young stars, fragmentation and formation of dense cores in the interstellar medium, and protoplanetary discs.

Key words: radiative transfer - scattering - hydrodynamics - methods: numerical

1 Introduction

Since we receive almost all the information from astrophysical objects in the form of photons, there has been a long tradition of radiative transfer in the astrophysical community. However, most of the time, radiation is used as a probe of the physical system under study. For that purpose radiative transfer is treated in great detail, using complex spectral opacities and without assuming local thermodynamical equilibrium (Kurucz 1996; Hauschildt et al. 1997). However, in many cases, radiation cannot be considered as a passive probe of the physical system, but it has to be considered as an important dynamical constituent.

For this reason, and thanks to the ever-increasing computing power available, there has been a rise in the interest in radiation-hydrodynamics, that is, in the dynamical coupling of gas and radiation. However, using the radiative transfer models developed in the context of spectral analysis for radiation-hydrodynamics is far out of the reach of present computers. Therefore, in order to study the dynamics of multi-dimensional systems including radiative transfer, one has to greatly simplify the full transfer equation, which depends upon six variables in a three-dimensional problem.

Our objective is to propose a radiative transfer model that can preserve and reproduce large angular anisotropies in the radiation field, which properly treats the time dependence of the radiation field and which couples naturally to the modern high-resolution numerical schemes used for hydrodynamics. All this, of course, has to be obtained at a reasonable numerical cost.

Different physical approximations have been developed to model radiative transfer in particular cases. Within the limit of large optical thickness, the diffusion approximation can be used (Turner & Stone 2001; Mihalas & Mihalas 1984; Dai & Woodward 1998). On the other hand, the transport limit is reached for transparent media. Particular methods have been implemented to describe this regime (Dai & Woodward 2000).

In many problems of physical interest, regions of large opacities are found next to transparent regions. One possibility for treating these problems is to couple two approximate models, i.e. one for each region. However, this introduces large drawbacks due to the domain partition, some loss in accuracy in the transition zone; and most of the time the semi-transparent regions would not be described correctly. Monte-Carlo codes that directly solve the transfer equation (Mihalas & Mihalas 1984; Pascucci et al. 2004, and references therein) can describe both regimes; however, they are difficult to couple to hydrodynamical codes and are very costly, especially in the diffusion regime.

Other techniques have therefore been developed based on discretization both in angles and space. One of them consists in choosing a set of discrete directions and in computing the integral over solid angle by a weighted sum over these directions. This technique, called the discrete ordinates method, can solve radiative transfer problems with relatively good accuracy and moderate computing cost. Although it has been greatly extended since Chandrasekhar first introduction, it still has two major drawbacks: ray effects and false scattering (Coelho 2002). False scattering (equivalent to the false diffusion in hydrodynamics) is due to the spatial discretization and may be reduced by refining the grid. The ray effect appears in free-streaming regions and consists in anomalous distortions of radiation intensity. Moreover, this method cannot properly treat the specularly reflected beams because of the choice of discrete directions.

To overcome this problem, both the short or the long characteristic methods choose the photons trajectories as privileged directions (Rukolaine et al. 2002). In these methods, the radiative transfer equation is integrated over ray propagation. In the long characteristic algorithm, the ray is projected from each grid point to the domain boundaries where specific intensity is known via boundary conditions. This method is unfortunately very time-consuming and not well-suited to three-dimensional problems. That is why, in the short characteristic method, the rays are projected only to the neighboring cells (van Noort et al. 2002). This is much less time-consuming, but the interpolations needed along grid lines make it more diffusive than the long characteristic one.

Instead of trying to solve the full transfer equation as in the Monte-Carlo, the discrete ordinates, or the characteristic methods, it is also possible to use another class of methods that consists in solving an approximate simplified model. This class is dedicated to model situations where radiative transfer is strongly coupled with other phenomena (fluid motion, chemical reactions, etc.) as is the case in our domain of interest. These methods consider the moments of the radiative transfer equation and consist in choosing a closure relation to solve them. The most common of these methods is the flux-limited diffusion, which solves the evolution of the first moment (radiative energy) and uses a closure relation valid in the diffusion limit, which is an isotropic radiative pressure tensor. In this scheme, the flux is always colinear and proportional to the gradient of radiative energy. In addition, the equation is modified with an ad-hoc function (the flux limiter) in order to ensure that the radiative flux remains in physically acceptable limits. This method is very useful in diffusive regions because it gives good results at a reasonable computational cost. Nevertheless, it should be used with care when dealing with free-streaming regions. Another method of closing the system is the variable tensor Eddington formalism (VTEF). It solves, in a first step, the moment equations with a fixed Eddington tensor and then computes the new tensor by solving the transfer equation locally with a fixed source function (Stone et al. 1992; Gehmeyr & Mihalas 1993) or with an approximate time dependence (Hayes & Norman 2003). If this procedure is iterated until the value of the source converges, one has a means of solving the transfer equation exactly. The VTEF methods give better results than the flux-limited diffusion but are much more complex because they require the local resolution of the transfer equation at each time step.

We have chosen in this work to use a moment model to describe radiative transfer but with a more general closure relation than the flux-limited diffusion, but one that is nevertheless analytical. This model allows treatment of radiative transfer from the diffusion to the free-streaming regime at, as we will show, a relatively low cost and with good accuracy. Furthermore, since we use a moment model, elastic scattering can be taken into account at almost no extra cost and with a reasonable accuracy, which is not the case for the other methods described before. In the next section, we present the moment model and our chosen closure. Section 3 describe the numerical scheme and its actual parallel implementation. The next section presents some multi-dimensional tests to verify our code both on purely radiative problems and on radiation hydrodynamics, as well as some performance analysis. Finally, the last section presents a summary and concluding remarks.

   
2 The physical model

For most problems and for the foreseeable future, solving the full transfer equation in three dimensions is much too costly both in time and memory. This problem is even greater if one wants to couple the radiative transfer with hydrodynamics to study multi-dimensional time-dependent problems. For this reason, we decided to use a moment model that is much less expensive than the full transfer equation and couples naturally to hydrodynamics.

The radiative transfer equation for the specific intensity is

 

$$\begin{array}{ccc} \left(\frac{1}{c} \frac{\partial}{\partial t}+... ...g(\vec{ n},\vec{ n'}) I(\vec{ x},t;\vec{ n'},\nu) {\rm d} \vec{ n'} \end{array}$$ (1)

where c is the speed of light, $\sigma^{\nu}_{\rm a}$ the absorption coefficient, $\sigma^{\nu}_{\rm s}$ the scattering coefficient, $\sigma^{\nu} =\sigma^{\nu}_{\rm s} + \sigma^{\nu}_{\rm a}$ the total cross section, B the black body specific intensity, and g the scattering angular redistribution function. $\vec{ n}$, $\vec{ r}$, and t are the angular, spatial, and temporal variables.

Integrating Eq. (1) and Eq. (1) $\times \vec{n}$over solid angle yields

 

$$\left\{ \begin{array}{lclcl} \partial_t E_{\rm r}^{\nu} &+& \nabl... ...a^{\nu} - g_1\sigma_{\rm s}^{\nu} ) c \vec{ F}_{\rm r}^{\nu} \end{array}\right.$$ (2)

where g₁ is the first moment of the angular redistribution function. Here, $E_{\rm r}, \vec{ F}_{\rm r}$, and ${\mathbb P}_{\rm r}$ are respectively the radiative energy density, the radiative energy flux, and the radiative pressure, which are defined in terms of the zeroth, first, and second moments of the specific intensity as:

 

$$\begin{array}{lcr} E_{\rm r}^{\nu} = \frac{1}{c}& \oint_{4\pi}& \... ...pi}& \vec{ n}\vec{ n}\; I(\vec{ x},t;\vec{ n},\nu) {\rm d}\vec{ n}. \end{array}$$ (3)

In this paper we concentrate on gray moment model, but the work we present can easily be generalized to multigroup radiative transfer (Turpault 2005). Assuming an isotropic scattering (g₁=0) and averaging the system (2) over frequency gives

 

$$\left\{ \begin{array}{lclcl} \partial_t E_{\rm r} &+& \nabla \cdo... ...r} &=& -c ( \sigma_{\rm F} +\sigma_{\rm s}) \vec{ F}_{\rm r} \end{array}\right.$$ (4)

where $\sigma_{\rm P}$ (respectively $\sigma_{\rm E}$ and $\sigma_{\rm F}$) are the means of $\sigma_{\rm a}^{\nu}$ weighted by the Planck function (respectively the radiative energy and flux) and $\sigma_{\rm s}$ the mean of the $\sigma_{\rm s}^{\nu}$ weigthed by the radiative flux. Therefore, if the frequency-dependent opacity coefficients ( $\sigma_{\rm a}^{\nu}$ and $\sigma_{\rm s}^{\nu}$) are known, as well as the underlying specific intensity (which is the case in our model, see below), these means can be computed precisely.

Unfortunately, the previous system is not closed since the radiative pressure is not known. The pressure is then often expressed as ${\mathbb P}_{\rm r} = {\mathbb D} E_{\rm r}$ where ${\mathbb D}$ is the Eddington tensor. There are many physical approximations that can be used to determine the Eddington tensor. The simplest one is to use the Eddington tensor given by the diffusion limit: ${\mathbb D} ={\mathbb I}/3$, where ${\mathbb I}$ is the identity matrix. This corresponds to an isotropic radiation field. We also chose an analytical Eddington tensor in order to keep the method simple and cost competitive. But, as we will see, the underlying photon distribution function (or equivalently the specific intensity) is not isotropic, which makes the method applicable to a wide range of conditions.

In order to close system (4), we want to express the Eddington tensor, and therefore the underlying specific intensity, only in terms of $E_{\rm r}$ and $\vec{ F}_{\rm r}$. It is then natural to suppose that the direction of the radiative flux is an axis of symmetry of the local specific intensity. With this assumption, the Eddington tensor is given by (Levermore 1984)

 

$${\mathbb D}=\frac{1-\chi}{2}{\mathbb I}+\frac{3\chi-1}{2} \vec{ n} \otimes \vec{ n}$$ (5)

where ${\mathbb I}$ is the identity matrix, $\chi$ is called the Eddington factor, and $\vec{ n}$ a unit vector aligned with the radiative flux.

We now need to specify the Eddington factor in order to close the system. For that purpose we assume that the specific intensity of our radiative transfer model is given by

$$B (\nu,f,T^*) = \frac{2 h \nu^3}{c^2} \left[\exp \left(\frac{h \n... ...\Vert^2}}{\Vert\vec{ f}\Vert^2}\vec{f}.\vec{\Omega}\right)\right)-1\right]^{-1}$$ (6)

where

$T^*\!=\!\frac{2}{\Vert\vec{ f}\Vert}\left(-1+\sqrt{4-3\Vert\vec{ f}\Vert^2}\ri... ...-3\Vert\vec{ f}\Vert^2}}\left(\frac{E_{\rm r}}{a_{\rm r}}\right)^{\frac{1}{4}}$ and $\vec{ f}=\frac{\vec{ F}_{\rm r}}{c E_{\rm r}}$ is the reduced flux. Note that by definition of $E_{\rm r}$ and $\vec{ F}_{\rm r}$, we have $\Vert\vec{ f}\Vert \le 1$, which means that radiative energy is transported at most at the speed of light. This distribution can either be obtained by applying a Lorentz transform to an isotropic one (Levermore 1984) or by minimizing the radiative entropy (Dubroca & Feugeas 1999). As we shall see, this simple assumption for the specific intensity allows us to compute the Eddington factor easily and analytically, but one has to keep in mind that there might be some cases where assuming such a simplified geometry may be a very poor approximation.

In the M₁ model (Ripoll et al. 2001; Dubroca & Feugeas 1999), the previous form of the specific intensity is used to compute the Eddington factor, therefore closing system (4):

$$\chi = \frac{3+4\Vert\vec{ f}\Vert^2}{5+2 \sqrt{4-3\Vert\vec{ f}\Vert^2}}\cdot$$ (7)

We can see that this closure relation recovers the two asymptotic regimes of radiative transfer well. In the free-streaming limit (i.e. transparent media), we have $\Vert\vec{ f}\Vert=1$, $\chi=1$ and ${\mathbb D}= \vec{ n} \otimes \vec{ n}$. On the other hand, in the diffusion limit, $\Vert\vec{ f}\Vert=0$, $\chi=1/3$ and ${\mathbb D}= \frac{1}{3} {\mathbb I}$, which corresponds to an isotropic radiation pressure. Equivalently, we can get these two limits directly with the specific intensity. When $\Vert\vec{ f}\Vert=0$, $T^*=\left(\frac{E_{\rm r}}{a_{\rm r}}\right)^{\frac{1}{4}}$ and the distribution is a Planck function, whereas when $\Vert\vec{ f}\Vert=1$, it tends to a Dirac in the direction of $\vec{ f}$.

   
3 Radiation transport in a static fluid

We first start by describing the interaction between a static fluid and radiation. The fluid can only be heated or cooled and its evolution is determined by the classical energy conservation equation with a source term characterizing the energy exchanges between the fluid and the radiation.

To ensure a good conservation of energy, we consider the equation for the total energy (radiation plus matter) instead of only matter energy. The system to be solved is:

 

$$\left\{ \begin{array}{lcccccl} \partial_t e &+& \partial_t E_{\rm... ...} & = & -(\sigma_{\rm F} +\sigma_{\rm s}) c \vec{ F}_{\rm r} \end{array}\right.$$ (8)

where e is the internal matter energy.

   
3.1 The Riemann solver

To numerically solve the system above, we use a second-order Godunov type algorithm for the hyperbolic subsystem formed by the last two equations. The first equation is then integrated using the flux obtained by the hyperbolic solver.

It is worth noticing that the wave speeds of this subsystem, which mathematically correspond to the eigenvalues of its Jacobian matrix, depend only on the norm of the reduced flux $\vec{ f}$ and on the angle $\theta$ of this flux with the considered interface. These wave speeds describe the speed at which the information is transported in the system, in the same way as the sound speeds in a fluid at rest. Figure 1 illustrates the behavior of the eigenvalues normalized by c for some characteristic values of $\theta$ and $\vec{ f}$.

	Figure 1: Eigenvalues of the Jacobian matrix normalized by c.
Open with DEXTER

The left plot corresponds to a flux perpendicular to the interface ( $\theta = 0$), which is similar to the mono-dimensional problem (cf. Fig. 1 of Audit et al. 2002). In particular, for a unit reduced flux (points A), the four eigenvalues are equal to c so the transport limit is correctly described. The middle plot represents the case where the flux is parallel to the interface. In that particular case, two eigenvalues are always equal to zero and the two others are equal in norms but of opposite sign. We can also notice that, when the reduced flux is unity (points C), the four eigenvalues are null. This is particularly interesting because it corresponds to the physical characteristic speeds of a transverse flux and because it means that there is no transport perpendicular to the radiative flux (cf. shadow test below). In all cases, we find that for $\Vert\vec{ f}\Vert=0$ (points B), the eigenvalues are { $-c/\sqrt{3},0,0,c/\sqrt{3}$}, which are the proper propagation speeds in the diffusion limit. These eigenvalues are then used in our Riemann solver, which is an HLLE (Harten-Lax-van Leer-Einfeldt) scheme (Einfeldt et al. 1991), see Appendix A for details.

   
3.2 The system solver

Our scheme is of order two in space. It means that all quantities evaluated at the interfaces are reconstructed by using a linear reconstruction in each cell. To insure numerical stability, the results presented in this paper used a minmod limiter, but another limiter could easily be used instead. HERACLES has been implemented in order to work in either Cartesian, cylindrical, or spherical geometry (cf. Appendix B for the specific divergence discretization in these geometries).

Noticing that the time step given by the Courant condition is much smaller for radiation than for hydrodynamics, and since we are most often interested in studying system over hydrodynamical time scales, we need to develop an implicit algorithm for the radiation. The time step is then given by the explicit Courant condition for hydrodynamics and by controlling the variation of the variables for the implicit radiative solver. When doing radiation-hydrodynamics, the smallest of these two time steps is taken (see Appendix C for details).

We now have to solve a non-linear set of equations (the non-linearities come from the Eddington factor and the T⁴ terms). This is done using a Raphson-Newton method, which can eventually be restricted to a single iteration if solving the linearized system is enough. Either way, we have to invert a $(2+d) N^d\times (2+d) N^d$matrix, where d is the dimensionality of the problem and N the number of cells along a direction. The typical scale of a simulation makes a direct inversion out of reach, so we need to resort to an iterative inversion method. We have implemented and tested two different methods that are relevant in different physical situations, each one efficient in one of the two limits of the transfer equation.

The Gauss-Seidel method (see Press et al. 1986, for further details) is very efficient in the transport limit. For example, in one dimension if we know a priori the direction in which the information propagates, the algorithm converges theoretically in only one iteration. To take advantage of this feature, we can also use a sweep method that inverts the cells ordering at each iteration. All the tests performed show that this method is very efficient both in time and in memory requirements. One could be tempted to use the improved method called successive overrelaxation (SOR), but it should be taken with care. Indeed, in some cases, too large extrapolation could lead to a non physical iterate that does not respect the condition $\Vert f\Vert \le 1$ and that breaks down the iterations because of the M₁closure.

On the other hand, the generalized minimal residual (GMRES) method (Saad & Schultz 1986), although more constraining in terms of memory requirement, is well-suited to the diffusion limit in which case the matrix is diagonally dominant. As each method dumps the residual in a different way, we also implemented a coupled version for which we switch between both methods. This version seems to have the best convergence rate (cf. Fig. D.2 and Appendix D for further details).

For both inversion methods, we need a convergence criterion. After trying various possibilities, we found that using the residual of the adimensioned equations gave the best results. In order to adimensionate the equation on the radiative flux, we divide it by the radiative energy times the speed of light and not by the radiative flux, which can be vanishing. In all the tests presented in this paper, the convergence is assumed to be achieved when this residual is less than 10^-5.

As our code is parallelized with the MPI library, we also quantified the scaling of these two methods as the number of processors increased. We found that both of them show good scaling; however, the Gauss-Seidel algorithm is nearly perfect, whereas GMRES could have a loss in performance reaching 20% (cf. Appendix D).

4 Radiation hydrodynamics in a moving fluid

We now generalize the previous analysis to a moving fluid. The fluid evolution is determined by the classical conservation equations (mass, momentum, and energy) with source terms characterizing the momentum and energy exchanges between the fluid and the radiation.

In order to write the radiation hydrodynamics equations, one has to choose the frame in which to evaluate the radiative quantities: laboratory frame or comoving frame (i.e. the frame moving with the fluid). The laboratory frame is convenient because the system remains globally conservative, which keeps the hyperbolic part (i.e. the left-hand side) of the system simple (Mihalas & Auer 2001). But in this frame, interactions with matter become complex because of Doppler and aberration effects that have to be incorporated in the source terms. On the other hand, using the radiative quantities expressed in the comoving frame (Lowrie et al. 2001) adds non-conservative terms to the equations. But the source terms remain unaffected by the fluid motions.

In HERACLES, we have chosen to express radiative quantities in the comoving frame. The equations of radiation hydrodynamics can then be written (Mihalas & Mihalas 1984) as

 

$$\left\{ \begin{array}{lclcl} \partial_t \rho & + & \nabla \cdot [... ...gma_{\rm F}+\sigma_{\rm s}}{c} \vec{ F}_{\rm r}\cdot\vec{ u} \end{array}\right.$$ (9)

 

$$\left\{ \begin{array}{lccccclll} \partial_t E_{\rm r} &+& \nabla ... ...&&&&& = -(\sigma_{\rm F} +\sigma_{\rm s}) c \vec{ F}_{\rm r} \end{array}\right.$$ (10)

where $\rho$ is the matter density, $\vec{ u}$ the velocity, E the total matter energy, P and T the pressure, and the temperature of the material. The previous equations are all Eulerian, but the radiative quantities are evaluated in the frame comoving with the fluid.

The resolution of the previous system is split into three steps. The first one updates the hydrodynamical quantities using a classical second-order Godunov type method. In the next step, the radiative quantities are updated implicitly, as described in the previous section. During this step, we use the velocity given by the hydro solver. In the third and final step, the $(\sigma_{\rm F} + \sigma_{\rm s}) F/c$source term is added to the gas momentum and total energy. Since in most cases this term is rather small, it is not necessary to treat it in the implicit system. This allows the reduction of the number of equations to be solved implicitly and makes the method more efficient:

 

$${\rm STEP \ 1 \;} \left\{ \begin{array}{lclclcl} \partial_t \rho ... ... \\ \partial_t E &+& \nabla \cdot [\vec{ u}(E+P)] & &&= & 0 \end{array}\right.$$ (11)

 

$${\rm STEP \ 2 \;}\left\{ \begin{array}{llllllllll} \partial_t e +... ...u & = & -(\sigma_{\rm F} +\sigma_{\rm s}) c \vec{ F}_{\rm r} \end{array}\right.$$ (12)

 

$${\rm STEP \ 3 \;}\left\{ \begin{array}{lcl} \partial_t(\rho \vec{... ...ma_{\rm F}+\sigma_{\rm s}}{c} \vec{ F}_{\rm r}\cdot\vec{ u}. \end{array}\right.$$ (13)

   
5 Verification tests

We performed a series of verification tests in order to better characterize the accuracy of HERACLES and to compare it with other codes. We reproduce here just the three that are particularly relevant to astrophysical conditions. Appendix E compiles others tests: beam test, pipe flow test, diffusion in a moving fluid, and a comparison with simple analytical models of matter-radiation coupling and Marshak wave.

5.1 Shadow test

We performed a 2D shadow test similar to the one presented in Hayes & Norman (2003). This test consists in lighting an oblate spheroid clump in a cylinder L=1 cm, R=0.12 cm. This spheroid is located on the symmetric axis and at the centre of the box width: $(z_{\rm c},r_{\rm c})=(0.5,0)$. Its extension is (z₀,r₀)=(0.1,0.06).

Initially, the medium is at equilibrium with radiation i.e. $T_0=T_{\rm r}=290 \; {\rm K}$, and the density is homogeneous ( $\rho_0=1 \ {\rm g}\ {\rm ~cm}^{-3}$) except for the clump with density $\rho_1$one thousand times greater. The boundary of this region is smoothed, such as $\rho(z,r)=\rho_0 +\frac{\rho_1-\rho_0}{1+\exp{\Delta}}$ with $\Delta=10 \left[ (\frac{z-z_{\rm c}}{z_0})^2 + (\frac{r-r_{\rm c}}{r_0})^2 -1 \right]$. The opacity of the medium is a function of density and temperature $\sigma=\sigma_0 (\frac{T}{T_0})^{-3.5}(\frac{\rho}{\rho_0})^2$ with $\sigma_0 = 0.1 \ {\rm ~cm}^{-1}$. Initially, the mean free path is then 10 cm in the cylinder and 10^-5 cm in the clump. At time t=0, a uniform source is lighted at the left boundary with $T_{\rm r}=1740 \; {\rm K}$. The mean free path of a photon being much smaller in the clump (by six orders of magnitude), a shadow develops behind it. Until the light has crossed the clump, the shadow should remain stable.

Figure 2 shows the radiative temperature for three runs performed on a $280\times 80$ grid. The first run (upper panel) was performed solving the diffusion equation. It corresponds to the time t=0.1 s, which means $3\times 10^9$ light crossing times. The other two runs used the M₁ model and differs by the method of computing the eigenvalues. In the first case (middle panel), we chose not to compute these eigenvalues and to set them arbitrarily equal to $\pm c$, whereas in the second run (lower panel), the eigenvalues were computed. To compare with the results obtained in Hayes & Norman (2003), we also plotted a radial profile of the radiative temperature at the outer boundary (cf. Fig. 3).

One can see that there is no longer a shadow in the first case. This is due to the fact that the diffusion equation is intrinsically isotropic ( $\mathbb{P}=1/3 E \mathbb{ I}$). Therefore, the photons can go around the obstacle immediately. In the two other cases, the shadow is better preserved even in the simple case without computing eigenvalues. As could be expected, the first method is more diffusive than the second one. The improvement obtained with the second method is easy to understand when looking at the real eigenvalues (cf. Sect. 3.1). The proper treatment of the propagation speed in the HLLE scheme is effective enough to inhibit large numerical diffusion between the shadowed and enlightened regions.

It is also interesting to compare the propagation speed of the radiation between the diffusion approximation and the M₁ model. For the diffusion, the characteristic length scales as the square root of time: $L \propto \sqrt{2 \alpha t}$ where $\alpha$ is the diffusion coefficient. In our example, it means that radiation will cross the box in a few 10^-9 s, whereas in reality as radiation propagates in a transparent medium, it goes at light speed and the crossing time is $3.33\times 10^{-11}$ s. In the two M₁ model simulations, we recover this value.

	Figure 2: Radiative temperature in the shadow test using the diffusion equation ( upper panel), M1 closure with fixed eigenvalues ( middle panel), and M1 closure with calculated eigenvalues ( lower panel).
Open with DEXTER

	Figure 3: Radial profiles of the radiative temperature in the shadow test using the diffusion equation (dotted line), M1 closure with fixed eigenvalues (dashed line), and M1 closure with calculated eigenvalues (solid line).
Open with DEXTER

5.2 Semi-transparent regime

5.2.1 Without scattering

We now present a test of HERACLES in a semi-transparent regime where the validity of the physical model used to treat the transfer has to be assessed. The following test can be seen as a star with solar luminosity illuminating its surrounding disc of gas at a distance approximately equal to 5 AU (which corresponds to the Sun-Jupiter distance).

The simulation box, whose dimensions are $L_x=7.48 \times 10^{12} \ {\rm ~cm}$ and $L_y=3.74 \times 10^{12} \ {\rm ~cm}$, is at equilibrium, and at t=0 an incoming horizontal radiative flux equal to $5.44 \times 10^{4} \ {\rm erg} \ {\rm s}^{-1} \ {\rm ~cm}^{-2}$ is set. The box is filled with matter at density  $10^{-15} \ {\rm g}\ {\rm ~cm}^{-3}$, while there is a vacuum on the outside. The width of the box corresponds to seven mean free paths (the absorption coefficient is $\kappa=10^{-12} \ {\rm ~cm}^{-1}$) and is sampled over 50 cells. This box is therefore a semi-transparent region, and this run will test the accuracy of the M₁ closure in this particular transient regime. Figure 4 shows the isotherms obtained at equilibrium by HERACLES and by a Monte-Carlo code (Dullemond & Natta 2003).

The two results agree with good precision. The differences could be due to the M₁ model itself or to the slightly different treatment of the boundary conditions. It is important to note that the Monte-Carlo code solves the transfer equation exactly and, therefore, the agreement between these two approaches is a very good test to verify our model, both physically and numerically, in this semi-transparent regime.

	Figure 4: Left panel: HERACLES temperature map. Right panel: isotherms for HERACLES (solid lines) and a Monte-Carlo code (dashed lines).
Open with DEXTER

5.2.2 With scattering

Scattering is a phenomenon that plays an important role in many astrophysical problems, such as in the interstellar medium where light is scattered by dust grains. The radiative model used in HERACLES allows us to treat scattering at no additional cost and with reasonable accuracy.

By definition, the importance of scattering compared to absorption is quantified by the albedo:

$$\omega = \frac{\sigma_{\rm s}}{\sigma_{\rm F} + \sigma_{\rm s}}\cdot$$

When $\omega = 0$, the medium is purely absorptive and, when $\omega =1$, it is purely scattering. We have already seen (Eq. (4)) that the scattering only appears in the equation on the flux and not in the one on the energy. This is understood because elastic scattering redistributes the photons without changing their frequencies.

We wanted to perform the same test as above but with scattering; but to have results to compare with, we changed the scale of the problem and borrowed initial conditions from the engineering literature (Crosbie & Schrenker 1984). We consider a two-dimensional box whose dimensions are $L_x={\tau_x}_0$ and $L_y=2 {\tau_y}_0$ where $\tau_x=(\sigma_{\rm F} +\sigma_{\rm s})x$ and $\tau_y=(\sigma_{\rm F} +\sigma_{\rm s})y$ are the optical coordinates. This box is lit from left by a source at $T_{\rm s}$ and the steady state is achieved.

We reproduce the results of only one test here, but several others were performed and show good agreement, too. We took ${\tau_x}_0=1$, ${\tau_y}_0=5$, and each direction is sampled over 100 cells. We chose this particular test because it reproduces a rectangular medium for which the influence of albedo is very important. Figure 5 shows the normalized radiative energy $(T_{\rm r}/T_{\rm s})^4$ isocontours at five values 0.3, 0.4, 0.5, 0.6, and 0.7, for $\omega =1$ and $\omega =0.9$. This figure is very similar to Fig. 12c of Crosbie & Schrenker (1984). As expected, the isocontours penetrate deeper in the domain when there is no absorption. Radiative energy is greater when absorption exists.

The differences in terms of distance covered between our results and those of Crosbie & Schrenker (1984) remain lower than 5%. But for HERACLES, contrary to other methods, the treatment of scattering is nearly free in terms of computational cost. Indeed, in our case, taking scattering into account implies slightly modifying the source term in the radiative flux equation, but it does not change the number of operations needed at each iteration of the implicit solver. On the other hand, when solving the transfer equation, treating scattering implies computing an integral over solid angles of the specific intensity, which is very costly.

	Figure 5: Normalized radiative energy isocontours for $\omega =1$ (solid lines) and $\omega =0.9$ (dashed lines).
Open with DEXTER

5.3 Radiative shocks

Radiative shocks might be encountered in various astrophysical systems: stellar accretion shocks, pulsating stars, interaction of either jets or supernovæ with the interstellar medium, etc. They are also reproduced in laboratory experiments, in a scaled manner, in order to better understand their underlying physics (González et al. 2006; Bozier et al. 1986).

There are two kinds of shocks depending on the initial fluid speed: subcritical and supercritical (Mihalas & Mihalas 1984). When the initial fluid speed is low, the matter temperature before the shock is much lower than the matter temperature behind it and the transition is sharp. The radiative temperature never overshoots the temperature behind the shock. Such shocks are called subcritical. On the other hand, when the initial fluid speed is high enough, the temperature on each side of the shock are equal. The transition is smoother, the radiative precursor is larger, and a matter temperature peak, which extends over roughly one mean free path, appears just behind the shock. This is a supercritical shock. The specific profiles and characteristics of these shocks are the result of the strong coupling between matter and radiation.

We consider a 1D homogeneous medium where the fluid moves uniformly from right to left and the left boundary is a wall. A shock is therefore generated at this boundary and travels the box from left to right. The initial conditions are such that $L_x=7 \times 10^{10} \; {\rm ~cm}$ (divided in 300 cells) and $\rho = 7.78 \times 10^{-10} \; {\rm g ~cm^{-3}}$. The gas is perfect with an adiabatic coefficient of 7/5 and temperature at 10 K in equilibrium with radiation. The medium is supposed to have a constant extinction coefficient: $\sigma = 3.1 \times 10^{-10} \; {\rm ~cm^{-1}}$.

In all the above figures, we plot the temperatures as function of x_i=x-u t to be able to compare our results with those obtained in Ensman (1994) and in Hayes & Norman (2003).

5.3.1 Subcritical shocks

We perform a first test with an initial speed $u= - 6 \; {\rm km ~ s^{-1}}$ so as to have a subcritical shock (cf. Fig. 6). As expected, radiation and matter are not in equilibrium: temperatures differ upstream and downstream. We perform two series of tests: one with the M₁ model and another one with a constant isotropic Eddington tensor ${\mathbb D} ={\mathbb I}/3$ (also known as the P₁ model). We can see that in both cases the shock is located at the same place and the temperatures behind the shock are very similar. However, some differences appear in the precursor that have the greatest effect on the radiative temperature. The radiative precursor appears larger in the M₁ model due to the M₁ ability to deal with large reduced flux and anisotropic photon distribution function. In the precursor, the reduced flux is very large with a value around 0.9. But, in the P₁ model, the reduced flux is limited by 1/$\sqrt{3}$. The M₁ model, in which the reduced flux can go from 0 to 1, is therefore more adequate for describing this region.

	Figure 6: Subcritical shock: gas (solid line) and radiative (dashed line) temperatures at 1.7e4, 2.8e4, and 3.8e4 s obtained with the M1 model (thick lines) compared to a constant and isotropic Eddington tensor model (thin lines).
Open with DEXTER

5.3.2 Supercritical shocks

After setting an initial speed to $u= - 20 \; {\rm km ~ s^{-1}}$, we obtain a supercritical shock (cf. Fig. 7). In that case, the radiative precursor is indeed larger than in the subcritical case. Radiative and matter temperatures are equal on both sides of the shock, showing that matter and radiation are in equilibrium in a larger zone. The peak behind the shock in the matter temperature profile is sampled over 4-5 cells, which means about three tenths of the photon mean free path. In the precursor, matter and radiation are in equilibrium over a large zone, implying a small reduced flux. It is only at the end of the radiative precursor that the reduced flux becomes large and that matter and radiation temperature differ significantly.

	Figure 7: Supercritical shock: temperatures (thick lines) of gas (solid line) and radiation (dashed line) at 4.0e3, 7.5e3, and 1.3e4 s and the corresponding reduced fluxes (thin lines).
Open with DEXTER

The ability of HERACLES to treat these radiative shocks properly shows that our resolution method divided in three steps (Eqs. (11)-(13)) does not alter the coupling between matter and radiation.

   
6 Summary and discussion

We have presented a new, parallelized, three-dimensional radiation hydrodynamics code called HERACLES. The radiative transfer is solved using a moment method and the M₁ closure relation. We have shown that, even though approximate, this model gives an accurate solution even when the radiation field presents large angular anisotropies. The numerical scheme proposed for integrating the M₁ model is implicit and allows us to numerically keep the correct physical properties of M₁.

The tests performed and presented in this paper reveal the accuracy of HERACLES both in the diffusion and transport limit, and despite all the approximations, it compares well to codes that solve the exact transfer equation. In particular, and in opposition to flux limited diffusion, HERACLES can deal with semi-transparent and free-streaming regions. The possibility of easily treating scattering with reasonable accuracy is also an advantage to our method. Finally, it naturally couples to the modern high-resolution numerical schemes used for hydrodynamics, and its reasonable numerical cost allows the treatment of three-dimensional radiation hydrodynamics problems.

The tests presented in this paper show that HERACLES will be very well-suited to studying various astrophysical problems from the study of protoplanetary discs to the formation and fragmentation of dense cores in the interstellar medium. It has already been used to study and interpret multidimensional effects in laboratory astrophysics radiative shocks experiments carried out at the PALS laser (González et al. 2006). And a study of the interaction of molecular jets with the interstellar medium is in preparation (González et al. 2007).

Some further developments of the code are also envisaged. The most useful for astrophysical applications are implementation of a multigroup M₁ model and coupling to an AMR grid to obtain a very high resolution while remaining cost competitive.

Acknowledgements

HERACLES is available at http://www-dapnia.cea.fr/Projets/COAST/heracles.htm or by sending an e-mail to edouard.audit@cea.fr
We would like to thank the CEA computing center, CCRT, where all the computations presented in this paper were done.

References

• Audit, E., Charrier, P., Chièze, J. P., & Dubroca, B. 2002, [arXiv:astro-ph/0206281] (In the text)
• Bozier, J. C., Thiell, G., Le Breton, J. P., Azra, S., & Decroisett, M. 1986, Phys. Rev. Lett., 57, 1304 [NASA ADS] [CrossRef]
• Coelho, P. J. 2002, JQSRT, 73, 231 [NASA ADS] (In the text)
• Crosbie, A. L., & Schrenker, R. G. 1984, JQSRT, 31, 339 [NASA ADS] (In the text)
• Dai, W., & Woodward, P. R. 1998, J. Comput. Phys., 142, 182 [CrossRef]
• Dai, W., & Woodward, P. R. 2000, J. Comput. Phys., 157, 199 [NASA ADS] [CrossRef] (In the text)
• Dubroca, B., & Feugeas, J. L. 1999, CRAS, 329, 915 (In the text)
• Dullemond, C. P., & Natta, A. 2003, A&A, 408, 161 [NASA ADS] [CrossRef] [EDP Sciences]
• Einfeldt, B., Munz, C. D., Roe, P. L., & Sjögreen, B. 1991, J. Comput. Phys., 92, 273 [NASA ADS] [CrossRef] (In the text)
• Ensman, L. 1994, ApJ, 424, 275 [NASA ADS] [CrossRef] (In the text)
• Gehmeyr, M., & Mihalas, D. 1993, A&AS, 25, 1366
• Gentile, N. A. 2001, J. Comput. Phys., 172, 543 [NASA ADS] [CrossRef] (In the text)
• González, M., Stehlé, C., Audit, E., et al. 2006, Laser Part. Beams, 24, 535
• González, M., Audit, E., & Lery, T. 2007, in preparation (In the text)
• Hauschildt, P. H., Baron, E., & Allard, F. 1997, ApJ, 483, 390 [NASA ADS] [CrossRef]
• Hayes, J. C., & Norman, M. L. 2003, ApJS, 147, 197 [NASA ADS] [CrossRef] (In the text)
• Kurucz, R. L. 1996, ASPC, 160
• Levermore, C. D. 1984, JQSRT, 31, 149 (In the text)
• Lowrie, R. B., Mihalas, D., & Morel, J. E. 2001, JQSRT, 69, 291 [NASA ADS] (In the text)
• Marshak, R. E. 1958, Phys. Fluids, 1, 24 [CrossRef] (In the text)
• Mihalas, D., & Auer, L. H. 2001, JQSRT, 71, 61 (In the text)
• Mihalas, D., & Mihalas, B. D. 1984, Foundation of Radiation Hydrodynamics (Oxford University Press)
• Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1986, Numerical recipes (Oxford University Press) (In the text)
• Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Ripoll, J. F., Dubroca, B., & Duffa, G. 2001, Combust. Theory Modell., 5, 261 [NASA ADS]
• Rukolaine, S. A., Vasilyev, M. G., Yuferev, V. S., & Galyukov, A. O. 2002, JQSRT, 73, 205 [NASA ADS] (In the text)
• Saad, Y., & Schultz, M. H. 1986, SIAM J. Sci. Stat. Comput., 7, 856 [CrossRef] (In the text)
• Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 [NASA ADS] [CrossRef]
• Su, B., & Olson, G. L. 1996, JQSRT, 56, 337 (In the text)
• Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef]
• Turpault, R. 2005, JQSRT, 94, 357 [NASA ADS] (In the text)
• van Noort, M., Hubeny, I., & Lanz, T. 2002, ApJ, 568, 1066 [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

   
Appendix A: The HLLE scheme

We consider the system:

$$\left\{ \begin{array}{ccccl} \partial_t E_{\rm r} &+&\nabla \cdot... ... & = & -(\sigma_{\rm F} +\sigma_{\rm s}) c \vec{ F}_{\rm r}. \end{array}\right.$$ (A.1)

Noting down $\mathcal U$ the vector of variables and $\mathcal F$ the corresponding flux vector, this system can be written as

$$\partial_t \mathcal U + \partial_x \mathcal F(\mathcal U) + \... ...) + \partial_z \mathcal H(\mathcal U) = \mathcal S(\mathcal U).$$

First, we have to compute the minimum and maximum eigenvalues (respectively $\lambda_{\rm min}$ and $\lambda_{\rm max}$) of the Jacobian matrix $\mathcal J(U)=\frac{\partial \mathcal F(\mathcal U)}{\partial \mathcal U}$. Then, the intercell fluxes are computed using the HLLE scheme

$$\mathcal F_{i+\frac{1}{2}}=\frac{\lambda^+_{i+\frac{1}{2}} \m... ...cal U_i)}{\lambda^+_{i+\frac{1}{2}}-\lambda^-_{i+\frac{1}{2}}}$$ (A.2)

where the index $i+\frac{1}{2}$ denotes the interface between cells iand i+1, $\lambda^+={\rm max}(0,\lambda_{\rm max})$ and $\lambda^-={\rm min}(0,\lambda_{\rm min})$.

The computation of these eigenvalues is rather time-consuming because we need them at each iteration in a time step and at each interface of the mesh. However, since they only depend on two parameters (i.e. $\Vert\vec{ f}\Vert$ and $\theta$), they can be tabulated easily. We therefore decided to compute them once for a set of $\theta$ and fand to interpolate the value needed. This method performs well because the eigenvalues have a smooth behavior. The maximum difference obtained between the exact eigenvalues and the interpolated ones never exceeds one percent using a $100\times 100$ interpolation grid.

   
Appendix B: The geometry

Table B.1: Divergence in the three different geometries handled by HERACLES.

In one Cartesian dimension, noting down the time index n and n+1, the cell index i, and the interfaces indexes $i\pm1/2$, the discretized equations we must solve are:

 

$$\left\{ \begin{array}{lllllll} \frac{e_i^{n+1}-e_i^n}{\Delta t} {... ...F}}_i^n + {\sigma_{\rm s}}_i^n ) c F_i^{n+1} {\Delta V}_i^n. \end{array}\right.$$ (B.1)

But HERACLES is able to work in either Cartesian, cylindrical, or spherical geometry, and in the two latter cases, the divergence involves terms due to geometrical effects. These terms are written in the Table B.1 where $\vec{ P}$ corresponds to the row vectors of the tensor ${\mathbb P}$, which is the radiative pressure or the total hydrodynamical pressure: $P{\mathbb I} + \rho \vec{ u}\otimes \vec{ u}$.

One should take particular care in the discretization of these terms in order to check three conditions. First, the discrete divergence of a constant must be null. Second, when r tends to infinity, the geometry effects vanish and one should recover the Cartesian divergence. Finally, when the resolution increases (i.e. dr tends to zero), the discrete scheme should not diverge. In order to verify these three conditions in HERACLES, we then discretize the "nabla'' terms with the classical formula $\int_V \nabla P {\rm d}V = \int_S P {\rm d}S= \Sigma P_S S$. For the derivative terms, we use a discretization over the direction considered involving interface values:  $\int_V \partial_{\rm r} P {\rm d}V = \frac{P_{i+1/2}-P_{i-1/2}}{r_{i+1/2}-r_{i-1/2}}V$and $\int_V \frac{1}{r} \partial_{\theta} P {\rm d}V = \frac{P_{i+1/2}-P_{i-1/2}}{r_i(\theta_{i+1/2}-\theta_{i-1/2})}V$. As for the additional terms, they are cell-centered.

This discretization ensures us that the momentum in the hydrodynamics equations is conserved. One can notice that the additive centered terms drop naturally when the diffusion limit is reached, because in such case the tensor of radiative pressure is isotropic: $P_{rr}=P_{\theta\theta}=P_{\phi\phi}$.

   
Appendix C: Time step control

The time step in HERACLES is the smallest of the hydrodynamical and the radiative time steps. The hydrodynamical time step is simply given by a Courant condition using a Courant factor of 0.8. For the radiation, we compute two time steps. An explicit one, also given by a Courant condition, and an implicit one controlled by the variation in the variables involved in the implicit step (gas energy, radiative energy, and radiative flux). We set a limit on the allowed variation for these variables during a time step, and the implicit time step is increased (resp. lowered) if the actual variations are lower (resp. higher) than this limit. The final radiative time step is the implicit one if it is at least 10 times greater than the explicit one. This last condition ensures that implicit integration is done only when it is profitable.

   
Appendix D: Performances

We present here the performances obtained with HERACLES. The code is parallelized with the MPI library in order to run on a large parallel supercomputer. The parallelization scheme used is a classical domain decomposition. This splitting ensures a good equilibrium of charge between all the processors. The communications at the boundary between each domain always take a negligible time in all the tests we have done.

To test the performances, we ran different simulations over 1, 16, 36, 64, and 121 processors such that each processor worked on $150\times 525$ cells. In Fig. D.1, we plot the CPU time per processor and per cell normalized by the Gauss-Seidel CPU time with one processor. We can clearly see that the Gauss-Seidel method has a nearly perfect scaling, whereas the GMRES method has a more complex behavior with a loss of performance that can reach 20%. This is due to the fact that the GMRES algorithm has a more complex communication scheme and is therefore more sensitive to small load unbalances. We can also see that one GMRES iteration nearly corresponds to three Gauss-Seidel iterations in terms of CPU time. Another drawback of the GMRES method is that it is very demanding in terms of memory space. In conclusion, the GMRES method is appropriate if it significantly accelerates the convergence (at least by a factor of 3) and if memory is not a crucial issue.

The convergence rate of the above algorithms can vary significantly from one test to the next. Figure D.2 shows the residual during one time step of a Marshak wave test. We can see that the Gauss-Seidel method in this particular test is far much efficient than the GMRES method. We have performed another Marshak wave test with higher opacities in order to be in a more diffusive regime. In that case, the results are inverted, with the GMRES method better than Gauss-Seidel. We recover the characteristics discussed in Sect. 3.2: Gauss-Seidel is more efficient in the transport limit and GMRES in the diffusion limit. However, we empirically discovered that coupling the two methods gives the best results. We first performed a few iterations with one method and then gave this result as initial guess to the other one. At the switch the residual drops by about one order of magnitude, due to the fact that for each method the first iterations are more efficient than what follows. Pushing this idea further, we tried the "yoyo'' method, which consists in switching back and forth between the two methods every 10 iterations. In all the tests we have done, this method was the most efficient (if the memory needed by GMRES can be afforded), even though the optimal number of iterations between each switch can vary somewhat.

In all our simulations, the CPU time needed for the radiative transfer was between 3 and 10 times the time taken by the hydrodynamics depending on the variation in the radiative quantities during one hydrodynamical time step and on the convergence criteria used. The memory needed by the radiative transfer is equal to the memory needed by the hydrodynamics if one uses the GS solver but it is higher when using the GMRES algorithm (depending on the size of the Krylov space used to find the solution).

Figure D.1: Normalized CPU time as a function of the number of processors for simulations with a constant cell number per processor.

Figure D.2: Residual in a Marshak wave test for the five following methods: 1. purely Gauss-Seidel (GS) iterations, 2. purely GMRES iterations, 3. 30 iterations of GS before GMRES, 4. 50 iterations of GMRES before GS, and 5. yoyo with a switch between GS and GMRES every 10 iterations.

   
Appendix E: Extra verification tests

E.1 Beam test

In the shadow test, the numerical diffusion is limited because radiation propagates along mesh axis. In order to quantify this efect, we performed a test where a narrow beam of radiation propagates with a certain angle (cf. Richling et al. 2001).

For this test, we computed a $128\times 128$ mesh that covers the domain x=[-1,1] and y=[-1,1]. The beam is introduced at x=-1 between y=[-0.875, -0.750] with an angle of $45^\circ$ and with a unit reduced flux. Initially, the temperature is at $T=T_{\rm r}=300 \; {\rm K}$and the beam is at $T=T_{\rm r}=1000 \; {\rm K}$. There is no scattering, absorption, or emission: $\sigma=0$. The beam should therefore cross the medium without dispersion.

If we plot a profile of the radiative energy, the beam is initially a step sampled over 8 cells. In the stationary state, we can see (cf. Fig. E.2) that, for fixed eigenvalues, the FWHM approximately corresponds to 30 cells. When the eigenvalues are computed, this FWHM falls to only 24 cells. The numerical diffusion has been therefore reduced by 20%. This test shows that the eigenvalue computation keeps the numerical diffusion under control even in a direction not along mesh axis.

Figure E.1: Radiative energy for the beam test: eigenvalues set to $\pm$c ( left) or calculated ( right).

Figure E.2: Horizontal cut in Fig. E.1 at the middle height. Solid line corresponds to calculated eigenvalues and dashed line to fixed ones.

E.2 Pipe flow test

This test, known as the pipe flow test (Gentile 2001), is not an astrophysical one, but it treats all the physics underlying the radiative transfer. It is very demanding because it involves at the same time regions of free-streaming, of diffusion, shadows, and the results depend strongly on the description of the heating and reemission of the walls.

In this problem, a cylindrical section of dense and opaque material is embedded in less dense and less opaque material (the pipe), which is itself embedded in a cylinder of the dense and opaque material. The geometry is sketched in Fig. E.3.

Figure E.3: Sketch of the geometry for the pipe flow test. The simulation box extends from r=0 to r=2 and from z=0 to z=7. The pipe can be divided in three regions. The first one is a cylinder with a radius of 0.5 which extends from z=0 to z=2.5 and from z=4.5 to z=7. The second is a cylinder with a radius of 1.5 extended from z=2.5 to z=3 and from z=4 to z=4.5. Finally the third one is a cylindrical shell with inner radius r=1 and outer radius r=1.5 which extends from z=3 to z=4. The five characteristic points are A(r=0 ; z=0.25), B(0 ; 2.75), C(1.25 ; 3.5), D(0 ; 4.25), E(0 ; 6.75).

Both materials have a heat capacity $c_v=10^{15} \; {\rm erg ~g^{-1} ~keV^{-1}}$, but the densities and opacities differ:

$$% latex2html id marker 1741 \left\{ \begin{array}{*{8}{l}l} ... ...&{\rm ~cm^{-1}} &\mbox{ inside the pipe.}\\ \end{array}\right.$$

The time unit is such that c=300. The end of the simulation is reached when $t_{\rm end}=100$. The simulation box is sampled over $400\times 1400$ cells.

Initially, the medium is at equilibrium $T=T_{\rm r}=0.05 \; {\rm keV}$everywhere. At time t=0, a uniform source with $T_{\rm s}=0.5 \; {\rm keV}$ is lighted on at the left side of the pipe.

In Fig. E.4, the problem is analyzed with the matter temperature at the five characteristic points A, B, C, D, and E. Points A and B are located in the pipe before the obstacle. Their temperature therefore depends on the code's ability to recover the free-streaming limit. These two plots show good agreement with the results published in Gentile (2001). The next three points depend closely on the physical equations solved. Indeed, in a diffusion model, photons would go round the obstacle very easily since the photon distribution function is isotropic, whereas in fact the photons have to be absorbed by the obstacle, which is heated, and then reemitted from the heated walls. This allows photons initially moving horizontally to go upward and then to go round the obstacle. For these three last points, our results are qualitatively good, although quantitatively some discrepancies exist: point E is heated faster, point C is heated a little bit later, etc. Only few people perform this test, so there are not many results to compare with. Nevertheless, we see that our model, although simple, compares well to exact codes with Monte-Carlo techniques.

Figure E.4: Matter temperature at the five characteristic points.

E.3 Diffusion in a moving fluid

This simple radiation hydrodynamics test compares diffusion in a moving fluid and in a fluid at rest. We use a Cartesian grid of $100\times 100$ cells. At t=0, we initialize a Gaussian pulse of energy at the center of the simulation box. First, we perform a test without radiation. The energy radiation is advected as a passive scalar in the fluid moving along the diagonal. Secondly, we look at a static diffusion: we turn on both hydrodynamics and radiation, but the fluid is at rest $\vec{u}=\vec{0}$. Finally, we conjugate the two precedent tests to look at dynamical diffusion: a Gaussian pulse is initialized in a moving fluid. This test will characterize the relative importance of advection and diffusion terms and will verify their coupling.

Figure E.5: Radiative energy for diffusion test: pure advection, static diffusion, diffusion in a moving fluid.

Figure E.5 shows the radiative energy in these three cases whereas Figs. E.6 and E.7 show slices with the pulse always recentered by a factor $\vec{ u}\Delta t$. In the first case, the Gaussian pulse should be advected without being affected at all. We first verify that the advection speed is correct because the two peaks are centered at the same position. Then we can see that the pulse is slightly diffused (the maximum has decreased by 8% after crossing over 12 cells). This diffusion could easily be lowered by using a more aggressive slope limiter. When we turn on hydrodynamics, all the pulses are again properly centered, which means that advection is still done at the right speed. Moreover, the diffusion of the radiative pulse in the moving and in the static fluids is very similar (the maximum discrepancy is less than 2%). Therefore, our treatment of the comoving terms preserves the good accuracy of the advection scheme and there is no problem in coupling advection with diffusion terms.

Figure E.6: Radiative energy for pure advection: initial pulse (dashed line) and pulse after crossing one eighth of the box (solid line). The pulse is recentered by a factor  $\vec{ u}\Delta t$.

Figure E.7: Radiative energy for: initial pulse (dotted line) and pulse at the end of the simulation with advection (dashed line) and without (solid line). Right panel is a zoom. All the pulses are recentered by a factor  $\vec{ u}\Delta t$.

E.4 Matter-radiation coupling

We now present a test of the matter-radiation coupling. If the radiation energy density is much higher than the gas one, the matter-radiation coupling can be solved analytically by assuming a constant radiative energy (Turner & Stone 2001). The matter energy then obeys the following equation where $E_{\rm r}$ is constant:

$$\partial_t e = - \sigma c (a_{\rm r} T^4 - E_{\rm r}).$$ (E.1)

The previous equation can be integrated by noticing that:

$$\begin{array}{ll} \int_{x_0}^{x_1}{\frac{{\rm d}x}{x^4-a^4}}=\fra... ...g(\frac{a-x_1}{a-x_0})-\frac{1}{4}\log(\frac{a+x_1}{a+x_0})\right). \end{array}$$ (E.2)

We have solved this matter-radiation coupling problem with HERACLES using various time steps. The time step was kept constant in units of $\tau(e)=\frac{e}{\sigma c (a_{\rm r} T^4 - E_{\rm r})}$. A time step equal to $\tau$ corresponds to a variation in e of about 100% at each time step.

Figure E.8: Gas energy versus time in the matter-radiation coupling test. The left (resp. right) subfigure corresponds to an initial matter energy of 102 (resp. 1010) erg cm-3. The thick solid line corresponds to the analytical solution, and squares (resp. diamonds, stars, triangles, crosses) to HERACLES simulation with a constant time step equal to 0.1 $\tau (e)$ (resp. 1, 10, 100, 105). Each point corresponds to a time step of the simulation.

We did two tests with a radiative energy $E_{\rm r} = 10^{12}$ erg cm^-3(i.e. $T_{\rm rad} = 3.39\times 10^6 \; {\rm K}$), an opacity $\sigma = 4 \times 10^{-8}$ cm^-1, and a heat capacity $c_v = 20.79 \; {\rm erg~cm^{-3} ~ K^{-1}}$. In the first one, the initial gas energy was 10² erg cm^-3 (i.e. $T_{\rm gas} = 4.81\;{\rm K} < T_{\rm rad}$), and in the second one 10¹⁰ erg cm^-3 (i.e. $T_{\rm gas} = 4.81\times 10^8\; {\rm K} > T_{\rm rad}$). The corresponding results are plotted in Fig. E.8. We can see that if the variation of the energy is well-sampled (i.e. the time step is less than $\tau$), then the analytical solution is properly reproduced. If the time step is large compared to the coupling time, the shape of the energy variation moves away from the analytical solution (as could be expected) but the proper value for the equilibrium is nevertheless found.

E.5 Marshak wave test

A Marshak wave is the propagation of a radiation heating front in a cold slab (Marshak 1958). An analytical solution of the Marshak wave problem was proposed by Su & Olson (1996) in the context of the diffusion equation and assuming that the gas heat capacity is proportional to T³ (i.e. $c_v=\alpha T^3$). We compared the solution obtained by HERACLES to their solution for the case $\epsilon \equiv 4 a_{\rm r}/\alpha=0.1$.

Figure E.9 represents the radiative and gas energies (normalized by their boundary values) obtained by HERACLES and by Su & Olson (1996) at times given by $\tau \equiv \epsilon c \sigma t = (0.1,1,10)$ (cf. their Fig. 3). At early times and at the foot of the Marshak wave, we can see large differences between the two methods. This was to be expected since in this region the radiation is instead in the transport regime and the diffusion approximation is then different from our M₁ model. At later times and inside the wave, the reduce flux is much smaller (i.e. the radiation is almost in the diffusion regime) and, accordingly, the results given by M₁ are in good agreement with the analytical solution.

Figure E.9: Marshak wave test: comparison of Heracles results (radiative energy with a solid line and gas energy with a dashed line) with the Su & Olson (1996) solution (radiative energy with stars and gas energy with crosses) at three different times.