© H. Bloch et al. 2021

1. Introduction

In many astrophysical situations, radiation is an important process which interacts with the surrounding gas; for example in (exo)planet atmospheres (e.g., Thomas & Stamnes 2002), massive stars (e.g., Kuiper et al. 2010; Mignon-Risse et al. 2020), and H II regions (e.g., Spitzer 1978) up to the cosmic reionization (e.g., Stiavelli 2009). In all these situations, the radiation can be absorbed, thus heating the surrounding gas. Photons could also be emitted by the gas or by another source. Photons could be scattered, which will change their direction of propagation and perhaps their frequency (Chandrasekhar 1960).

Two main regimes can arise, depending on the mean free path of photons compared to the characteristic length of the system (Mihalas & Mihalas 1984). On one hand, in the diffusive limit, the medium is optically thick (mean free path of photons much smaller than the characteristic length), and the radiation and the matter strongly interact with each other. On the other hand, in the free-streaming regime, the radiation does not affect the gas, and the medium is optically thin (mean free path of photons greater than the characteristic length). Numerically, it is difficult to accurately capture both limits.

Because of the high number of degrees of freedom (i.e., the time, the position, the direction of propagation, and the frequency of photons), only a few problems can be solved analytically (Chandrasekhar 1960) and direct simulations are out of reach for modern computers. Different models have been developed to reduce the computational cost. Here, we focus on the moment models (Levermore 1984), in which the specific intensity is averaged over the direction of propagation of photons. These models present several advantages: the computational cost is lower than other methods such as the Monte-Carlo method, and in general it is easy to couple these models with a grid-based hydrodynamics code.

One can only consider the moment of zero order (the radiative energy), leading to the flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981). Because the moment models consider only the moment of zero order, their computational cost is quite low, but they are very diffusive in the free-streaming regime. To tackle this issue, one can use a two-moment model (radiative energy and radiative flux), with the M₁ closure relation (Dubroca & Feugeas 1999). However, this method can suffer from artifacts when multiple beams cross in the free-streaming regime (González 2006). One can solve this issue by using a three-moment model (radiative energy, radiative flux, and radiative pressure) with the M₂ closure relation (Pichard et al. 2016). However, because of the increase of unknowns, the computational cost also increases. In this work, we have chosen to use the two-moment model with the M₁ closure relation because the computational cost remains affordable, and in our applications we do not encounter the problem of beams crossing in the free-streaming regime.

Even though the M₁ model is accurate in both free-streaming and diffusive regimes at the continuous level, numerical schemes also need to properly capture both limits. Several approaches have been developed. For example, Berthon & Turpault (2011) presented a scheme based on an HLL solver with source terms modified with a free parameter. Following this idea, we propose a new so-called asymptotic preserving scheme also based on an HLL solver. Nevertheless, we have chosen another parameter to recover the asymptotic behavior in the diffusive limit. Furthermore, our integration of source terms is different. In many physical applications (e.g., clouds), optically thick regions are found next to optically thin zones. We propose a well-balanced modification of the source term, which allows us to accurately reach steady states in the presence of sharp transitions.

As a first application, we are interested in the development of H II regions in massive pre-stellar dense cores (Churchwell 2002). We focus on solving the radiation hydrodynamics equations. An explicit solver for the radiative transfer would be restricted by a Courant-Friedrichs-Lewy (CFL) condition, limited by the speed of light. This will result in a very low time-step compared to the hydrodynamics time-step, which is limited by the speed of sound through the fluid. Several methods have been developed to get around this problem; we have chosen a time-implicit solver (e.g., González et al. 2007). The temporality of the radiative transfer is preserved, which is not the case with the reduced-speed-of-light approximation (RSLA, e.g., Gnedin & Abel 2001). However, a time-implicit solver method is costly because it requires solving large sparse linear systems.

Fortunately, progress has been made in linear algebra for high-performance computing (HPC). Because the linear system we have to solve is large and sparse, direct methods are out of reach. Iterative solvers with preconditioners have been developed to tackle this issue (e.g., Saad 2003). In this work, we use the library Trilinos (Heroux et al. 2005) because it allows us to target different architectures, such as multi-core, many-core, and GP-GPU. It also provides, among others, algebraic multigrid (AMG) preconditioners.

The paper is organized as follows. In the following section, we present the moment model and the M₁ closure relation in more detail. In Sect. 3, we outline our new numerical scheme, which is well balanced and asymptotic preserving in the diffusive limit. In Sect. 4, we provide details of some implementation features of our code ARK-RT, especially the Kokkos and Trilinos libraries used for shared memory parallelism and linear algebra for high-performance computing. We also show some performance results. In Sect. 5, we present some numerical test cases to show the importance of the asymptotic preserving and well-balanced properties. Finally, in Sect. 5.6 we present a physical application of our method where we study the stability of the ionization front in the dense cores of a H II region.

2. Model

2.1. Physical model

In this work, we only consider gray radiative transfer, that is, we have computed the average over the frequency. We also assume local thermodynamic equilibrium (LTE) and we do not consider scattering. The mathematical description of radiative transfer was formalized by Chandrasekhar (1960):

$$\begin{array}{c}\hfill (\frac{1}{c}\frac{\partial }{\partial t}+\mathbf{\Omega }·\mathrm{\nabla })I(\mathit{x},t,\mathbf{\Omega })=\sigma B(\mathit{x},t)-\sigma I(\mathit{x},t,\mathbf{\Omega }),\end{array}$$(1)

where I, the quantity of interest, is the specific intensity, c is the speed of light, σ is the opacity of the medium, and B is the black-body-specific intensity. Also, x, t, and Ω are the spatial, temporal, and angular variables. This model can be generalized to multigroup radiative transfer (e.g., Turpault 2005).

Because the specific intensity I depends on six variables in the three-dimensional case, the numerical treatment of Eq. (1) becomes rapidly costly. To reduce the computational cost, we use a moment method, which we present in detail in the following section.

2.2. M₁ model

Let us consider the three first moments of the specific intensity: the radiative energy E_r, the radiative flux F_r, and the radiative pressure P_r defined as:

$$\begin{array}{c}\hfill \begin{array}{c}{E}_r=\frac{1}{c}{\displaystyle {\oint }_{4\pi }}I(\mathit{x},t,\mathbf{\Omega })\mathrm{d}\mathbf{\Omega }\hfill \\ {\mathit{F}}_r={\displaystyle {\oint }_{4\pi }}I(\mathit{x},t,\mathbf{\Omega })\mathbf{\Omega }\mathrm{d}\mathbf{\Omega }\hfill \\ {\mathtt{P}}_r=\frac{1}{c}{\displaystyle {\oint }_{4\pi }}I(\mathit{x},t,\mathbf{\Omega })\mathbf{\Omega }\otimes \mathbf{\Omega }\mathrm{d}\mathbf{\Omega }.\hfill \end{array}\end{array}$$(2)

The mean over solid angles of Eq. (1) and its product by Ω give the following system:

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\partial }_t{E}_r+\mathrm{\nabla }·{\overrightarrow{F}}_r=\hfill & c\sigma (a_r{T}_g^4-E_r)(3\mathrm{a})\hfill \\ \hfill \\ {\partial }_t{\overrightarrow{F}}_r+c^2\mathrm{\nabla }·{\mathtt{P}}_r\hfill & =-c\sigma {\overrightarrow{F}}_r,(3\mathrm{b})\hfill \end{array}\end{array}$$

where T_g is the gas temperature and a_r is the radiation constant.

The fluid and the radiation exchange energy and momentum through emission and absorption. To ensure the conservation of the total energy when the hydrodynamics is frozen, the energy exchange term is given by

$$\begin{array}{c}\hfill {\partial }_t(\rho c_v{T}_{\mathrm{g}})=-c\sigma (a_r{T}_{\mathrm{g}}^4-E_r).\end{array}$$(4)

Here, ρc_vT_g is the gas internal energy, with ρ the density of the fluid and c_v the heat capacity, which is defined by $c_{\mathit{v}}=\frac{k_b}{\mu m_{\mathrm{H}}(\gamma -1)}$ for a perfect gas, where k_b is the Boltzmann constant, μ is the mean molecular weight, m_H is the mass of hydrogen, and γ is the adiabatic index of the gas.

We still have to specify a closure relation, that is, a way to express P_r as a function of E_r and F_r. The one chosen here is the M₁ model (Levermore 1984). From Dubroca & Feugeas (1999), we write P_r = DE_r where D is the Eddington tensor, defined by $\mathtt{D}=\frac{1-\chi }{2}\mathtt{I}+\frac{3\chi -1}{2}\mathit{n}\otimes \mathit{n}$ with χ the Eddington factor, I the identity matrix, and $\mathit{f}=\frac{{\mathit{F}}_r}{c{E}_r}=f\mathit{n}$ is the reduced flux. Here, χ can be specified either by applying a Lorentz transform to an isotropic distribution of photons (Levermore 1984), or by minimizing the radiative entropy (Dubroca & Feugeas 1999). In both cases, we have $\chi =\frac{3+4f^2}{5+2\sqrt{4-3f^2}}$. Let us note that f ≤ 1, which ensures that the radiative energy cannot be transported faster than the speed of light.

The M₁ model preserves both free-streaming and diffusive limits. On one hand, if f = 1, then P_r = DE_rn ⊗ n, and only the transport regime remains. On the other hand, if f = 0, the model in the diffusive regime simplifies to the P₁ model, with ${\mathtt{P}}_r=\frac{1}{3}{E}_r\mathtt{I}$. The radiative pressure tensor becomes isotropic. We look into the diffusion limit more precisely in Sect. 2.4.

2.3. Radiation hydrodynamics

We now consider the radiation hydrodynamics equations. The fluid evolution is described by the Euler equations expressing conservation of mass, balance of momentum, and balance of energy. Because photons are relativistic particles, we have to evaluate the quantities in the laboratory frame or the comoving frame, which is moving with the fluid. On one hand, using the comoving frame introduces nonconservative terms in the left-hand side of the equations. On the other hand, the hyperbolic part of the system remains simple in the laboratory frame but some source terms have to be incorporated to describe the interactions between the matter and the radiation. We have chosen the second approach. These source terms characterize the momentum and energy exchanges between the fluid and the radiation:

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill {\partial }_t\rho +\mathrm{\nabla }·(\rho \overrightarrow{u})& =0(5\mathrm{a})\hfill \\ \hfill {\partial }_t(\rho \overrightarrow{u})+\mathrm{\nabla }·(\rho \overrightarrow{u}\otimes \overrightarrow{u}+p\mathtt{I})& =\rho \overrightarrow{g}+\frac{\sigma }{c}{\overrightarrow{F}}_r-\frac{1}{c}{\overrightarrow{S}}_{{\overrightarrow{F}}_r}(\overrightarrow{u})(5\mathrm{b})\hfill \\ \hfill {\partial }_t(\rho E)+\mathrm{\nabla }·((\rho E+p)\overrightarrow{u})& =\rho \overrightarrow{g}·\overrightarrow{u}\hfill \\ \hfill & -c\sigma (a_r{T}_g^4-E_r)-S_{E_r}(\overrightarrow{u})(5\mathrm{c})\hfill \\ \hfill {\partial }_t{E}_r+\mathrm{\nabla }·{\overrightarrow{F}}_r& =c\sigma (a_r{T}_g^4-E_r)+S_{E_r}(\overrightarrow{u})(5\mathrm{d})\hfill \\ \hfill {\partial }_t{\overrightarrow{F}}_r+c^2\mathrm{\nabla }·{\mathtt{P}}_r& =-c\sigma {\overrightarrow{F}}_r+c{\overrightarrow{S}}_{{\overrightarrow{F}}_r}(\overrightarrow{u}),(5\mathrm{e})\hfill \end{array}\end{array}$$

where u is the material velocity, p is the pressure of the fluid, g is the external gravitational field, and $\rho E=\rho e+\frac{1}{2}\rho {\mathit{u}}^2$ is the density of total matter energy with e the specific internal energy. The terms S_Er(u) and S_Fr(u) depend on the velocity u. Using Eqs. (29) and (31) from Lowrie et al. (1999), we have

$$\begin{array}{c}\hfill \{\begin{array}{c}{S}_{E_r}(\mathit{u})=\frac{\sigma }{c}\mathit{u}·{\mathit{F}}_r+\frac{\sigma }{c}{E}_r\mathit{u}·\mathit{u}+\frac{\sigma }{c}\mathit{u}·(\mathit{u}·{\mathtt{P}}_r)\hfill \\ {\mathit{S}}_{{\mathit{F}}_r}(\mathit{u})=\sigma \mathit{u}·{\mathtt{P}}_r+\sigma a_r{T}_{\mathrm{g}}^4\mathit{u}+\frac{\sigma }{c^2}\mathit{u}·(\mathit{u}·{\mathit{F}}_r),\hfill \end{array}\end{array}$$(6)

at first order in $\frac{\mathit{u}}{c}$. To close the system, we also add the equation of state of an ideal gas: p = ρe(γ − 1).

2.4. Diffusive limit in a static fluid

Let us now focus on the diffusive regime with the hydrodynamics frozen, that is, the limit of large opacity, long timescale, and u = 0. We consider the P₁ closure relation

$$\begin{array}{c}\hfill \{\begin{array}{c}{\partial }_t{E}_r+\mathrm{\nabla }·{\mathit{F}}_r=c\sigma (a_r{T}_{\mathrm{g}}^4-E_r)\hfill \\ {\partial }_t{\mathit{F}}_r+\frac{c^2}{3}\mathrm{\nabla }{E}_r=-c\sigma {\mathit{F}}_r\hfill \\ {\partial }_t(\rho c_v{T}_{\mathrm{g}})=-c\sigma (a_r{T}_{\mathrm{g}}^4-E_r).\hfill \end{array}\end{array}$$(7)

Following Berthon & Turpault (2011), we introduce a rescaling parameter ε to write the time and opacity as t̃ = εt and $\stackrel{\sim }{\sigma }=\epsilon \sigma$, respectively. The radiative energy, radiative flux, and gas temperature are expanded with ε; for example E_r = E_r, 0 + εE_r, 1 + 𝒪(ε²). System (7) becomes

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\epsilon }^2{\partial }_{\stackrel{\sim }{t}}{E}_r+\epsilon \mathrm{\nabla }·{\overrightarrow{F}}_r\hfill & =c\stackrel{\sim }{\sigma }(a_r{T}_g^4-E_r)(8\mathrm{a})\hfill \\ {\epsilon }^2{\partial }_{\stackrel{\sim }{t}}{\overrightarrow{F}}_r+\epsilon \frac{c^2}{3}\mathrm{\nabla }{E}_r\hfill & =-c\stackrel{\sim }{\sigma }{\overrightarrow{F}}_r(8\mathrm{b})\hfill \\ {\epsilon }^2{\partial }_{\stackrel{\sim }{t}}(\rho c_v{T}_g)\hfill & =-c\stackrel{\sim }{\sigma }(a_r{T}_g^4-E_r).(8\mathrm{c})\hfill \end{array}\end{array}$$

By expanding Eqs. 8a and 8b at zero order, we have

$$\begin{array}{c}\hfill \{\begin{array}{c}{E}_{r,0}=a_r{T}_{\mathrm{g},0}^4\hfill \\ {\mathit{F}}_{r,0}=\mathbf{0}.\hfill \end{array}\end{array}$$(9)

Expanding Eq. (8b) at first order leads to

$$\begin{array}{c}\hfill {\mathit{F}}_{r,1}=-\frac{c}{3\stackrel{\sim }{\sigma }}\mathrm{\nabla }{E}_{r,0}.\end{array}$$(10)

Finally, expanding the sum of Eqs. (8a) and (8c) at second order gives

$$\begin{array}{c}\hfill {\partial }_{\stackrel{\sim }{t}}(E_{r,0}+\rho c_v{T}_{\mathrm{g},0})-\mathrm{\nabla }\left(\frac{c}{3\stackrel{\sim }{\sigma }}\mathrm{\nabla }{E}_{r,0}\right)=0.\end{array}$$(11)

In the diffusive limit, the total energy E_r + ρc_vT_g at zero order obeys the diffusion equation given by Eq. (11). Similar development of the radiation hydrodynamics case is presented in Appendix A. Section 3 presents our design of an asymptotic preserving scheme, that is, a numerical scheme that will degenerate to the discretization of Eq. (11) in the diffusive regime.

3. Numerical scheme and algorithm

3.1. Radiation transport in a static fluid

Let us first introduce some notations: we use Δx to denote the step along the x-direction. here, Δt is the time interval between the current time tⁿ and t^n + 1. We write x_i the center of the cell i and $x_{i+\frac{1}{2}}$ the interface between the cell i and the cell i + 1. We use the notation $u_i^n$ to represent the averaged quantity associated with the field u at time tⁿ in the cell i (finite volume). Finally, we use $u_{i+\frac{1}{2}}^n$ to represent the quantity associated with the field u at time tⁿ and at the interface between cells i and i + 1.

The development of the numerical scheme is presented only in the one-dimensional case, but its extension to higher dimensions is straightforward. To ease notations, we drop the indices r for all radiative variables.

We present the development of the numerical scheme using a time-implicit integration, but a similar development can be done with a semi-implicit solver: source terms remain implicit, but the hyperbolic part is time-explicit.

3.1.1. Hyperbolic system

Following González et al. (2007), we discretize the hyperbolic part of Eq. (3) using a first-order Godunov-type solver (Toro 2009). From Berthon & Turpault (2011), we also introduce an extra parameter α which is specified in Sect. 3.1.3:

$$\begin{array}{c}\hfill \{\begin{array}{c}{E}_i^{n+1}=E_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({\alpha }_{i+\frac{1}{2}}{\mathcal{F}}_{i+\frac{1}{2}}^{\ast }-{\alpha }_{i-\frac{1}{2}}{\mathcal{F}}_{i-\frac{1}{2}}^{\ast })\hfill \\ +c{\sigma }_i\mathrm{\Delta }t(a_r{\left(T_i^{n+1}\right)}^4-E_i^{n+1})\hfill \\ F_i^{n+1}=F_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({\mathcal{P}}_{i+\frac{1}{2}}^{\ast }-{\mathcal{P}}_{i-\frac{1}{2}}^{\ast })-c\mathrm{\Delta }t{\{\sigma F\}}_i^{n+1}\hfill \\ \rho c_v{T}_i^{n+1}=\rho c_v{T}_i^n-c{\sigma }_i\mathrm{\Delta }t(a_r{\left(T_i^{n+1}\right)}^4-E_i^{n+1}),\hfill \end{array}\end{array}$$(12)

where ${\mathcal{F}}_{i+\frac{1}{2}}^{\ast }$ and ${\mathcal{P}}_{i+\frac{1}{2}}^{\ast }$ are the numerical fluxes given by

$$\begin{array}{c}\hfill \begin{array}{c}{\mathcal{F}}_{i+\frac{1}{2}}^{\ast }=\frac{{\lambda }_{i+\frac{1}{2}}^{+}{F}_i^{n+1}-{\lambda }_{i+\frac{1}{2}}^{-}{F}_{i+1}^{n+1}}{{\lambda }_{i+\frac{1}{2}}^{+}-{\lambda }_{i+\frac{1}{2}}^{-}}+\frac{{\lambda }_{i+\frac{1}{2}}^{+}{\lambda }_{i+\frac{1}{2}}^{-}}{{\lambda }_{i+\frac{1}{2}}^{+}-{\lambda }_{i+\frac{1}{2}}^{-}}(E_{i+1}^{n+1}-E_i^{n+1})\hfill \\ {\mathcal{P}}_{i+\frac{1}{2}}^{\ast }=c^2\frac{{\lambda }_{i+\frac{1}{2}}^{+}{P}_i^{n+1}-{\lambda }_{i+\frac{1}{2}}^{-}{P}_{i+1}^{n+1}}{{\lambda }_{i+\frac{1}{2}}^{+}-{\lambda }_{i+\frac{1}{2}}^{-}}+\frac{{\lambda }_{i+\frac{1}{2}}^{+}{\lambda }_{i+\frac{1}{2}}^{-}}{{\lambda }_{i+\frac{1}{2}}^{+}-{\lambda }_{i+\frac{1}{2}}^{-}}(F_{i+1}^{n+1}-F_i^{n+1}),\hfill \end{array}\end{array}$$(13)

with ${\lambda }_{i+\frac{1}{2}}^{+}=max(0,{\lambda }_{max})$ and ${\lambda }_{i+\frac{1}{2}}^{-}=min(0,{\lambda }_{min})$, where λ_max and λ_min are the eigenvalues of Eq. (3). From Berthon et al. (2007), we have

$$\begin{array}{c}\hfill {\lambda }_{max,min}=c(\frac{f_x}{\xi }\pm \frac{\sqrt{2}\sqrt{(\xi -1)(\xi +2)(2(\xi -1)(\xi +2)+3f_y^2)}}{\sqrt{3}\xi (\xi +2)}),\end{array}$$(14)

with $\xi =\sqrt{4-3f^2}$. See Fig. 1 of González et al. (2007) for more details about the structure of the eigenvalues. $\{\sigma F{\}}_i^{n+1}$ is a well-chosen discretization of the term σF_r in the cell i and at time t^n + 1, and is specified in the following section.

3.1.2. Well-balanced modification of the source term

From Berthon et al. (2015), a well-balanced scheme catches the correct steady regime. The steady state, if it exists, is given by

$$\begin{array}{c}\hfill \{\begin{array}{cc}{E}_r\hfill & =a_r{T}_g^4(15\mathrm{a})\hfill \\ \mathrm{\nabla }·{\overrightarrow{F}}_r\hfill & =0(15\mathrm{b})\hfill \\ c\mathrm{\nabla }·{\mathtt{P}}_r\hfill & =-\sigma {\overrightarrow{F}}_r(15\mathrm{c}).\hfill \end{array}\end{array}$$

Equation (15c) is discretized by $c\frac{P_{i+1}^{n+1}-P_{i-1}^{n+1}}{2\mathrm{\Delta }x}=-{\{\sigma F\}}_i^{n+1}$. An obvious choice for $\{\sigma F{\}}_i^n$ is

$$\begin{array}{c}\hfill {\{\sigma F\}}_i^{n+1}={\sigma }_i{F}_i^{n+1}.\end{array}$$(16)

However, using this formulation, Eq. (15c) is discretized as

$$\begin{array}{c}\hfill -\frac{c}{2}({(\mathrm{\nabla }·P)}_{i+\frac{1}{2}}^{n+1}+{(\mathrm{\nabla }·P)}_{i-\frac{1}{2}}^{n+1})={\sigma }_i{F}_i,\end{array}$$(17)

with ${(\mathrm{\nabla }·P)}_{i+\frac{1}{2}}^{n+1}=\frac{P_{i+1}^{n+1}-P_i^{n+1}}{\mathrm{\Delta }x}$. The radiative flux remains cell-centered and is equal to the divergence of radiative pressure, which itself is defined at the interfaces of the cells. This can create some spurious flux at the interface when looking for a steady state with a constant flux in the box (see Sect. 5.2). Inspired by well-balanced schemes for hydrodynamics (e.g., Padioleau et al. 2019) which preserve the hydrostatic balance between the pressure forces and the gravitational force (and the similarity between this balance and that between radiative pressure and the radiative flux source term in Eq. (15c)), we choose to use an average of a face discretization of the radiative flux source term:

$$\begin{array}{c}\hfill {\{\sigma F\}}_i^{n+1}=\frac{1}{2}({\sigma }_{i+\frac{1}{2}}{F}_{i+\frac{1}{2}}^{n+1}+{\sigma }_{i-\frac{1}{2}}{F}_{i-\frac{1}{2}}^{n+1}),\end{array}$$(18)

with

$$\begin{array}{c}\hfill \{\begin{array}{c}{\sigma }_{i+\frac{1}{2}}=\frac{1}{2}({\sigma }_i+{\sigma }_{i+1})\hfill \\ F_{i+\frac{1}{2}}^{n+1}=\frac{1}{2}(F_i^{n+1}+F_{i+1}^{n+1}).\hfill \end{array}\end{array}$$(19)

One way to interpret this equation is to remember that

$$\begin{array}{c}\hfill {\{\sigma F\}}_i^{n+1}=\frac{1}{\mathrm{\Delta }x}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}}\sigma (x){\mathit{F}}_r(t^{n+1},x)\mathrm{d}x.\end{array}$$(20)

Equation (16) is obtained with the rectangle rule for numerical integration of Eq. (20):

$$\begin{array}{cc}& {\{\sigma F\}}_i^{n+1}=\frac{x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}}{\mathrm{\Delta }x}\left(\sigma \left(\frac{x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}}{2}\right){\mathit{F}}_r(t^{n+1},\frac{x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}}{2})\right)\hfill \\ \hfill & =\sigma (x_i){\mathit{F}}_r(t^{n+1},x_i)\hfill \\ \hfill & ={\sigma }_i{F}_i^{n+1},\hfill \end{array}$$(21)

whereas Eq. (18) is given by the trapezoidal rule:

$$\begin{array}{cc}& {\{\sigma F\}}_i^{n+1}=\frac{x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}}{2\mathrm{\Delta }x}(\sigma \left(x_{i-\frac{1}{2}}\right){\mathit{F}}_r(t^{n+1},x_{i-\frac{1}{2}})\hfill \\ \hfill & +\sigma \left(x_{i+\frac{1}{2}}\right){\mathit{F}}_r(t^{n+1},x_{i+\frac{1}{2}}))\hfill \\ \hfill & =\frac{1}{2}({\sigma }_{i-\frac{1}{2}}{F}_{i-\frac{1}{2}}^{n+1}+{\sigma }_{i+\frac{1}{2}}{F}_{i+\frac{1}{2}}^{n+1}).\hfill \end{array}$$(22)

To have

$$\begin{array}{c}\hfill {\sigma }_{i+\frac{1}{2}}{F}_{i+\frac{1}{2}}^{n+1}=-c\frac{P_{i+1}^{n+1}-P_i^{n+1}}{\mathrm{\Delta }x}\end{array}$$(23)

in the whole domain, we also impose it as boundary condition:

$$\begin{array}{c}\hfill {\sigma }_{\frac{1}{2}}{F}_{\frac{1}{2}}^{n+1}=-c\frac{P_1^{n+1}-P_0^{n+1}}{\mathrm{\Delta }x},\end{array}$$(24)

where $P_0^{n+1}$ is the radiative pressure given by the boundary condition. In that way, the radiative flux is centered at the interfaces of the cells, as is the divergence of radiative pressure. A von Neumann stability analysis of the modified scheme is presented in Appendix B.

3.1.3. Asymptotic preserving scheme

We still have to specify our choice for ${\alpha }_{i+\frac{1}{2}}$. Here, ${\alpha }_{i+\frac{1}{2}}=1$ corresponds to a classic HLL scheme. However, the solution given by an asymptotic preserving scheme has to approximate the solution of Eq. (11) as soon as the asymptotic regime is reached, that is, it must incorporate large opacities and long timescales. Unfortunately, a standard HLL scheme does not have this property (see Sect. 5.1). To tackle this issue and obtain an asymptotic preserving scheme, we choose

$$\begin{array}{c}\hfill {\alpha }_{i+\frac{1}{2}}=\frac{1}{1-3{\sigma }_{i+\frac{1}{2}}\mathrm{\Delta }x(1-f_{i+\frac{1}{2}}^2)\frac{{\lambda }_{i+\frac{1}{2}}^{+}{\lambda }_{i+\frac{1}{2}}^{-}}{c({\lambda }_{i+\frac{1}{2}}^{+}-{\lambda }_{i+\frac{1}{2}}^{-})}},\end{array}$$(25)

with $f_{i+\frac{1}{2}}=\frac{1}{2}(f_i^n+f_{i+1}^n)$. The derivation of Eq. (25) is presented in Appendix C. Other choices can be made (see e.g., Berthon & Turpault 2011). If ${\sigma }_{i+\frac{1}{2}}\mathrm{\Delta }x$ goes to zero, ${\alpha }_{i+\frac{1}{2}}$ goes to 1, and we recover a standard HLL scheme. Considering the diffusive limit, we prove that the scheme is asymptotic preserving in Appendix C. We show that

$$\begin{array}{c}\hfill \{\begin{array}{c}{E}_{i,0}^{n+1}=a_r{\left(T_{i,0}^{n+1}\right)}^4\hfill \\ F_{i,0}^{n+1}=0\hfill \\ {\sigma }_{i+\frac{1}{2}}{F}_{i+\frac{1}{2},1}^{n+1}=-\frac{c}{3\mathrm{\Delta }x}(E_{i+1,0}^{n+1}-E_{i,0}^{n+1})\hfill \\ E_{i,0}^{n+1}+\rho c_v{T}_{i,0}^{n+1}=E_{i,0}^n+\rho c_v{T}_{i,0}^n\hfill \\ +\frac{c\mathrm{\Delta }t}{3\mathrm{\Delta }{x}^2}(\frac{E_{i+1,0}^{n+1}-E_{i,0}^{n+1}}{{\sigma }_{i+\frac{1}{2}}}-\frac{E_{i,0}^{n+1}-E_{i-1,0}^{n+1}}{{\sigma }_{i-\frac{1}{2}}}),\hfill \end{array}\end{array}$$(26)

which is a standard discretization of Eqs. (9)–(11).

Unfortunately, we cannot prove that this scheme will preserve the admissible states f <  1 and, indeed, numerical experiments with this scheme have shown that we can get f >  1 when we are close to the free-streaming regime. In these situations, we can either enforce f <  1 (Sect. 5.6) or come back to a centered discretization of the source term (Sect. 5.4). Furthermore, the development of the asymptotic preserving scheme with the well-balanced modification of the source term is only done in the case of a static fluid. Using the asymptotic correction, Equation (25) is only the first step to obtaining an asymptotic preserving scheme in the case of a moving fluid (see Sect. 6.2).

3.2. Coupling to hydrodynamics

Following González et al. (2007), the resolution of the whole of the system shown in Eq. (5) describing radiation hydrodynamics is split into three steps:

1. Update the hydrodynamics quantities (Eqs. (5a)–(5c) without the terms of energy and momentum exchange) using the well-balanced and all-regime solver developed in Padioleau et al. (2019);

2. Update the radiative quantities and gas temperature (Eqs. 3 and (4)) using the solver developed in Sect. 3.1. During this step, the hydrodynamics quantities are frozen;

3. Add source terms S_Er(u) and S_Fr(u). For simplicity, all source terms which depend on the velocity are treated explicitly. The term $\frac{\sigma }{c}{\mathit{F}}_r$ in Eqs. (5b) and (5e) is discretized using the well-balanced scheme proposed in Sect. 3.1.2. All the other terms remain cell-centered.

This splitting allows us to reduce the number of equations solved implicitly, making the method more efficient.

3.3. Algorithm for nonlinear implicit solver

The time-step given by the CFL condition is much smaller for radiation than for hydrodynamics. Indeed, for the radiation, it is limited by the speed of light, whereas it is limited by the speed of sound in the fluid for the hydrodynamics. Because we turn to radiation hydrodynamics, we consider a long timescale for the radiative transfer. Therefore, we use a time-implicit integration for the radiative transfer.

3.3.1. Newton–Raphson method and linear solver

Because of the Eddington tensor, the eigenvalues in the numerical fluxes, and the $a_r{T}_{\mathrm{g}}^4$ factor, the system is nonlinear and is solved using a Newton-Raphson method. At each iteration, we have to solve a linear system. Because the system is large ((2 + d)N unknowns, where d is the number of dimensions and N the total number of cells) and sparse, it cannot be solved using a direct method. Because of the numerical fluxes, the matrix is not symmetric, and we use the biconjugate gradient stabilized method (Van der Vorst 1992).

3.3.2. Preconditioner

Using large time-steps for the radiative transfer, the matrix is ill-conditioned and iterative methods might not converge. One way to deal with this issue is to use a preconditioner. Instead of solving the original linear system Ax = b, we solve the right preconditioned system AK⁻¹Kx = b via solving AK⁻¹y = b to compute y and then Kx = y. As long as the matrix K is invertible, this gives the same solution as the original system. If K is well chosen, the condition number of the matrix AK⁻¹ is lower than that of A. For example, the Jacobi or diagonal preconditioner is given by

$$\begin{array}{c}\hfill {\mathtt{K}}_{\mathit{ij}}=\{\begin{array}{cc}{\mathtt{A}}_{\mathit{ij}}\hfill & \mathrm{if}i=j\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\end{array}$$(27)

See Saad (2003) for more details.

Numerous preconditioners have been developed over the years; some of them perform well and some do not, depending on the problem considered. Among different preconditioners, we explore the algebraic multigrid (AMG) technique. We present this method in the following section, and we compare it with other preconditioners in Sect. 4.4: a standard ILU(k) factorization (Saad 2003), a slightly modified variant of the standard ILU factorization (Saad 1994), an additive Schwarz domain decomposition (Cai & Sarkis 1999), and a classical relaxation method, Jacobi with damping (Saad 2003).

3.3.3. Algebraic multigrid preconditioner

The AMG methods were first developed as linear solvers for symmetric positive definite matrices arising from the discretization of scalar elliptic PDEs. For such a matrix, classical iterative methods are able to efficiently compute the high frequencies of the solution, but lack efficiency when computing its low frequencies. However, the computation is easier on a coarser grid with fewer unknowns. The idea of the multigrid solver is to build a coarser grid, solve the problem on this coarse grid, and then finally interpolate the solution on the fine grid. We can subsequently define a restriction operator R which transfers vectors from the fine grid to the coarse grid and an interpolation operator P used to return to the finer grid. P and R are nonsquared matrices. From Saad (2003), the main steps of the method are as follows:

1. Pre-smoothing: a few iterations of a simple method such as Jacobi or an incomplete factorization are performed to get the value $\stackrel{\mathbf{\sim }}{\mathit{x}}$.

2. The residual $\stackrel{\mathbf{\sim }}{\mathit{r}}=\mathit{b}-\mathtt{A}\stackrel{\mathbf{\sim }}{\mathit{x}}$ is projected over the coarse grid with the restriction operator R, to get the residual equation $\mathtt{RAP}\mathit{y}=\mathtt{R}\stackrel{\mathbf{\sim }}{\mathit{r}}$.

3. This equation is solved, possibly with a direct solver.

4. The solution y is interpolated over the fine grid with the interpolation operator P and then $\overline{\mathit{x}}=\stackrel{\mathbf{\sim }}{\mathit{x}}+\mathtt{P}\mathit{y}$.

5. Post-smoothing: a few iterations of a simple method are again performed to get the solution $\stackrel{\mathbf{\sim }}{\overline{\mathit{x}}}$.

The solution $\stackrel{\mathbf{\sim }}{\overline{\mathit{x}}}$ is used as a preconditioner result. If the coarse grid has too many unknowns to be solved directly, this process is applied recursively: the coarse grid becomes the fine grid and a coarser grid is built. Therefore, we have a hierarchy of grids. With a geometric multigrid solver, the restriction and interpolation operators are determined by the mesh, whereas with an algebraic multigrid solver they are automatically generated using data from the matrix.

Before describing our numerical tests, we present details of our implementation and parallelization in the following section.

4. Implementation and parallelization

We used the code ARK-RT¹ for implementation, which is a fork of the code ARK developed in Padioleau et al. (2019). The hydrodynamics and gravity part of our solver is similar to ARK and is solved with a well-balanced and all-regime solver. Because the solver for the radiative transfer equations is time-implicit, we have to solve large, sparse linear systems. This is done using the library Trilinos, described in the following section.

4.1. Linear algebra

We use the second generation of packages of the Trilinos framework (Heroux et al. 2005), namely Kokkos (Edwards et al. 2014) for shared memory computation, Tpetra (Baker & Heroux 2012) for distributed vectors and matrices, Belos (Bavier et al. 2012) for linear solvers, Ifpack2 (Prokopenko et al. 2016) for classical preconditioners, and MueLu (Berger-Vergiat et al. 2019a,b) for the AMG preconditioner. Let us first focus on the package Kokkos.

4.2. Shared memory computation

As new architectures have increasing numbers of cores, the distributed memory model is insufficient to take advantage of all the computational power available. Therefore, we need to use a shared memory model inside the nodes. Furthermore, computational nodes are increasingly heterogeneous; for example, multi-core, many-core, or accelerators such as GP-GPUs. Each architecture requires its own interface, such as OpenMP or C++11 threads for multi-core and many-core processors and CUDA or OpenACC for NVIDIA GPUs. This raises the problem of portability and performance portability: many HPC codes are optimized for some specific architectures, and so running the code on a different architecture will result in poor performance.

The library Kokkos (Edwards et al. 2014) tackles this issue. The user has a unique code which can be compiled with different shared memory models such as OpenMP or CUDA. Kokkos is based on different abstractions, like execution spaces (where a function is executed), memory spaces (where the data are stored), and execution policies (how the function is executed). Also, Kokkos provides execution patterns such as parallel loops and multidimensional arrays, and the storage is optimized according to the architecture.

4.3. Implementation

For the hydrodynamics step, Kokkos is used as an independent library for shared memory computation. Communications between the nodes are handled by the Message Passing Interface (MPI) programming model through a regular domain decomposition. Following Kestener (2017), inside each node, the domains are endowed with ghost cells used to implement physical boundary conditions, but also to contain values from neighbor domains. The code is organized with computational kernels, and each kernel is a C++ functor; see Padioleau et al. (2019) for more details.

The second step is the time-implicit solver for radiative transfer. The values of the matrix and the right-hand side of the linear system have to be updated at each iteration of the Newton–Raphson method, but the graph of the matrix does not change. Using Trilinos, the rows of the matrix are distributed across the MPI processes, each of them is associated with a unique global index. Each global index has a matching local index on the owning process.

Trilinos provides several methods to update the coefficients of the matrix. One of them uses only global indices, which is recommended for Trilinos. However, this function can only be called by the host. This has two main consequences: First, because all the rows of the matrix have to be updated, we use a sequential loop over the rows of the matrix. Second, in the case where we are using GP-GPUs, we have to update the matrix with data coming from the device. Therefore, we have to transfer some data from the device to the host, update the matrix with these data, and then transfer the matrix from the host to the device; this last step is done implicitly by Trilinos. This process increases the computational cost.

Another way is to use local indices. The package KokkosKernels, which is part of the Kokkos ecosystem, provides several ways to update the coefficients of the matrix through a kernel. This allows the use of a parallel loop (via Kokkos::parallel_for) and we avoid data transfers between the host and the device. For performance reasons, we use local indices.

4.4. Performances

Thanks to Trilinos, we can use many preconditioners. Unfortunately, they do not behave in the same way when the size of the system increases. All tests are performed on Poincare, our local cluster at Maison de la Simulation. Each node consists of two Sandy Bridge E5-2670 at 2.60 GHz (2 × 8 cores, 32 Go RAM) processors. We use a hybrid configuration MPI/OpenMP, with one MPI process per socket to avoid NUMA effects.

We first performed a weak scaling test, where we consider a two-dimensional case with periodic boundary conditions and a hot source located at the center of each domain. Each MPI process provides a piece of the whole domain of 1500² cells, and therefore the size of the system increases with the number of MPI tasks. The resolution is close to our target for three-dimensional simulations. As shown in the left panel of Fig. 1, for all preconditioners, the number of iterations remains constant, at around 10 iterations for the AMG preconditioner, around 20 iterations for both incomplete factorizations and the additive Schwarz domain decomposition, and around 250 iterations for the relaxation. As shown in the right panel, the speed-up reaches a plateau of 80% to 90% of maximum performance depending on the preconditioner.

Fig. 1. Weak scaling test. Each MPI process treats 15002 cells. We tested different preconditioners: Jacobi with damping (Relaxation), AMG, standard ILU(k) factorization (RILUK), a variant of the standard ILU factorization (ILUT), and additive Schwarz domain decomposition (Schwarz). Left panel: number of iterations to solve the linear system as a function of the number of cells. Right panel: speed-up as a function of the number of MPI processes.

Figure 2 shows the number of iterations and the speed-up as a function of the number of MPI processes for a strong scaling test. We now consider a Marshak wave propagation in the diffusive limit. The global resolution remains constant as the number of processes increases, and is set to 2048² cells. Because the global resolution is constant, one can expect the number of iterations to also remain constant when the number of MPI processes increases. However, using the AMG (orange curve) and the incomplete factorizations (RILUK, ILUT; green curve in left panel), when four MPI processes or more are used, the number of iterations is twice the number of iterations reached with one or two MPI processes. Therefore, the computational time is the same when using two or four MPI processes. Furthermore, all tested preconditioners and the linear solver require several communications per iteration, which likely becomes the main cost when the local resolution decreases.

Fig. 2. Strong scaling test. The global resolution is 20482 cells. The same preconditioners as those listed in Fig. 1 were tested. Left panel: number of iterations to solve the linear system as a function of the number of MPI processes. Right panel: speed-up as a function of the number of MPI processes.

Thanks to Kokkos, we can use exactly the same code on different architectures, namely Sandy Bridge processors and NVIDIA GP-GPUs (e.g., K80). Unfortunately, the memory required by the AMG preconditioner with a 1500² simulation is greater than the memory available on a K80 GPU. For the next tests, we use a lower resolution of 1000² cells. Table 1 summarizes the computational time for a fixed problem with different schemes and different architectures. As the explicit solver is restricted by a CFL condition, it requires several thousand time-steps whereas the implicit solver only needs a few time-steps to reach the same final time. Updating the matrix in parallel allows for a 25% reduction in computational time required. On a CPU, the implicit solver is around 160 times faster than the explicit solver, whereas on a GPU it is only 11 times faster.

Figure 3 compares the computational time with different preconditioners, on both a CPU and a GPU. Except for the implicit solver using the AMG preconditioner, all solvers are faster on a GPU than on a CPU; up to three times faster for the relaxation preconditioner. Part of the AMG algorithm probably remains sequential. On a CPU, the AMG preconditioner is faster than the relaxation preconditioner. The other preconditioners are slower, by up to a factor eight between the relaxation and the additive Schwarz domain decomposition on a GPU.

Fig. 3. Computational time for the implicit solver with different preconditioners (as listed for Fig. 1) on different architectures (Sandy Bridge CPU and K80 NVIDIA GPU). The resolution is 10002 cells.

Figure 4 compares the memory consumption with different preconditioners. Using GP-GPU, most of the data are located on the device, but Trilinos still allocates some memory, namely between 0.125 GB and 0.208 GB, on the host, unlike the explicit solver. Using the relaxation as a preconditioner, the amount of memory allocated is lower than with the other preconditioners (7.3 GB for the relaxation against 11.5 GB for the AMG).

Fig. 4. Memory consumption for the implicit solver with different preconditioners (as listed in Fig. 1) on different architectures (Sandy Bridge CPU and K80 NVIDIA GPU). The resolution is 10002 cells.

Choosing a well-suited preconditioner can be challenging and problem dependent. Once the preconditioner is chosen, it depends on many parameters. For example, Trilinos allows the user to choose the damping factor ω for the relaxation method or the smoother and the coarse solver for the AMG. Performances and stability can largely depend on these choices. For example, the relaxation method seems to be well suited for this problem with low computational time, and memory consumption, but in many other test cases, the linear solver will not converge. The AMG preconditioner performs well on a CPU but is less efficient on a GPU. Both incomplete factorizations and the additive Schwarz domain decomposition are slightly less efficient than the AMG preconditioner. Overall, we find the AMG preconditioner and the relaxation method to produce a good compromise between stability and performance.

The performances we obtain thanks to Kokkos and Trilinos are encouraging for the study of astrophysical problems. The time-step given by the hydrodynamics can be written as $\mathrm{CFL}\frac{\mathrm{\Delta }x}{c}$. Using a relaxation as a preconditioner, we need CFL ≥ 50 on a CPU and CFL ≥ 100 on a GPU to save computational time, whereas when using an incomplete factorization, we need CFL ≥ 250 on a CPU and CFL ≥ 1000 on a GPU. We need a larger CFL number on a GPU because the explicit solver is more efficient on a GPU than on a CPU.

In the following section, we use several numerical tests to show that the scheme developed in Sect. 3 is well suited for the study of radiation hydrodynamics problems.

5. Numerical results

We performed a series of verification tests to validate different properties of the scheme: the asymptotic correction with a Marshak wave, the well-balanced property to reach a steady state with a jump of opacity, the properties of the M₁ model with a beam test and a shadow test, and the coupling to the hydrodynamics with radiative shocks. We also present a physical application, investigating the stability of the ionization front in the dense core of a H II region. To ease notations, we define the radiative temperature as $T_r={\left(\frac{E_r}{a_r}\right)}^{\frac{1}{4}}$.

5.1. Marshak wave

From Mihalas & Mihalas (1984), a Marshak wave is the propagation of hot radiation into a cold medium. We consider a one-dimensional case in the diffusive limit in order to test the asymptotic preserving scheme developed in Sect. 3.1.3.

The length of the computational domain is 1 cm; it is discretized with 400 points. Initially, the medium is at equilibrium with the radiation: T₀ = T_r = 300 K, the initial radiative flux is F_r = 0. We consider a perfect gas with $\gamma =\frac{5}{3}$. The hydrodynamics is frozen. The density is constant, such that ρc_v = 1 J K⁻¹ cm⁻³, and the opacity is also constant, with σ = 10 000 cm⁻¹, and therefore σΔx = 25. At time t = 0, a source is lit at the left boundary with T_r = 1000 K.

The results are shown in Fig. 5 at time t_f = 2 × 10⁻⁴ s. We compare different solutions: a reference solution, the solution given by our asymptotic preserving scheme, and the solution given by a standard scheme. The reference solution is given by a standard discretization of Eq. (11). From Audit et al. (2002), when a scheme that is not asymptotic preserving is used with σΔx ≫ 1, the grid does not sample the mean free path of the photons and the solution is dominated by numerical diffusion. The relative L² error between the reference solution and the solution with ${\alpha }_{i+\frac{1}{2}}$ given by Eq. (25) is 1.1%, whereas with the standard HLL scheme the relative L² error is 84%. Using the asymptotic correction, we recover the correct behavior in the asymptotic limit.

Fig. 5. Marshak wave simulation. This figure shows a snapshot of the gas temperature at time tf = 2 × 10−4 s, with and without the asymptotic correction and the reference solution. Spatial resolution is n = 400 and the opacity is σ = 10 000 cm−1.

5.2. Steady state with a jump of opacity

In the previous case, the opacity is constant, we now consider a test with a jump of opacity, still in the one-dimensional case. We use this test to highlight the need for a well-balanced modification of the source term.

The length of the computational domain is 1 cm; it is discretized with 100 points. Initially, the medium is at equilibrium with the radiation: T₀ = T_r = 300 K, the initial radiative flux is F_r = 0. The opacity σ is now a function of space:

$$\begin{array}{c}\hfill \sigma (x)=\{\begin{array}{ccc}10000{\mathrm{cm}}^{-1}\hfill & \mathrm{if}\hfill & x<0.5,\hfill \\ 0\hfill & \mathrm{if}\hfill & x>0.5.\hfill \end{array}\end{array}$$(28)

At time t = 0, a source is lit at the left boundary with T_r = 1000 K.

Figure 6 shows the radiative flux at time t_f = 10⁻³ s. From Eq. (15b), when the steady state is reached, we expect the radiative flux to be constant in the box. Using a standard discretization of the source term, such as Eq. (16), a spurious peak located at the discontinuity of opacity is observed (orange curve). The value taken by the radiative flux is more than 20 times the expected value. This seems to be caused by a numerical instability. This can result in f >  1 during the iterations of our Newton–Raphson implicit scheme, which is not physically admissible. However, using the well-balanced modification of the source term proposed by Eq. (18) (blue curve), the spurious peak no longer appears and the constant steady state is reached.

Fig. 6. Simulation of a steady state with a jump of opacity. The opacity is piecewise constant, and a jump is located at x = 0.5 cm (gray line). This figure shows a snapshot of the radiative flux at time tf = 10−3 s.

Using the standard discretization of the source term in Eq. (16), one can show that the numerical scheme is unconditionally stable in that the error between the numerical solution and the exact solution goes to zero as Δx and Δt go to zero. The spurious peak seems to be due to a lack of precision in the integration of the source term. Using Eq. (18), the source term is defined at the interfaces of the cells and balances the divergence of radiative pressure, which is also defined at the interfaces.

5.3. Beam

We now perform the same two-dimensional test as in Richling et al. (2001), González et al. (2007). The domain [−1, 1] × [−1, 1] is discretized with 128 × 128 cells. The initial temperature is T₀ = T_r = 300 K, the initial radiative flux is F_r = 0, and the opacity is σ = 0. At time t = 0, a beam with T = T_r = 1000 K is introduced with an angle of 45°. The beam is located at x = −1 and y ∈ [−0.875, −0.75]. Because we are in the free-streaming regime, we cannot use large time-steps. For performance reasons, we use the semi-implicit scheme.

Because there is no opacity, the beam propagates in the vacuum, and we expect it to cross the box without dispersion. The direction of the beam is not along the mesh axis; we use this test to quantify the numerical diffusion. Figure 7 shows the radiative energy at steady state, with the eigenvalues fixed to ±c and with computed eigenvalues. Because there is no opacity, neither the asymptotic correction nor the well-balanced source term affect the result, and we recover the same result as in González et al. (2007). Figure 8 shows the horizontal cut at the middle height. The beam introduced at the boundary is sampled over 8 cells. Ideally, without any numerical diffusion, we would expect the width of the beam to remain exactly 8 cells. With the computed eigenvalues, we can keep the numerical diffusion under control. Using fixed eigenvalues, the width of the beam at middle height is approximately 30 cells, whereas it is only 24 cells with calculated eigenvalues.

Fig. 7. Beam simulation, showing radiative energy. The eigenvalues are fixed to ±c (upper panel) or calculated with Eq. (14) (lower panel).

Fig. 8. Beam simulation, showing a horizontal cut in Fig. 7 at the middle height.

The main difference in this test between ARK-RT and HERACLES (González et al. 2007) is the computation of the eigenvalues. We use the exact eigenvalues given by Eq. (14) from Berthon et al. (2007), whereas in González et al. (2007), to save computational time, the eigenvalues are computed once at the beginning of the simulation and then interpolated. However, this approximation does not impact the result.

5.4. Shadow

Let us now consider a two-dimensional test with source terms. Following Hayes & Norman (2003), González et al. (2007), we consider a shadow test. The computational domain is a cylinder of length L = 1 cm and radius R = 0.12 cm, and is discretized with 280 × 80 cells. A spheroid clump is located at the center of the box, on the symmetric axis: (z_c, r_c)=(0.5, 0). The extension of the clump is (z₀, r₀) = (0.1, 0.06). Initially, the medium is at equilibrium with the radiation, with T₀ = T_r = 290 K. We consider a homogeneous gas, with ρ₀ = 1 g cm⁻³, except for the clump with density ρ₁ = 100ρ₀. The boundary of the clump is smoothed: $\rho (r,z)={\rho }_0+\frac{{\rho }_1-{\rho }_0}{1+exp\mathrm{\Delta }}$ with $\mathrm{\Delta }=10({\left(\frac{z-z_c}{z_0}\right)}^2+{\left(\frac{r-r_c}{r_0}\right)}^2-1)$. The opacity in the medium is $\sigma ={\sigma }_0{\left(\frac{T}{T_0}\right)}^{-3.5}{\left(\frac{\rho }{{\rho }_0}\right)}^2$ with σ₀ = 0.1 cm⁻¹. At time t = 0, a source is lighted at the left boundary with T_r = 1740 K and the reduced flux is set to f = 1. Because f is close to 1 in the free-streaming regime, we encounter f >  1 in the simulation. To tackle this issue, we use the nonwell-balanced scheme: the radiative flux source term is discretized using Eq. (16). Because we are in the free-streaming regime, we cannot use large time-steps. For performance reasons, we use the semi-implicit scheme. To recover the same result as in González et al. (2007), we use ${\lambda }_{i+\frac{1}{2}}^{+}=max(0.1\times c,{\lambda }_{max})$ and ${\lambda }_{i+\frac{1}{2}}^{-}=min(-0.1\times c,{\lambda }_{min})$, where λ_max and λ_min are given by Eq. (14).

Figure 9 shows the radiative temperature at the final time t_f = 10⁻¹⁰ s with different closure relations: the P₁ model, the M₁ model with fixed eigenvalues, and the M₁ model with computed eigenvalues. Because of the high opacity in the clump, the light does not cross it and we expect the shadow behind it to remain stable.

Fig. 9. Shadow simulation, showing snapshots of the radiative temperature at time tf = 10−10 s with different closure relations: P1 model (upper panel), M1 model with fixed eigenvalues (middle panel), and M1 model with computed eigenvalues (lower panel).

As in Hayes & Norman (2003), González et al. (2007), we plot the radial profile of the radiative temperature at the right boundary (Fig. 10). Using the P₁ model, the radiative pressure is isotropic, and therefore the photons go around the obstacle immediately, heating the whole domain. Using the M₁ closure relation, the shadow is better preserved, and the temperature behind the obstacle remains at its initial value, 290 K. As the opacity remains rather low outside of the clump and the light has not crossed the obstacle, the asymptotic correction has no impact on the result. Because the boundary of the clump is smoothed, the transition between the optically thick and thin medium is less sharp than in Sect. 5.2 and the well-balanced modification of the source term is not necessary.

Fig. 10. Shadow simulation, showing the radial profiles of the radiative temperature at time tf = 10−10 s with different closure relations: P1 model, M1 model with fixed eigenvalues, and M1 model with computed eigenvalues.

5.5. Radiative shocks

We finally consider radiative shocks: the gas and the radiation exchange energy and momentum. Following Ensman (1994), Hayes & Norman (2003), González et al. (2007), we consider a one-dimensional homogeneous medium, with ρ = 7.78 × 10⁻¹⁰ g cm⁻³ and σ = 3.1 × 10⁻¹⁰ cm⁻¹. We consider a perfect gas with an adiabatic coefficient $\gamma =\frac{7}{5}$ and a mean molecular weight μ = 1. The length of the domain is 7 × 10¹⁰ cm. It is discretized with 400 cells. The initial temperature at the left boundary is set to 10 K and is increased by 0.25 K per cell. Initially, the radiation is at equilibrium with the gas. The left boundary condition is reflective, and the initial velocity of the fluid is set to u₀. According to the value of u₀, the shock will be subcritical or supercritical; see González et al. (2007) for more details. To compare our results with those of Ensman (1994), Hayes & Norman (2003), González et al. (2007), we plot the temperature as a function of x_i = x − u₀t.

5.5.1. Subcritical shock

We first consider a subcritical shock, the initial velocity is set to u₀ = −6 km s⁻¹. Figure 11 shows the gas temperature, the radiative temperature, and the reduced flux at three different times: 1.7 × 10⁴ s, 2.8 × 10⁴ s, and 3.8 × 10⁴ s. As expected, the gas and the radiation are not at equilibrium, before nor after the shock. The gas temperature reaches 1135 K, as in González et al. (2007), whereas it is only 850 K in Ensman (1994).

Fig. 11. Subcritical shock simulation, showing snapshots of gas temperature, radiative temperature, and reduced flux at different times: 1.7 × 104 s, 2.8 × 104 s, and 3.8 × 104 s.

5.5.2. Supercritical shock

We now consider a supercritical shock, where the initial velocity is set to u₀ = −20 km s⁻¹. Figure 12 shows the gas temperature, the radiative temperature, and the reduced flux at three different times: 4 × 10³ s, 7.5 × 10³ s, and 1.3 × 10⁴ s. As in González et al. (2007), the radiative temperature is the same as the matter temperature on both sides of the shock. The gas and the radiation are therefore at equilibrium. The radiative precursor is larger than that of the subcritical shock, as intended, and the temperature reaches 5000 K, as in Ensman (1994). We also recover the Zel’dovich spike.

Fig. 12. Supercritical shock simulation, showing snapshots of gas temperature, radiative temperature, and reduced flux at different times: 4.0 × 103 s, 7.5 × 103 s, and 1.3 × 104 s.

5.6. Expansion of H II region

Now that we have confirmed the good behavior of the numerical scheme with both the asymptotic preserving and the well-balanced properties, we can apply it to a physical situation: the propagation of the ionization front in a massive pre-stellar dense core.

5.6.1. Model

We consider the early stage of the development of an H II region in a massive pre-stellar dense core (Churchwell 2002). We focus on a region of the dense core at about 100 AU from the massive young stellar object (YSO). This region has been heated by the YSO during the pre-main sequence phase. The temperature reached at this location by infrared heating is of the order of 1000 K and the transport of energy in this region can be dominated by convection. We have inferred the convective state of this region by computing thermal and adiabatic gradients based on observations of Herpin et al. (2009) (Fig. 7). High-energy photons emitted by the YSO when entering the main sequence start to ionize the surrounding gas. This triggers the propagation of an ionization front in a convective medium, and we are interested in the stability of such a front perturbed by the pre-existing convective motions.

The interaction of the ionizing photons with the gas is described by Eq. (5). The only photons able to ionize the gas are emitted by the YSO, that is, there is no local source of ionizing photons. Following Tremblin (2012), we need to modify this model to take into account photochemistry and thermal balance. We define the fraction of ionization X = n_H⁺/n_H where n_H = n_H⁺ + n_H⁰, n_H⁺ is the number of ionized atoms and n_H⁰ is the number of cold atoms. The evolution of the number of ionized atoms is simply the number of incoming photons that interact with the gas minus the number of ionized atoms that recombine (on the spot approximation, see Lesaffre 2002). Therefore,

$$\begin{array}{c}\hfill {\partial }_t(\rho X)+\mathrm{\nabla }·(\rho X\mathit{U})={\sigma }_{\gamma }{F}_{\gamma }{n}_{\mathrm{H}}(1-X)-\beta X^2{n}_{\mathrm{H}}^2,\end{array}$$(29)

where F_γ is the number of incoming photons per unit of surface and time, σ_γ is the average cross-section at the temperature of the star, and β gives the recombination rate: β = 2 × 10⁻¹⁰T^−0.75 with T the temperature of thermodynamic equilibrium (Black 1981).

The thermal balance is the difference between the heating rate and the cooling rate. The extra energy of the absorbed photons is converted into kinetic energy of electrons; this is the only source of heating during the ionization, and therefore the heating rate is given by (1 − X)n_HF_γσ_γe_γ. In this simplified model, the equilibrium temperature is obtained from the balance between the heating from the ionization and the cooling from the recombination. We do not consider any other effects, such as metal cooling. Therefore, we take e_γ = 1 eV (Lesaffre 2002) to recover the observed temperature around 1000 K. From Tremblin (2012), the cooling rate is given by $\beta X^2{n}_{\mathrm{H}}^2{k}_{b}T/(\gamma -1)$. We also add a term of Newtonian forcing, ${\partial }_t{T}_{\mathrm{g}}=\frac{T_{\mathrm{g}}-T_{\mathrm{forcing}}}{{\tau }_{\mathrm{forcing}}}$, to trigger convection. T_forcing is the equilibrium temperature profile, depending on space, and τ_forcing is the relaxation timescale. The gas temperature will relax toward the equilibrium temperature profile T_forcing.

By writing cE_r = F_γe_γ, ρ = n_Hm_H and σ = σ_γn_H, we finally have to solve the following system:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\partial }_t\rho +\mathrm{\nabla }·(\rho \mathit{u})=0\hfill \\ {\partial }_t(\rho \mathit{u})+\mathrm{\nabla }·(\rho \mathit{u}\otimes \mathit{u}+p\mathtt{I})=\rho \mathit{g}+\frac{\sigma (1-X)}{c}{\mathit{F}}_r\hfill \\ {\partial }_t(\rho E)+\mathrm{\nabla }·((\rho E+p)\mathit{u})=\rho \mathit{g}·\mathit{u}+c\sigma (1-X)E_r\hfill \\ -\beta \frac{{\rho }^2{X}^2}{m_{\mathrm{H}}^2}\frac{k_b{T}_{\mathrm{g}}}{\gamma -1}-\frac{T_{\mathrm{g}}-T_{\mathrm{forcing}}}{{\tau }_{\mathrm{forcing}}}\hfill \\ {\partial }_t{E}_r+\mathrm{\nabla }·{\mathit{F}}_r=-c\sigma (1-X)E_r\hfill \\ {\partial }_t{\mathit{F}}_r+\mathrm{\nabla }·{\mathtt{P}}_r=-c\sigma (1-X){\mathit{F}}_r\hfill \\ {\partial }_t(\rho X)+\mathrm{\nabla }·(\rho X\mathit{u})=\frac{\sigma (1-X)c{E}_r{m}_{\mathrm{H}}}{e_{\gamma }}-\frac{\beta {\rho }^2{X}^2}{m_{\mathrm{H}}}.\hfill \end{array}\end{array}$$(30)

In this test, we use the M₁ solver with the asymptotic correction presented in Sect. 3.1.3, but we do not use the well-balanced discretization of the source term because of stability issues that we discuss in Sect. 6.1.

5.6.2. Setup

We consider a square domain of 5 AU on all sides. We use a setup similar to Padioleau et al. (2019) for compressible convection simulations. The temperature is set to 500 K and 1000 K at the top and bottom of the box, respectively. The initial temperature is a linear interpolation between the top and the bottom of the box. It is also the forcing temperature profile T_forcing. We take τ_forcing = 10⁷ s. These parameters are chosen to trigger the initial convective motions. We also set the pressure at the bottom: 10⁻³ dyne ⋅ cm⁻² (Herpin et al. 2009). The density and the pressure are linked by the ideal gas law: $p=\frac{\rho k_b{T}_{\mathrm{g}}}{m_{\mathrm{H}}\mu }$, where μ is the mean molecular weight. The nonionized medium is made of hydrogen, with μ₁ = 1. When the medium is fully ionized, it is made of atomic nucleus and electrons, and therefore has twice as many particles for the same mass. Because the distribution of nucleus and electrons is homogeneous, the mean molecular weight is μ₂ = 0.5. When the medium is partially ionized, we take μ as the mean of the previous values balanced by the fraction of ionization, that is, μ = (1 − X)μ₁ + Xμ₂. The density is initialized with the recursive formula $p_{i+1}-p_i=\frac{1}{2}({\rho }_i+{\rho }_{i+1})g\mathrm{\Delta }z$, which is the discrete version of the hydrostatic balance ∇p = −ρg.

We impose Neumann boundary conditions for the temperature. The pressure and density are imposed by an extrapolation of the hydrostatic balance. Because the hydrodynamics solver is well-balanced for the gravity, this configuration will remain static, even if the initial condition is unstable. The hydrostatic equilibrium is destabilized with a velocity mode perturbation of the form

$$\begin{array}{c}\hfill \begin{array}{c}u(x,y)=2·{10}^{-4}{c}_{s}sin\left(2\pi \frac{x-x_{\mathrm{mid}}}{L_x}\right)sin\left(\pi \frac{y-y_{\mathrm{mid}}}{L_y}\right)\hfill \\ v(x,y)=2·{10}^{-4}{c}_{s}cos\left(2\pi \frac{x-x_{\mathrm{mid}}}{L_x}\right)cos\left(\pi \frac{y-y_{\mathrm{mid}}}{L_y}\right),\hfill \end{array}\end{array}$$(31)

with c_s the speed of sound, x_mid = y_mid = 2.5 AU and L_x = L_y = 5 AU. Without any interaction with the ionizing photons, the convective motions are stationary.

The opacity is set to $\sigma =\frac{{\sigma }_{\gamma }\rho }{m_{\mathrm{H}}}$ with σ_γ = 6 × 10⁻¹⁸ cm² (Lesaffre 2002). The radiative energy and flux are set to zero and the medium is not ionized (X = 0). We initialize the hydrodynamics variables with the steady state described above. At time t = 0, the bottom boundary of the region is ionized: the reduced flux is set to a value of 1 and the radiative energy is set to $\frac{F_{\ast }{e}_{\gamma }}{c}$ with F_* = 3 × 10¹⁷ cm⁻² s⁻¹ in the boundary. The boundary conditions for the hydrodynamics variables remain unchanged.

5.6.3. Results

As the initial condition is such that E_r is close to zero, this can easily create some spurious values such that f >  1. This is clearly a numerical artifact induced by the very low value of the radiative energy in regions where no ionizing photons are present. Even with a centered discretization of the radiative flux source term and without the asymptotic correction, we still encounter f >  1 during the simulation. Because of this problem, we impose f = 1 in the computation of the Eddington tensor in the cells where f >  1.

Figure 13 shows the evolution of the position of the ionization front as a function of time. With and without the initial convective rolls, the position of the ionization front oscillates around an equilibrium position, between 0.3 AU and 0.4 AU. The oscillations around the equilibrium are expected and have been observed with simpler models (Tremblin et al. 2012).

Fig. 13. Evolution of the position of the ionization front as a function of time, with and without the initial velocity perturbation.

Figures 14 and 15 show the ionization front at time t = 6 × 10⁸ s and at the final time t_f = 10¹⁰ s. With and without the initial convective rolls, a numerical noise appears as a consequence of the long timescales. Because of the numerical noise, some lack of symmetry can appear, such as is seen in the left panel of Fig. 14. The fraction of ionization, which is always between 0 and 1, reaches values between 10⁻¹² and 10⁻⁶. The effect of the preconditioner and the MPI domain decomposition is discussed in Appendix D. However, the numerical noise does not affect the position of the ionization front.

Fig. 14. Snapshots of the fraction of ionization and the velocity field at time t = 6 × 108 s without the initial velocity perturbation (left panel) and with it (right panel).

Fig. 15. Snapshots of the fraction of ionization and the velocity field at the final time tf = 1010 s without the initial velocity perturbation (left panel) and with it (right panel).

The stability of the ionization front is an issue that has been discussed for a long time in the literature (Mizuta et al. 2008). For example, 3D simulations of the expansion of a spherical ionization front in 3D Cartesian grids have shown instabilities either on the axis of the grid or in the diagonal depending on the numerical scheme (see Fig. A3 Bisbas et al. 2015). The dependence of the instability on the grid casts doubt on the possibility of a physical regime. Our test case shows that even with convective motions of large amplitude, the ionization front remains stable.

6. Discussion and conclusion

We present a new radiation hydrodynamics code. The radiative transfer is described with a two-moment model and uses the M₁ closure relation. We first discuss some limitations of our work before presenting our conclusions.

6.1. Well-balanced discretization of the source term

In Sect. 3.1.2 we propose a well-balanced discretization of the source term on the radiative flux equation. This discretization allows us to properly capture the steady state with constant flux and with a discontinuity of opacity. However, this discretization can lead to spurious oscillations in the radiative flux, a problem that we have encountered in the test case for the expansion of H II regions. Although we have changed the discretization of the source term to achieve a well-balanced property, our integration of the hyperbolic part on the source term is still split into two steps. Such a splitting strategy might be unstable if the source term is not taken into account in the hyperbolic part. A possible solution to this problem would be to incorporate the source term in a Lagrange-projection-like scheme such as Buet & Despres (2008). Such a strategy might be necessary to treat the radiation hydrodynamics problem of cloud interfaces.

6.2. Asymptotic limit for radiation hydrodynamics

In Sect. 3.1.3, we present an asymptotic correction which allows us to capture the asymptotic behavior, whereas the solution given by a standard scheme is dominated by numerical diffusion. The asymptotic correction uses the numerical diffusion to recover the physical one in a static fluid. Nevertheless, this scheme does not capture the asymptotic regime in a moving fluid, as presented in Appendix A. Most of the schemes proposed in the literature do not preserve this asymptotic regime (e.g., González et al. 2007; Berthon & Turpault 2011). The diffusive regime depends on the material velocity; our scheme cannot reach it. A possible solution would be to limit the numerical diffusion with a correction similar to a low Mach correction, as in Chalons et al. (2016), in conjunction with a cell-centered discretization of the source term, as proposed by the Lagrange-projection scheme of Buet & Despres (2008).

6.3. Conclusion

The model for radiation hydrodynamics proposed in this paper has the correct behavior in both free-streaming and diffusive limits. It is discretized with an asymptotic preserving scheme. This asymptotic correction is important to capture the correct behavior in the diffusive limit and to preserve the free-streaming regime in a static fluid.

We take advantage of the libraries Kokkos and Trilinos to reach high-performance computing and to solve linear systems. This approach allows us to take advantage of different architectures and to use large time-steps for radiative transfer. Using the implicit solver is profitable as soon as the time-step given by the hydrodynamics is 50−100 times larger than the explicit time-step for radiative transfer, depending on the preconditioner and the architecture.

The solver is then coupled with the hydrodynamics code that implements an all-regime and well-balanced solver for hydrodynamics with gravity. The tests performed show that ARK-RT is well suited to studying many astrophysical problems. We used this code to study the development of H II regions in massive pre-stellar dense cores, especially the propagation of the ionization front in the presence of convection. We show that even with the destabilizing effect of convection, the ionization front is strongly stable against perturbations. A linear stability analysis similar to the Rayleigh-Taylor instability but including a source term could provide more insight into this behavior.

Further work will consist of the development of a numerical scheme that preserves the admissible states while keeping the asymptotic preserving and well-balanced properties. We will then be able to take full advantage of the largest next-generation, massively parallel architectures to study atmospheric physics with opacity interfaces, among other things.

Acknowledgments

P. T. acknowledges supports by the European Research Council under Grant Agreement ATMO 757858. The authors thank the referee for his useful comments, especially about the test on the expansion of H II regions. H. B. thanks also Rodolphe Turpault for his welcome in Bordeaux and the helpful discussion.

References

Appendix A: Diffusive limit for radiation hydrodynamics

As in Sect. 2.4, we consider the diffusive limit with the P₁ closure relation. We introduce a rescaling parameter ε to write the time and opacity as t̃ = εt and $\stackrel{\sim }{\sigma }=\epsilon \sigma$, respectively. Because the velocity of the fluid is smaller than the speed of light ($\frac{u}{c}\ll 1$), we also rescale it as $\stackrel{\mathbf{\sim }}{\mathit{u}}=\frac{\mathit{u}}{\epsilon }$. Let us focus on the equations describing the evolution of the radiative variables:

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\epsilon }^2{\partial }_{\stackrel{\sim }{t}}{E}_r+\epsilon \mathrm{\nabla }·{\overrightarrow{F}}_r\hfill & =c\stackrel{\sim }{\sigma }(a_r{T}_g^4-E_r)\hfill \\ \hfill & +\epsilon \frac{\stackrel{\sim }{\sigma }}{c}\overrightarrow{\stackrel{\sim }{u}}·{\overrightarrow{F}}_r+{\epsilon }^2\frac{4}{3}\frac{\stackrel{\sim }{\sigma }}{c}{E}_r\overrightarrow{\stackrel{\sim }{u}}·\overrightarrow{\stackrel{\sim }{u}}(\mathrm{A}.1\mathrm{a})\hfill \\ {\epsilon }^2{\partial }_{\stackrel{\sim }{t}}{\overrightarrow{F}}_r+\epsilon \frac{c^2}{3}\mathrm{\nabla }{E}_r\hfill & =-\stackrel{\sim }{\sigma }c{\overrightarrow{F}}_r+\epsilon \frac{c}{3}\stackrel{\sim }{\sigma }\overrightarrow{\stackrel{\sim }{u}}{E}_r\hfill \\ \hfill & +\epsilon c\stackrel{\sim }{\sigma }\overrightarrow{\stackrel{\sim }{u}}{a}_r{T}_g^4+{\epsilon }^2\frac{\stackrel{\sim }{\sigma }}{c}\overrightarrow{\stackrel{\sim }{u}}·(\overrightarrow{\stackrel{\sim }{u}}·{\overrightarrow{F}}_r).(\mathrm{A}.1\mathrm{b})\hfill \end{array}\end{array}$$

By expanding Eqs. (A.1a) and (A.1b) at zero order, we have

$$\begin{array}{c}\hfill \{\begin{array}{c}{E}_{r,0}=a_r{T}_{\mathrm{g},0}^4\hfill \\ {\mathit{F}}_{r,0}=\mathbf{0}.\hfill \end{array}\end{array}$$(A.2)

Expanding Eq. (A.1b) at first order leads to

$$\begin{array}{c}\hfill {\mathit{F}}_{r,1}=-\frac{c}{3\stackrel{\sim }{\sigma }}\mathrm{\nabla }{E}_{r,0}+\frac{4}{3}{E}_{r,0}{\stackrel{\mathbf{\sim }}{\mathit{u}}}_0.\end{array}$$(A.3)

Finally, looking at the radiative energy and the gas internal energy at second order, source terms cancel each other, and only the divergence of the radiative flux at first order remains, and we have

$$\begin{array}{c}\hfill {\partial }_{\stackrel{\sim }{t}}(\rho c_v{T}_{\mathrm{g},0}+E_{r,0})-\mathrm{\nabla }·\left(\frac{c}{3\stackrel{\sim }{\sigma }}\mathrm{\nabla }{E}_{r,0}\right)=-\frac{4}{3}\mathrm{\nabla }\left(E_{r,0}{\stackrel{\mathbf{\sim }}{\mathit{u}}}_0\right).\end{array}$$(A.4)

One can also look at Eq. (A.1a) at second order. Expanding Eq. (A.1a) at second order gives

$$\begin{array}{cc}& {\partial }_{\stackrel{\sim }{t}}{E}_{r,0}-\mathrm{\nabla }\left(\frac{c}{3\stackrel{\sim }{\sigma }}\mathrm{\nabla }{E}_{r,0}\right)=c\stackrel{\sim }{\sigma }(6a_r{T}_{\mathrm{g},0}^2{T}_{\mathrm{g},1}^2+4a_r{T}_{\mathrm{g},0}^3{T}_{\mathrm{g},2}-E_{r,2})\hfill \\ \hfill & -\frac{4}{3}\mathrm{\nabla }·\left(E_{r,0}{\mathit{u}}_0\right)\hfill \\ \hfill & -\frac{1}{3}{\mathit{u}}_0·\mathrm{\nabla }{E}_{r,0}+2\frac{4}{3}\frac{\sigma }{c}{E}_{r,0}{\mathit{u}}_0^2.\hfill \end{array}$$(A.5)

We recover Eq. (43) of Krumholz et al. (2007). The term $c\stackrel{\sim }{\sigma }(6a_r{T}_{\mathrm{g},0}^2{T}_{\mathrm{g},1}^2+4a_r{T}_{\mathrm{g},0}^3{T}_{\mathrm{g},2}-E_{r,2})$ is the development at second order of the term κ₀(4πB − cE). Because we do not neglect any terms $O\left(\frac{\stackrel{\mathbf{\sim }}{\mathit{u}}}{c}\right)$, some coefficients are slightly different. See the discussion in Krumholz et al. (2007) for the importance of the term $\frac{4}{3}\frac{\stackrel{\sim }{\sigma }}{c}{E}_{r,0}{\stackrel{\mathbf{\sim }}{\mathit{u}}}_0^2$.

Appendix B: Von Neumann stability analysis for the well-balanced modification of the source term

For simplicity, we split Eq. 3b into a pure hyperbolic problem ∂_tF_r + c²∇ · P_r = 0 and a source problem ∂_tF_r = −cσF_r. We focus on the one-dimensional source problem, with periodic boundary conditions on the domain [0, T] × [0, 1] with T the final time. The following can easily be extended to an arbitrary space interval. Because we use periodic boundary conditions, we can apply the von Neumann stability analysis (see e.g., Anderson 1995), based on the decomposition of the numerical solution into Fourier series. Let us recall that, using the well-balanced modification of the source term, the source problem is discretized as

$$\begin{array}{c}\hfill F_j^{n+1}+r_{j-\frac{1}{2}}{F}_{j-1}^{n+1}+(r_{j-\frac{1}{2}}+r_{j+\frac{1}{2}})F_j^{n+1}+r_{j+\frac{1}{2}}{F}_{j+1}^{n+1}=F_j^n,\end{array}$$(B.1)

with $r_{j+\frac{1}{2}}=\frac{c{\sigma }_{j+\frac{1}{2}}\mathrm{\Delta }t}{4}$. We define the function Fⁿ, piecewise constant, such that

$$\begin{array}{c}\hfill F^n(x)=\{\begin{array}{cc}{F}_j^n\hfill & \mathrm{if}{x}_{j-\frac{1}{2}}<x<x_{j+\frac{1}{2}}\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\end{array}$$(B.2)

This function is then extended to ℝ by periodicity. Fⁿ can now be expanded in a Fourier series:

$$\begin{array}{c}\hfill F^n(x)={\displaystyle \sum _{k\in \mathbb{Z}}}{\widehat{F}}^n(k)e^{2ik\pi x},\end{array}$$(B.3)

with

$$\begin{array}{c}\hfill {\widehat{F}}^n(k)={\displaystyle {\int }_0^1}{F}^n(x)e^{-2ik\pi x}\mathrm{d}x.\end{array}$$(B.4)

We can define the 2-norm of the function Fⁿ:

$$\begin{array}{cc}& ||F^n{||}_2={\left({\displaystyle {\int }_0^1}{\left(F^n(x)\right)}^2\mathrm{d}x\right)}^{\frac{1}{2}}\hfill \\ \hfill & ={\left({\displaystyle \sum _{k\in \mathbb{Z}}}{|{\widehat{F}}^n(k)|}^2\right)}^{\frac{1}{2}}\hfill \\ \hfill & ={\left({\displaystyle \sum _{j=1}^J}\mathrm{\Delta }x{\left(F_j^n\right)}^2\right)}^{\frac{1}{2}}.\hfill \end{array}$$(B.5)

We apply the Fourier transform to Eq. (B.1):

$$\begin{array}{c}\hfill {\widehat{F}}^{n+1}(k)(1+r_{j-\frac{1}{2}}{e}^{-2ik\pi \mathrm{\Delta }x}+r_{j-\frac{1}{2}}+r_{j+\frac{1}{2}}+r_{j+\frac{1}{2}}{e}^{2ik\pi \mathrm{\Delta }x})={\widehat{F}}^n(k).\end{array}$$(B.6)

We define the amplification factor A(k) as

$$\begin{array}{c}\hfill A(k)=\frac{1}{1+r_{j-\frac{1}{2}}{e}^{-2ik\pi \mathrm{\Delta }x}+r_{j-\frac{1}{2}}+r_{j+\frac{1}{2}}+r_{j+\frac{1}{2}}{e}^{2ik\pi \mathrm{\Delta }x}},\end{array}$$(B.7)

and we then have ${\widehat{F}}^{n+1}(k)=A(k){\widehat{F}}^n(k)$. By induction, we have ${\widehat{F}}^n(k)={(A(k))}^n{\widehat{F}}^0(k)$. The coefficient ${\widehat{F}}^n(k)$ remains bounded if and only if |A(k)| ≤ 1. In this case, for all k ∈ ℤ, $|{\widehat{F}}^{n+1}(k)|\le |{\widehat{F}}^n(k)|$. Therefore, ||F^n + 1||₂ ≤ ||Fⁿ||₂ ≤ ||F⁰||₂ and the scheme is unconditionally stable. We now have to prove that |A(k)| ≤ 1: