Edouard Audit

Pascal Tremblin

Matthias González

H Bloch

S Fromang

T Padioleau

P Kestener

S Kokh

Radiative Rayleigh

Michael A Heroux

Roscoe A Bartlett

Vicki E Howle

Robert J Hoekstra

Jona- Than J Hu

Tamara G Kolda

Richard B Lehoucq

Kevin R Long

Roger P Paw- Lowski

Eric T Phipps

Andrew G Salinger

Heidi K Thornquist

Je veux remercier toutes les personnes grâce à qui ces trois années se sont aussi

Le transfert radiatif décrit les interactions entre le rayonnement et le milieu dans lequel il se propage. Il apparaît dans une grande gamme d'applications physiques, de l'astrophysique à la médecine. Dans de nombreuses situations astrophysiques, le rayonnement est un processus important qui interagit avec le gaz environnant, par exemple dans les atmosphères d'(exo)planètes [Thomas and [START_REF] Gary | Radiative Transfer in the Atmosphere and Ocean[END_REF], dans les étoiles massives [START_REF] Kuiper | CIRCUMVEN-TING THE RADIATION PRESSURE BARRIER IN THE FORMATION OF MASSIVE STARS VIA DISK ACCRETION[END_REF]Mignon-Risse et al., 2020], les régions H ii (Spitzer 1978 et fig. 1.1), jusqu'à la réionisation cosmique [Stiavelli, 2009]. La figure 1.1 montre la nébuleuse d'Orion en couleurs visibles. Au centre, une étoile massive émet de la lumière. Les photons les plus énergisants ionisent la région environnante, la région H ii. Les photons moins énergisants traversent cette zone sans interaction et dissocient les molécules de dihydrogène dans la région de photodissociation (Photodissociation Region, PDR) extérieure. L'expansion du gaz chaud dans la région H ii comprime la PDR ; les échanges d'énergie entre le rayonnement et la matière affectent l'hydrodynamique du système.

L'hydrodynamique radiative est le cas particulier où le rayonnement interagit avec un fluide. La lumière peut être absorbée, ce qui va chauffer le fluide. Du rayonnement peut également être émis par le fluide ou par une autre source. Le rayonnement peut être diffusé, ce qui va changer sa direction de propagation et éventuellement sa fréquence [Mihalas and Mihalas, 1984].

La distance parcourue par les photons avant d'être absorbés ou diffusés est modulée par l'opacité du milieu environnant. D'une part, le milieu peut être optiquement épais, le rayonnement et l'hydrodynamique sont fortement couplés. D'autre part, le milieu peut être optiquement mince, la lumière et le fluide n'interagissent pas.

La théorie décrivant l'interaction entre le rayonnement et la matière a été formalisée par Chandrasekhar [1960]. L'intensité spécifique obéit à une équation intégrodifférentielle, qui ne peut pas être résolue analytiquement dans la plupart des cas.

Comme les expériences et les observations, les simulations numériques peuvent améliorer la compréhension du système physique simulé. Parfois, les expériences ne peuvent pas être réalisées, car elles seraient trop dangereuses (accident), trop longues (climatologie), inatteignables (astrophysique) ou interdites (essais nucléaires). Les simulations numériques peuvent aussi être plus simples : par exemple, il est plus facile de réaliser des simulations en mécanique des fluides plutôt que des expériences en soufflerie.

Une étape supplémentaire est nécessaire, pour passer du modèle continu à une solution approchée. C'est la discrétisation de l'équation, pour obtenir un schéma numérique. Nous présentons brièvement quelques méthodes parmi les nombreuses techniques qui existent pour résoudre l'équation du transfert radiatif, ainsi que leurs limitations.

Le tracé de rayons résout le problème en suivant la propagation de rayons à travers le fluide [Wise and Abel, 2011]. Bien que cette méthode soit très précise, elle est très coûteuse quand elle est couplée à un code d'hydrodynamique. Le nombre de degrés de liberté est proportionnel au nombre de cellules spatiales multiplié par le nombre de sources de rayonnement.

En utilisant la méthode de Monte-Carlo, nous suivons la propagation de paquets de photons et leurs interactions avec le fluide [Roth and Kasen, 2014]. Elle est basée sur un processus stochastique, rendant la méthode précise, mais difficile à coupler avec un code d'hydrodynamique basé sur une grille. Quand le nombre de paquets de photons devient trop petit dans une région, un bruit numérique peut apparaître et polluer les résultats de la simulation, rendant la méthode inutilisable pour étudier les problèmes de stabilité d'interface.

Une autre approximation largement utilisée est la méthode des moments de l'intensité spécifique (en moyennant sur les directions de propagation des photons, Levermore 1984). Dans cette thèse, nous considérons les trois premiers moments de l'intensité spécifique : l'énergie radiative E r , le flux radiatif F r et la pression radiative P r . Une relation de fermeture est nécessaire pour exprimer la pression radiative en fonction de l'énergie radiative et du flux radiatif. Celle utilisée est le modèle M 1 [Levermore, 1984 ;Dubroca and Feugeas, 1999] pour ses bonnes propriétés dans les milieux optiquement épais et optiquement minces, ainsi que le compromis entre temps de calcul et précision. Remarquons que l'énergie radiative est positive E r ą 0 et que le flux réduit vérifie ||Fr|| cEr ď 1, où c est la vitesse de la lumière. L'énergie radiative ne peut pas être transportée plus vite que la vitesse de la lumière.

Le modèle M 1 est un système hyperbolique, il peut être discrétisé avec la méthode des volumes finis. Tous les schémas numériques présentés dans cette thèse sont inspirés des schémas numériques développés pour les équations d'Euler compressibles. Cependant, les spécificités du transfert radiatif et du modèle M 1 complexifient les schémas.

On peut montrer que, dans la limite des longues échelles de temps et des hautes opacités (milieu optiquement épais), le modèle M 1 dégénère vers une équation de diffusion. Un schéma qui préserve le régime asymptotique capture ce comportement.

Un solveur explicite pour le transfert radiatif serait restreint par une condition de Courant-Friedrichs-Lewy (CFL), limité par la vitesse de la lumière. On obtiendrait un pas de temps très petit devant celui de l'hydrodynamique, qui est limité par la vitesse du son du fluide. Plusieurs méthodes ont été développées pour contourner ce problème.

Nous avons choisi d'utiliser un solveur implicite [González et al., 2007]. La temporalité du transfert radiatif est préservée, ce qui n'est pas le cas avec l'approximation de la vitesse de la lumière réduite (Reduced Speed of Light Approximation, RSLA, Gnedin and Abel 2001).

Notre but est de développer un schéma numérique avec ces trois propriétés :

-le schéma doit préserver le régime asymptotique pour capturer le comportement en temps long et à haute opacité ; -la discrétisation en temps doit être implicite, pour éviter que le pas de temps ne soit limité par la vitesse de la lumière ; -la solution doit être admissible, c'est-à-dire les états admissibles E r ą 0 et ||Fr|| cEr ď 1 doivent être préservés. Ceci est nécessaire pour des raisons physiques et numériques. De nombreuses heures de calcul peuvent être nécessaires pour étudier les problèmes astrophysiques. L'usage d'outils pour le calcul haute performance (High Performance Computing, HPC) est nécessaire. Ils permettent de réduire le coût de calcul et d'augmenter la taille du problème à résoudre. Plusieurs processeurs traitent différentes parties d'un programme ou différentes données : des problèmes plus grands peuvent être résolus plus vite.

Au cours des dernières décennies, différentes technologies ont émergé, telles que les processeurs multi-coeurs et les GPUs. Ceci soulève le problème de portabilité de performance. Un code peut être optimisé pour une architecture spécifique et l'exécuter sur une autre architecture donnera lieu à de mauvaises performances. Certaines bibliothèques, comme Kokkos, proposent des solutions.

Les solveurs implicites peuvent nécessiter la résolution de grands systèmes linéaires. Cela peut représenter la majorité du temps de calcul. Heureusement, l'algèbre linéaire pour le calcul haute performance a été étudiée au fil des années. Comme les systèmes linéaires que nous devons résoudre sont grands et creux, les méthodes directes sont inutilisables. Des solveurs itératifs préconditionnés ont été développés pour contourner ce problème [Saad, 2003]. La plupart de ces méthodes nécessitent de stocker la matrice. Dans cette thèse, nous utilisons la bibliothèque Trilinos [START_REF] Michael | An overview of the trilinos project[END_REF] car elle nous permet d'utiliser différentes architectures telles que les processeurs multi-coeurs et les GPUs. Elle fournit, entre autres, des préconditionneurs multigrille algébrique (Algebraic Multigrid, AMG). D'autres méthodes, qui ne nécessitent pas le stockage de la matrice, existent également. C'est notamment le cas des méthodes multigrille géométrique (Geometric Multigrid, GMG) qui ont été développées pour la résolution de systèmes, linéaires ou non linéaires, venant de la discrétisation d'équations différentielles [START_REF] Briggs | A Multigrid Tutorial[END_REF].

Cette thèse se concentre sur le développement et l'implémentation de schémas numériques pour l'hydrodynamique radiative, et non sur leurs applications à l'astrophysique.

Description des travaux 1.Transfert radiatif

Le premier chapitre est dédié à la présentation du modèle utilisé pour décrire le transfert radiatif. Nous établissons d'abord l'équation du transfert radiatif en effectuant un bilan d'énergie. Elle décrit l'évolution de l'intensité spécifique.

Ces deux schémas (deux étapes et une étape) préservent le régime asymptotique si nous utilisons les flux donnés par Buet and Despres [2008]. Ces flux prennent en compte les termes sources. Nous montrons également que ces schémas sont entropiques et préservent les états admissibles, E r ą 0 et ||Fr|| cEr ď 1. Les tests numériques standards montrent le bon comportement des schémas dans les limites de transport et de diffusion.

Méthode de Jacobi non linéaire et multigrille géométrique

Dans ce chapitre, nous présentons une étape supplémentaire vers le développement d'un schéma implicite qui préserve le régime asymptotique et les états admissibles. Nous discrétisons le modèle M 1 avec un solveur HLL implicite sans terme source. Au lieu de résoudre le système non linéaire avec la méthode de Newton-Raphson, nous utilisons la méthode de Jacobi présentée par [START_REF] Pichard | Martin Mathematiques appliquees et calcul scientifique Bordeaux[END_REF]. Nous montrons que la solution construite avec ce processus itératif préserve les états admissibles.

Cependant, cette méthode peut converger lentement quand la résolution augmente. Pour s'attaquer à ce problème, nous utilisons un algorithme de multigrille géométrique [START_REF] Briggs | A Multigrid Tutorial[END_REF]. Le problème est résolu sur une grille plus grossière, où la méthode de Jacobi converge plus vite. Nous introduisons un pseudo-temps [START_REF] Kifonidis | On multigrid solution of the implicit equations of hydrodynamics -experiments for the compressible euler equations in general coordinates[END_REF] pour préserver les états admissibles.

Les premiers tests numériques dans le régime de transport montrent l'efficacité de la méthode.

Bruno Dubroca and Jean-Luc Feugeas. Etude théorique et numérique d'une hiérarchie de modèles aux moments pour le transfert radiatif. Comptes Rendus de l'Académie des Sciences -Series I -Mathematics, 329 (10) :915 -920, 1999. ISSN 0764-4442. doi : https://doi.org/10.1016/S0764-4442(00)87499-6. URL http: //www.sciencedirect.com/science/article/pii/S0764444200874996. 3 Radiative transfer describes the interactions between the radiation and the surrounding medium. It appears in a wild range of physical situations, from astrophysics, up to medicine. In many astrophysical situations, radiation is an important process that interacts with the surrounding gas, e.g., in (exo) planet's atmospheres (e.g., Thomas and Stamnes 2002), massive stars (e.g., [START_REF] Kuiper | CIRCUMVEN-TING THE RADIATION PRESSURE BARRIER IN THE FORMATION OF MASSIVE STARS VIA DISK ACCRETION[END_REF]; Mignon-Risse et al. 2020), H ii regions (e.g., Spitzer 1978 and fig. 1.2), up to the cosmic reionization (e.g., Stiavelli 2009). Figure 1.2 shows the Orion nebula in visible colors. In the middle, a massive star emits light. Photons with the highest energy ionize the surrounding zone, the H ii region. Less energetic photons travel across this area without interaction and dissociate dihydrogen molecules in the outer photodissociation region (PDR). The expanding hot gas in the H ii region compresses the PDR; energy exchanges between light and matter affect the hydrodynamics of the system.

Introduction (english version)

Radiation hydrodynamics is the special case where the radiation interacts with a fluid. Light can be absorbed, thus heating the fluid. The radiation could also be emitted by the fluid or by another source. Radiation could be scattered, this will change its direction and perhaps its frequency [Mihalas and Mihalas, 1984].

The distance traveled by the photons before being absorbed or scattered is shaped by the opacity of the surrounding medium. On one hand, the medium can be optically thick, the radiation and the hydrodynamics are strongly coupled. On the other hand, the medium can be optically thin, the light and the fluid do not interact with each other.

The theory describing the interaction between radiation and matter was formalized by Chandrasekhar [1960]. The specific intensity obeys an integro-differential equation, which cannot be solved analytically in most cases.

As experiments and observations, numerical simulations can improve the understanding of the simulated physical system. Sometimes, experiments cannot be performed, should they be too dangerous (crash), too long (climatology), unreachable (astrophysics), or forbidden (nuclear test). Numerical simulations can also be simpler: for example, it is easier to perform fluid simulations instead of experiments in a wind tunnel.

Another step is required, to go from the continuous model to an approximate solution. It is the discretization of the equation, to obtain a numerical scheme. Let us now briefly present some methods and their limitations, among the numerous techniques that exist to solve the radiative transfer equation.

Ray tracing solves the problem by following the propagation of beams through the fluid (e.g., Wise and Abel 2011). Although this method is very precise, it is highly costly when coupled to hydrodynamics codes. The number of degrees of freedom scales with the number of spatial cells multiplied by the number of radiation sources.

Using a Monte-Carlo method, we follow the propagation of "photon energy packets" and their interactions with the fluid (e.g., Roth and Kasen 2014). It is based on a stochastic process, making this method accurate, but difficult to couple with a grid based hydrodynamics code. When the number of photon packets becomes too small in a region, a numerical noise can arise and pollute the simulation results, making the method unusable to study interface stability problems.

Another widely used approximation is moment models [Levermore, 1984] of the specific intensity (by averaging over the direction of propagation of photons). In this thesis, we 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 . A closure relation is needed to express the radiative pressure as a function of the radiative energy and the radiative flux. The one we use is the M 1 model [Levermore, 1984;Dubroca and Feugeas, 1999] for its good properties in both optically thick and thin media and the compromise between computational cost and precision. Let us notice that the radiative energy is positive E r ą 0 and the reduced flux verifies ||Fr|| cEr ď 1, where c is the speed of light. It ensures that the radiative energy cannot be transported faster than the speed of light.

The M 1 model is a hyperbolic system, it can be discretized with a finite volume method. All numerical schemes presented in this thesis are inspired by numerical schemes developed for the compressible Euler equations. However, the specific features of the radiative transfer and the M 1 model increase the complexity of the schemes.

One can show that, in the limit of long timescale and high opacity (optically thick medium), the M 1 model degenerates towards a diffusion equation. An asymptotic preserving scheme captures this behavior.

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 one, which is limited by the speed of sound of the fluid. Several methods have been developed to get around this problem. The one we have chosen is a time-implicit solver (e.g., González et al. 2007). The temporality of the radiative transfer will be preserved, which is not the case with the reduced speed of light approximation (RSLA, e.g, Gnedin and Abel 2001).

We aim at building a numerical scheme with these three properties:

-the scheme should be asymptotic preserving to capture the behavior in long timescale and high opacity; -the discretization in time should be implicit, to avoid a time step restricted by the speed of light; -the solution should be admissible, i.e., the admissible states E r ą 0 and ||Fr|| cEr ď 1 should be preserved. This is needed for both physical and numerical reasons.

To study astrophysical problems, many computational hours can be needed. This requires the use of High Performance Computing (HPC) tools. It allows reducing the computational cost and increasing the size of the problems to be solved. Multiple processing units handle different parts of the same program or data simultaneously: bigger problems can be solved faster.

In the last decades, different architectures have emerged, such as multi-cores, manycores, and GP-GPUs. This raises the problem of performance portability. A code can be optimized for a specific architecture and running it on a different architecture will result in bad performance. Some libraries, such as Kokkos, address this issue.

Implicit solvers can require solving large linear systems. It can represent most of the computational time. Fortunately, linear algebra for HPC has been investigated over the years. Because the linear system we will 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). Most of these methods require storing the matrix. In this work, we use the library Trilinos [START_REF] Michael | An overview of the trilinos project[END_REF] because it allows us to target different architectures, such as multi-core, many-core, and GP-GPUs. It provides, among others, Algebraic Multigrid (AMG) preconditioners. Other methods that do not require to store the matrix also exist. For example, Geometric Multigrid (GMG) methods have been developed to solve linear and nonlinear systems arising from the discretization of differential equations (e.g., [START_REF] Briggs | A Multigrid Tutorial[END_REF].

This thesis focuses on the development and the implementation of numerical schemes for radiation hydrodynamics, not their application to astrophysics.

Description of the work 1.4.1 Radiative transfer

The first chapter is dedicated to the presentation of the model used to describe the radiative transfer. We first derive the radiative transfer equation by computing the energy balance. It describes the evolution of the specific intensity.

This equation can be costly to solve numerically when the radiative transfer is coupled with the hydrodynamics. To reduce the computational cost, we use a moment method, in which the specific intensity is averaged over the direction of propagation and the frequency of the photons. This leads to solving a hyperbolic system involving radiative energy, radiative flux, and radiative pressure. This model requires specifying a closure relation. The one we use is the M 1 model. We express the radiative pressure as a function of the radiative energy and radiative flux by maximizing the radiative entropy.

This model has the correct behavior in both free-streaming and diffusive limits. We look at this system in the diffusive limit, i.e., large opacity and long timescale.

For physical applications, the model has to take into account the exchange of energy between the radiation and the hydrodynamics. We model this phenomenon by adding nonconservative terms. These terms depend on the opacity, we investigate again the behavior of the model in the diffusive limit, when the coupling to the hydrodynamics is taken into account.

A first asymptotic preserving solver

In this chapter, we present a first numerical scheme, time-implicit, and asymptotic preserving. It is based on an HLL solver, numerical fluxes are corrected to capture the asymptotic behavior. We also suggest a well-balanced modification of the source terms to reach a constant steady state with a jump of opacity.

Using a time-implicit solver, a nonlinear system has to be solved. It is done with the Newton-Raphson method and linear systems are solved using the library Trilinos.

The resolution of the radiation hydrodynamics problem is split into three steps. The first one is the update of the hydrodynamics quantities. The second one is the update of the radiative quantities. The third one is the addition of coupling source terms.

We perform standard tests in both free-streaming and diffusive limits, including purely radiative tests and radiation hydrodynamics ones. We also study the propagation of an ionization front in a massive pre-stellar dense core, and we show that the ionization front is strongly stable against perturbations even with destabilizing convective motions.

An all-regime-like asymptotic preserving scheme

Inspired by numerical schemes for hydrodynamics [Chalons et al., 2016], we propose in this chapter a second numerical scheme based on an all-regime solver. The M 1 model is first written in a form similar to the Euler equations. Diffusion and transport phenomena are then split into two subsystems. We use numerical fluxes given by Buet and Despres [2008] to discretize the diffusion system. We use an upwind scheme to discretize the transport step while ensuring a conservative update.

However, this leads to a scheme of stencil 2, with two steps. We also propose an all-regime scheme with only one step.

Both schemes (two steps and one step) are asymptotic preserving if we use numerical fluxes given by Buet and Despres [2008]. These fluxes take into account source terms. We also show that these schemes are entropic and preserve the admissible states, E r ą 0 and ||Fr|| cEr ď 1. Standard numerical tests show the good behavior of the schemes in both freestreaming and diffusive limits.

Nonlinear Jacobi method and geometric multigrid

In this chapter, we present an additional step towards the development of an asymptotic preserving and time-implicit solver that preserves the admissible states. We discretize the M 1 model with a time-implicit HLL solver without source terms. Instead of solving the nonlinear system with the Newton-Raphson method, we use the Jacobi method presented in [START_REF] Pichard | Martin Mathematiques appliquees et calcul scientifique Bordeaux[END_REF]. We show that the solution built with this iterative process preserves the admissible states.

However, this method can be slow to converge when the resolution increases. To tackle this issue, we use a geometric multigrid algorithm [START_REF] Briggs | A Multigrid Tutorial[END_REF]]. The problem is solved on a coarser grid, where the Jacobi method converges faster. Introducing a pseudo-time [START_REF] Kifonidis | On multigrid solution of the implicit equations of hydrodynamics -experiments for the compressible euler equations in general coordinates[END_REF] allows us to preserve the admissible states.

First numerical tests in the free-streaming regime show the efficiency of the method.

High performance computing and linear algebra

The scheme developed in chapter 2 requires solving large sparse ill-conditioner linear systems. We use the library Trilinos to do it. We present the different packages we used and we show how they allow us to target different architectures, such as multicores, many-cores, and GP-GPUs. We also show some performance results and we discuss the impact of the preconditioner on a physical application.

Publications and communications

The work presented in this thesis led to the following publications and oral communications:

-H. Bloch, P. Furthermore, the work presented in chapters 3 and 4 is the subject of articles in preparation. This chapter is dedicated to the presentation of radiative transfer equation (section 1.1). The quantity of interest is the specific intensity, which describes the rate of radiative transfer of energy, at a point x and time t. It is used in the study of processes involving propagation of radiation and its interaction with the surrounding medium. It relies on classical electromagnetic radiation, whereas it is conceptually different from the description of Maxwell electromagnetic fields.

Bibliography

Due to 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, if coupled with hydrodynamics, direct simulations are out of reach for modern computers. Different models have been developed to reduce the computational cost. We will focus on the moment models [Levermore, 1984] in section 1.2: the specific intensity is averaged over the direction of propagation of photons. It presents several advantages: the computational cost is lower than other methods such as Monte-Carlo method and, mostly, it is easy to couple it with a grid-based hydrodynamics code. Furthermore, we derive a specific closure relation, the M 1 model, by maximizing the radiative entropy (section 1.3). This closure relation allows us to treat accurately both free-streaming and diffusive regimes. In particular, we investigate the diffusive limit in section 1.4. Finally, we enhance this model for a moving fluid in section 1.5 to derive the radiation hydrodynamics equations.

Radiative transfer equation

In the vacuum, photons propagate in a straight line at velocity cΩ, where c is the speed of light and Ω is the direction of propagation of photons. At a time t, a photon is determined by its position x, its direction of propagation Ω, and its frequency ν.

To derive the radiative transfer equation, we follow Chandrasekhar [1960]. Let us define the specific intensity I such that the flux of radiative energy Φ crossing an oriented elementary surface dS around the point x in the elementary solid angle dΩ around the direction Ω with frequency in rν, ν `dνs, at time t, is Φ px, t, Ω, νq " I px, t, Ω, νq Ω ¨dSdΩdν.

(1.1) by the gas or by another source. Photons could be scattered, which will change their direction of propagation and perhaps their frequency. Let us now focus on each of these phenomena.

We first investigate absorption. We consider the flux of radiative energy propagating in an absorbing medium. A part dΦ a of the radiative energy is absorbed:

dΦ a px `dl, t, Ω, νq " ´σa ν Φ px, t, Ω, νq dl, (1.2)
with σ a ν the absorption coefficient. Its inverse, 1 σ a ν , is the absorption mean free path. We now look at thermal emission. A part

dΦ e px `dl, t, Ω, νq " η ν Ω ¨dSdΩdνdl, (1.3)
with η ν the emission coefficient, is emitted by the surrounding matter. From Chandrasekhar [1960], at local thermodynamic equilibrium (LTE),

η ν " σ a ν B pν, T g q " σ a ν 2hν 3 c 1 exp ´hν k b Tg ¯´1 , (1.4)
where B is the black body specific intensity, h is the Planck constant, k b is the Boltzmann constant and T g is the temperature of the matter. From now on, we always assume LTE. At global thermodynamic equilibrium, Φ does not depend on x, therefore 0 " BΦ Bl " dΦ e `dΦ a , (1.5) which leads to I px, t, Ω, νq " B pν, T g q .

(1.6)

We will recover eq. 1.6 in section 1.B. Finally, let us consider scattering. On one hand, part of the energy propagating initially in the direction Ω is scattered into a different direction Ω 1 . This will decrease the flux of radiative energy in the direction Ω. By writing σ s ν the scattering coefficient, the energy lost by scattering is ´σs ν Φ px, t, Ω, νq dl.

(1.7)

The scattering mean free path is 1 σ s ν . On the other hand, energy propagating initially in a direction Ω 1 , different from Ω, can be scattered and then is propagated in the direction Ω. This will increase the flux of radiative energy in the direction Ω. Let us introduce the phase function p ν px, Ω, Ω 1 q which represents the part of the flux of radiative energy that is scattered from the direction Ω 1 into the direction Ω at the point x. Then, the gain is

σ s ν ż S 2 p ν px, Ω, Ω 1 q Φ px, t, Ω 1 , νq dΩ 1 dl.
(1.8)

If the medium is homogeneous, which we assume from now on, then p ν does not depend on the position x. Furthermore, we assume that p ν depends only on the angle between Ω and Ω 1 . Therefore, we write p ν pΩ ¨Ω1 q. With our choice of normalization, we have

ż S 2
p ν pΩ ¨Ω1 q dΩ " 1, (1.9) but other choices are possible.

The difference of energy in the elementary volume between the moment t and the moment t `dt is then

pI px `dl, t `dt, Ω, νq ´I px, t, Ω, νqq Ω ¨dSdΩdν " ´pσ a ν `σs ν q I px, t, Ω, νq Ω ¨dSdΩdνdl `σa ν B pν, T g q Ω ¨dSdΩdνdl `σs ν ż S 2 p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 Ω ¨dSdΩdνdl.
(1.10)

Let us notice that dt " dl c . Then, using a Taylor expansion at first order with dl that goes to 0, one has

I px `dl, t `dt, Ω, νq " I px, t, Ω, νq `dl ˆBI Bl `1 c BI Bt ˙px, t, Ω, νq `opdlq. (1.11)
Furthermore, BI Bl " Ω ¨∇I and

ˆ1 c B t `Ω ¨∇˙I px, t, Ω, νq " ´pσ a ν `σs ν q I px, t, Ω, νq `σa ν B pν, T g q `σs ν ż S 2 p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 .
(1.12) Equation 1.12 is an integro-differential equation that describes the transport of the specific intensity. One can notice the analogy between eq. 1.12 and the Boltzmann equation that describes the statistical behavior of a thermodynamic system out of equilibrium.

Two main regimes can arise, depending on the mean free path of photons compared to the characteristic length of the system, written L [Mihalas and 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, 1 σ a ν `σs ν ! L), 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, 1

σ a ν `σs ν " L).
Because solving eq. 1.12 within the scope of radiation hydrodynamics will be costly, we present now the moment method used in this work.

Moment model

One can consider only the moment of order 0 (the radiative energy), leading to the flux-limited diffusion (FLD) approximation [START_REF] Levermore | A flux-limited diffusion theory[END_REF]. Because this model considers only the moment of order 0, its computational cost is quite low, but it is 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 1 closure relation [Dubroca and Feugeas, 1999]. However, this method can suffer from artifacts when multiple beams cross in the free-streaming regime [START_REF] González | Contribution to the numerical study of radiation hydrodynamics: from radiative shocks experiments to astrophysical jets[END_REF]. One can solve this issue by using a three-moment model (radiative energy, radiative flux, and radiative pressure) with the M 2 closure relation [START_REF] Pichard | The m 2 model for dose simulation in radiation therapy[END_REF]. 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 1 closure relation because the computational cost remains affordable, and we do not encounter in our applications the problem of beams crossing in the free-streaming regime.

We limit ourselves to gray radiative transfer, but the following can be generalized to multigroup radiative transfer (e.g., [START_REF] Turpault | A consistent multigroup model for radiative transfer and its underlying mean opacities[END_REF].

Let us consider the three first moments of the specific intensity: the gray radiative energy E r , the gray radiative flux F r and the gray radiative pressure P r defined as

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % E r " 1 c ż 8 0 ż S 2 I px, t, Ω, νq dΩdν F r " ż 8 0 ż S 2 I px, t, Ω, νq ΩdΩdν P r " 1 c ż 8 0 ż S 2
I px, t, Ω, νq Ω b ΩdΩdν.

( dΩdν " cE r .

(1.14)

This condition rewrites f ď 1, where f " ||f || " ||Fr|| cEr is the reduced flux. It ensures that the radiative energy cannot be transported faster than the speed of light.

The mean over solid angles and frequency of eq. 1.12 and its product by Ω give the following system:

$ ' & ' % B t E r `∇ ¨Fr " cσ a `ar T 4 g ´Er Bt F r `c2 ∇ ¨Pr " ´cσ s F r .
(1.15a)

(1.15b) The computation of the source terms is made explicit in section 1.A. σ a and σ s are weighted means of σ a ν and σ s ν respectively. More details can be found in section 1.A and [START_REF] González | Contribution to the numerical study of radiation hydrodynamics: from radiative shocks experiments to astrophysical jets[END_REF].

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 B t pρc v T g q " ´cσ a `ar T 4 g ´Er ˘.

(1.16) ρc v T g is the gas internal energy, with ρ the density of the fluid and c v the heat capacity, defined by c v "

k b
µm H pγ´1q for a perfect gas, where µ is the mean molecular weight, m H is the mass of hydrogen and γ is the adiabatic index of the gas.

In the following section, we specify the closure relation, which is a way to express P r as a function of E r and F r .

Closure relation

We have chosen here to use the M 1 closure relation for its good properties. To derive it, Levermore [1984] applies a Lorentz transform to an isotropic distribution of photons. The same relation can also be obtained by maximizing the radiative entropy [Dubroca and Feugeas, 1999;[START_REF] Buet | Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics[END_REF]. We present the latter because some results will be used in chapter 3.

We only give here the main idea of the derivation of the model, detailed computations can be found in section 1.B.

From [START_REF] Buet | Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics[END_REF], we first define the gray radiative entropy as

S r " ´2k b c 3 ż 8 0 ż S 2
ν 2 pn log n ´pn `1q log pn `1qq dΩdν, (1.17)

where n " c 2 2h I ν 3 is the photon occupation number at the frequency ν. At a given time t, in a unit volume centered in x of volume d 3 x, there is ş 8 0 ş S 2 n px, t, Ω, νq dΩdνd 3 x photons. S r is the entropy of a photon gas. The derivation of eq. 1.17 relies on Gibbs statistical mechanics (e.g., ter Haar 1995) and is beyond the scope of this work.

S r is the physical entropy, it reaches a maximum. One can also define the mathematical entropy, with the opposite sign (e.g., Dubroca and Feugeas 1999). In that case, the radiative entropy reaches a minimum.

One can show that the radiative entropy S r is a strictly concave function of the specific intensity I. We write n the occupation number that maximizes the radiative entropy. Using the Lagrange multiplier method (e.g., [START_REF] Boyd | [END_REF]), there exists θ r and β such that

n " 1 exp ´hν k b θr `1 ´β¨Ω c ˘¯´1
.

(1.18) θ r is homogeneous to a temperature, and we have θ r ě 0. Let us emphasize that θ r is not the gas temperature T g . Similarly, one can show that β is homogeneous to a velocity and β " ||β|| ď c. Let us notice that if β " 0, then I pnq is the black body specific intensity at temperature θ r . After computation (see section 1.B), one has

$ ' ' ' ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' ' ' ' % E r " c 4 3c 2 `β2 3 pc 2 ´β2 q 3 a r θ 4 r F r " 4 3 c 6 a r θ 4 r pc 2 ´β2 q 3 β " pq `Er q β P r " β b F r c 2 `qI " ˆ1 ´χ 2 I `3χ ´1 2 n b n ˙Er , (1.19)
with q " c 2 ´β2 3c 2 `β2 E r , I the identity matrix, and n " f f is a unit vector aligned with the radiative flux and χ " 3`4f 2 5`2 a 4´3f 2 is the Eddington factor [Levermore, 1984].

Using eq. 1.19 to express P r as a function of E r and F r is the closure relation of the M 1 model. It preserves both free-streaming and diffusive limits. On one hand, if f " 1, then P r " E r n b n, only the transport regime remains. On the other hand, if f " 0, the model in the diffusive regime simplifies into the P 1 model, with P r " 1 3 E r I. The radiative pressure tensor becomes diagonal and isotropic. We will look into the diffusion limit more precisely in section 1. 4.

Finally, one can also show that the maximum radiative entropy obeys the following hyperbolic equation

B t S r `∇ ¨pβS r q " cσ a θ r `ar T 4 g ´Er ˘`σ s cθ r β ¨Fr . (1.20)
This equation will be used in chapter 3. It allows us to prove that the numerical scheme developed there preserves the admissible states E r ą 0 and f ď 1.

Diffusive limit

Let us now focus on the diffusive regime with the hydrodynamics frozen, i.e., the limit of large opacity and long timescale. From [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF], when the grid does not sample the mean free path of the photons, the numerical solution is dominated by the numerical diffusion. In order to understand this phenomena, we explore here the diffusive limit at the continuous level. Computations done here will lead to the development of so-called asymptotic preserving schemes in chapters 2 and 3. These schemes will capture the behavior described in the section, independently of the resolution.

Following [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF], we introduce a rescaling parameter ε to write the time (resp. the opacity) as t " εt (resp. σa " εσ a and σs " εσ s ). The radiative energy, the radiative flux, and the gas temperature are expanded with ε, e.g., E r " E r,0 `εE r,1 `O pε 2 q. System 1.15 and eq. 1.16 become

$ ' ' & ' ' % ε 2 B r t E r `ε∇ ¨Fr " c r σ a `ar T 4 g ´Er ε2 B r t F r `εc 2 ∇ ¨Pr " ´c r σ s F r ε 2 B r t pρc v T g q " ´c r σ a `ar T 4 g ´Er ˘.
(1.21a)

(1.21b)

(1.21c) By expanding eqs. 1.21a and 1.21b at order 0, we have

# E r,0 " a r T 4 g,0 F r,0 " 0. (1.22)
Because F r,0 " 0, let us consider the P 1 closure relation instead of the M 1 closure relation. We use P r " 1 3 E r I instead of eq. 1.19. This model is known to have the correct behavior only in the diffusive regime. It does not capture properly the free-streaming regime and should be used carefully. However, we use it here for its simplicity.

Expanding eq. 1.21b at order 1 leads to F r,1 " ´c 3 r σ s ∇E r,0 .

(1.23)

Finally, expanding the sum of eqs. 1.21a and 1.21c at order 2 gives B r t pE r,0 `ρc v T g,0 q ´∇ ˆc 3 r σ s ∇E r,0 ˙" 0.

(1.24)

In the diffusive limit, the total energy E r `ρc v T g at order 0 obeys the diffusion equation given by eq. 1.24. Chapters 2 and 3 present our design of two asymptotic preserving schemes, that is, numerical schemes that will degenerate to the discretization of eq. 1.24 in the diffusive regime.

In this section, we only have considered the case where the hydrodynamics is frozen. In the next one, we present the radiation hydrodynamics equations, which take into account the evolution of the fluid.

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 the hydrodynamics is no longer frozen and both photons and the gas are moving, 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 non-conservative 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: S Fr puq " σ s u ¨Pr `σa a r T 4 g u `2σ a c 2 u ¨pu ¨Fr q `pσ s ´σa q E r u,

$ ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' % B t ρ `∇ ¨pρuq " 0 B t pρuq
(1.26) at first order in u c . To close the system, we also add the equation of state of an ideal gas: p " ρepγ ´1q.

An asymptotic development for the radiation hydrodynamics case is presented in section 1.C.

Bibliography

1.A Computation of source terms

The mean over solid angles and frequency of eq. 1.12 and its product by Ω give the following system:

$ & % B t E r `∇ ¨Fr " S E 1 c B t F r `c∇ ¨Pr " S F , (1.27) where $ ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' % S E " ´ż 8 0 ż S 2 pσ a ν `σs ν q I px, t, Ω, νq dΩdν `ż 8 0 ż S 2 σ a ν B pν, T g q dΩdν `ż 8 0 ż S 2 ˆσs ν ż S 2 p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 ˙dΩdν S F " ´ż 8 0 ż S 2 pσ a ν `σs ν q I px, t, Ω, νq ΩdΩdν `ż 8 0 ż S 2 σ a ν B pν, T g q ΩdΩdν `ż 8 0 ż S 2 ˆσs ν ż S 2
p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 ˙ΩdΩdν

(1.28) are the source terms. We now make them explicit. Let us consider first the term for the radiative energy equation. We compute each term individually.

First,

ż 8 0 ż S 2 σ s ν ż S 2 p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 dΩdν " ż 8 0 σ s ν ż S 2 ˆżS 2 p ν pΩ ¨Ω1 q dΩ loooooooooooomoooooooooooon "1,@Ω 1 I px, t, Ω 1 , νq dΩ 1 dν " ż 8 0 ż S 2
σ s ν I px, t, Ω, νq dΩdν.

( B pΩ, ν, T g q ΩdΩ " 0.

(1.37)

Finally, S F " ´pσ F `p1 ´gq σ R q F r .

(1.38)

The means of opacity σ E , σ P , σ F and σ R can be connected with Rosseland and Planck mean opacities. See [START_REF] González | Contribution to the numerical study of radiation hydrodynamics: from radiative shocks experiments to astrophysical jets[END_REF] for more details. In the following, we always make the approximation σ E " σ P and we write σ E " σ P " σ a . Furthermore, we note σ s " σ F `p1 ´gq σ R to write eq. 1.27 as

# B t E r `∇ ¨Fr " cσ a `ar T 4 g ´Er Bt F r `c2 ∇ ¨Pr " ´cσ s F r .
(1.39)

1.B Maximization of radiative entropy 1.B.1 Radiative entropy

Let us recall that we have defined the radiative entropy as

S r " ´2k b c 3 ż 8 0 ż S 2
ν 2 pn log n ´pn `1q log pn `1qq px, t, Ω, νq dΩdν.

(1.40)

We also define the entropy flux as

Q r " ´2k b c 2 ż 8 0 ż S 2
ν 2 pn log n ´pn `1q log pn `1qq px, t, Ω, νq ΩdΩdν.

(1.41)

Let us first compute the equation of evolution for the radiative entropy. From now on, we assume that eq. 1.12 has a solution smooth enough. This hypothesis allows us to compute the derivatives without further justification. Then, because n "

" k b h ż 8 0 ż S 2 1 ν log ˆn n `1 ˙px, t, Ω, νq pσ a ν `σs ν q I px, t, Ω, νq dΩdν ´kb h ż 8 0 ż S 2 1 ν log ˆn n `1 ˙px, t, Ω, νq σ a ν B pν, T g q dΩdν ´kb h ż 8 0 ż S 2 1 ν log ˆn n `1 ˙px, t, Ω, νq σ s ν ż S 2
p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 dΩdν.

(1.44) The radiative entropy obeys a hyperbolic equation. The source terms will be made explicit in section 1.B.4.

1.B.2 Maximization of entropy

Let us now compute n that maximizes the radiative entropy. We first show that the radiative entropy S r is a strictly concave as a function of the specific intensity I. Let us recall that a function is strictly concave if and only if its second derivative is strictly negative.

After some computations similar to what has been done in the previous section, one can show that

$ ' ' & ' ' % BS r BI " ´kb ch ż 8 0 ż S 2 1 ν log ˆn n `1 ˙px, t, Ω, νq dΩdν B 2 S r BI 2 " ´ck b 2h 2 ż 8 0 ż S 2 1 ν 4
1 pn pn `1qq px, t, Ω, νq dΩdν.

(1.45)

The second equation shows that S r is strictly concave. Because max S r pnq " ´min p´S r pnqq, computing the maximum of radiative entropy leads to solve the convex optimization problem ´min p´S r pnqq under the constraints ´Er pnq ď 0 and ||F r pnq||´cE r pnq ď 0.

Using the Karush-Kuhn-Tucker (KKT) conditions (e.g., Boyd et Let us introduce two variables, θ r P R and β P R 3 defined as

$ ' ' & ' ' % θ r " 1 λ 1 `cλ 2 β " c 2 θ r λ 2 F r pnq ||F r pnq|| .
(1.47)

θ r is obviously non-negative and homogeneous to a temperature. Some easy computations show that β is homogeneous to a velocity and β " ||β|| ď c:

λ 1 ě 0 ô β " c 2 λ 2 λ 1 `cλ 2 ď c 2 λ 2 cλ 2 " c.
( n given by eq. 1.52 maximizes the radiative entropy S r . At global thermodynamic equilibrium, β " 0 and the specific intensity given by eq. 1.52 is the black body specific intensity (eq. 1.6), by identifying θ r with the gas temperature, T g .

1.B.3 M 1 closure relation

To derive P r as a function of E r and F r , we express the three moments of the specific intensity as functions of θ r and β.

1.B.3.1 Radiative energy

Let us start with the radiative energy. Using

I " 2hν 3 n c 2 " 2h c 2 ν 3 exp `´ν α p1`A¨Ωq ˘´1 with α " k b θr h and A " ´β c , one has E r " 2h c 3 ż 8 0 ż S 2 ν 3 exp `´ν α p1 `A ¨Ωq ˘´1 dΩdν. (1.53)
Using the change of variable r " ´ν α p1 `A ¨Ωq, one has

E r " 2hα 4 c 3 ż 8 0 r 3 exp prq ´1 dr ż S 2 1
p1 `A ¨Ωq 4 dΩ.

(1.54)

The integral over the frequency has already been made explicit in section 1.A. Let us now focus on the integral over the unit sphere. We choose a reference frame such that the z-axis is parallel to the vector A and θ is the angle between A and Ω (see fig. 1.3). Therefore, A ¨Ω " A ||Ω|| lo omo on

"1 cos θ, (1.55)
with A " ||A||. Furthermore, we write the units vector e x , e y and e z . In particular, e z " A A . We can now compute the integral over the unit sphere:

ż S 2 1 p1 `A ¨Ωq 4 dΩ " ż 2π 0 dϕ ż π 0 sin θ p1 `A cos θq 4 dθ " 4π p3 `A2 q 3 p1 ´A2 q 3 .
(1.56)

Finally, E r " c 4 3c 2 `β2 3 pc 2 ´β2 q 3 a r θ 4 r .

(1.57)

1.B.3.2 Radiative flux

Likewise, one can compute F r :

F r " 2h c 2 ż 8 0 ż S 2 ν 3 exp `´ν α p1 `A ¨Ωq ˘´1 ΩdΩdν " 2hα 4 c 2 π 4 15 ż S 2 Ω p1 `A ¨Ωq 4 dΩ.
(1.58)

Again, let us focus on the integral over the sphere. Using spherical coordinates, we can write Ω " ¨sin θ cos ϕ sin θ sin ϕ cos θ ' .

(1.59) Therefore, the x-coordinate of ş (1.60)

Likewise, the y-coordinate of the integral is also 0. The z-coordinate is

ż 2π 0 dϕ ż π 0 cos θ sin θ p1 `A cos θq 4 dθ " ´16π A 3 p1 ´A2 q 3 .
(1.61)

Finally,

F r " ´π4 15
2hα 4 c 2 16π A 3 p1 ´A2 q 3 e z " 4 3
c 6 a r θ 4 r pc 2 ´β2 q 3 β " pq `Er q β, (1.62) with q " c 2 ´β2 3c 2 `β2 E r " 1´χ 2 E r
, where χ is the Eddington factor, defined by χ " 3`4f 2 5`2 a 4´3f 2 . The last expression for F r will be used in chapter 3.

Equations 1.57 and 1.62 are equivalent to

β " c f 3χ ´1 2 f θ r " 2 f ´´1 `a4 ´3f 2 ¯1 4 b f 2 ´2 `a4 ´3f 2 ˆEr a r ˙1 4 .
(1.63)

Finally, using I " hν 3 n c 2 and eq. 1.52, the distribution of photons associated to the M 1 model is

2hν 3 c 2 1 exp ´hν k b θr ´1 ´2´a4´3f 2 f 2 f ¨Ω¯¯.
(1.64)

1.B.3.3 Radiative pressure

Finally,

P r " 2h c 3 ż 8 0 ż S 2 ν 3 exp `´ν α p1 `A ¨Ωq ˘´1 Ω b ΩdΩdν " 2hα 4 c 3 π 4 15 ż S 2 Ω b Ω p1 `A ¨Ωq 4 dΩ.
(1.65)

As before, we focus on the integral over the unit sphere. Ω b Ω is a symmetric 3 ˆ3 matrix and, using eq. 1.59, we have

Ω b Ω " ¨sin 2 θ cos 2 ϕ
sin 2 θ cos ϕ sin ϕ sin θ cos θ cos ϕ sin 2 θ cos ϕ sin ϕ sin 2 θ sin 2 ϕ sin θ cos θ sin ϕ sin θ cos θ cos ϕ sin θ cos θ cos ϕ cos 2 θ ' .

(1.66)

Because P r is a tensor, we write its components P i,j r for i, j P tx, y, zu. With this notation, P x,y r " P y,x r " P x,z r " P z,x r " P y,z r " P z,y r " 0. One can easily compute individually the non-zero coefficients:

P x,x r " P y,y r " 2hα 4 c 3 π 4 15 4π 3 p1 ´A2 q 2 P z,z r " 2hα 4 c 3 π 4 15 4π p1 `3A 2 q 3 p1 ´A2 q 3 " 2hα 4 c 3 π 4 15 ˆ4π 3 p1 ´A2 q 2 `4π ˆ4A 2 3 p1 ´A2 q 3
˙.

(1.67) Therefore, P r " c 4 a r θ 4 r 3 pc 2 ´β2 q 2 pe x b e x `ey b e y `ez b e z q `4c 4 a r θ 4 r β 2 3 pc 2 ´β2 q 3 e z b e z " c 4 a r θ 4 r 3 pc 2 ´β2 q 2 I `4c 4 a r θ 4 r β 2 3 pc 2 ´β2 q 3 e z b e z .

(1.68)

Nevertheless,

e z " A A " ´β β " ´3c 2 `β2 4cβ f . (1.69) Finally, P r " ˆc2 ´β2 3c 2 `β2 I `3c 2 `β2 4c 2 f b f ˙Er " β b F r c 2 `qI.
(1.70)

One can also check that eq. 1.70 is equivalent to P r " DE r , where D is the Eddington tensor, defined by D " 1´χ 2 I `3χ´1 2 n b n, with n " f f a unit vector aligned with the radiative flux [Levermore, 1984].

1.B.4 Back to the radiative entropy

In the previous section, we have expressed the moments of the specific intensity as functions of θ r and β. Following Buet and Despres [2004, Appendix B], we do the same for the maximum of radiative entropy. Using eq. 1.17, the maximum of entropy is given by

S r " ´2k b c 3 ż 8 0 ż S 2
ν 2 pn log n ´pn `1q log pn `1qq dΩdν.

(1.71)

Using the change of variable r " ´ν α p1 `Ω ¨Aq, with M α "

S r " ´8πk 4 b θ 3 r c h 3 pc 2 ´β2 q 2 ż S 2
1 exp `r α ˘´1 to write S r " ´8πk 4 b θ 3 r c h 3 pc 2 ´β2 q 2 hp1q.
(1.74)

We now have to compute hpαq. Let us notice that Mα Mα`1 " exp `ξ α ˘. We then have

h 1 pαq " ´α2 ż 8 0 ξ 3 exp pξq pexp pξq ´1q
dξ.

(1.75)

By integrating by parts, we recognize the Bose-Einstein function:

h 1 pαq " ´α2 ż 8 0 4ξ 3 exp pξq ´1 dξ " ´4π 4 α 2 15 . (1.76)
Because the entropy is defined up to a constant, integrating this equation leads to

hpαq " ´4π 4 α 3 45 . (1.77) Finally, S r " 4c 4 a r θ 3 r 3 pc 2 ´β2 q 2 .
(1.78)

Likewise, we compute the entropy flux, defined by eq. 1.41:

Q r " ´2k b c 2 ż 8 0 ż S 2 ν 2 Ω pn log n ´pn `1q log pn `1qq dΩdν " ´2k 4 b θ 3 r c 2 h 3 gp1q ż S 2 Ω p1 `Ω ¨Aq 3 dΩ " 4c 4 a r θ 3 r 3 pc 2 ´β2 q 2 β " S r β.
(1.79) Equation 1.44 simplifies into

B t S r `∇ ¨pS r βq " k b h ż 8 0 ż S 2 1 ν log ˆn n `1˙p σ a ν `σs ν q I px, t, Ω, νq dΩdν ´kb h ż 8 0 ż S 2 1 ν log ˆn n `1 ˙σa ν B pν, T g q dΩdν ´kb h ż 8 0 ż S 2 1 ν log ˆn n `1 ˙σs ν ż S 2
p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 dΩdν.

(1.80) We can now specify the source term of eq. 1.80. Let us recall that

k b hν log `n n`1 ˘" β¨Ω cθr ´1 θr . Therefore, B t S r `∇ ¨pS r βq " ż 8 0 ż S 2 ˆβ ¨Ω cθ r ´1 θ r ˙pσ a ν `σs ν q I px, t, Ω, νq dΩdν ´ż 8 0 ż S 2 ˆβ ¨Ω cθ r ´1 θ r ˙σa ν B pν, T g q dΩdν ´ż 8 0 ż S 2 ˆβ ¨Ω cθ r ´1 θ r ˙σs ν ż S 2
p ν pΩ ¨Ω1 q I px, t, Ω 1 , νq dΩ 1 dΩdν.

(1.81) Because θ r and β do not depend on Ω and ν, we recognize computations done in section 1.A and

B t S r `∇ ¨pS r βq " ´β ¨SF cθ r `SE θ r " cσ a θ r `ar T 4 g ´Er ˘`σ s cθ r β ¨Fr .
(1.82)

1.C Diffusive limit for radiation hydrodynamics

As in section 1.4, we consider the diffusive limit with the P 1 closure relation, but with hydrodynamics. We introduce a rescaling parameter ε to write the time (resp. the opacity) as t " εt (resp. σa " εσ a and σs " εσ s ). Because the velocity of the fluid is smaller than the speed of light ( u c ! 1), we also rescale it as ũ " u ε . Let us focus on the equations describing the evolution of the radiative variables:

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % ε 2 B r t E r `ε∇ ¨Fr " c r σ a `ar T 4 g ´Er ˘`ε 2 r σ a ´r σ s c r u ¨Fr `ε2 4 3 r σ s c E r r u ¨r u ε 2 B r t F r `ε c 2 3 ∇E r " ´r σ s cF r `εc ˆ4 3 r σ s ´r σ a ˙r uE r `εc r σ a r ua r T 4 g `ε2 2 r σ a c r u ¨pr u ¨Fr q.
(1.83a)

(1.83b) By expanding eqs. 1.83a and 1.83b at order 0, we have

# E r,0 " a r T 4 g,0
F r,0 " 0.

(1.84)

Expanding eq. 1.83b at order 1 leads to

F r,1 " ´c 3 r σ s ∇E r,0 `4 3 E r,0 Ă u 0 . (1.85)
Finally, looking at the radiative energy and the gas internal energy at order 2, source terms cancel each other, only the divergence of the radiative flux at order 1 remains, and we have

B r t pρc v T g,0 `Er,0 q ´∇ ¨ˆc 3 r σ s ∇E r,0 ˙" ´4 3 ∇ pE r,0 Ă u 0 q . (1.86)
One can also look at eq. 1.83a at order 2. Expanding eq. 1.83a at order 2 gives

B r t E r,0 ´∇ ˆc 3 r σ s ∇E r,0 ˙" c r σ a `6a r T 2 g,0 T 2 g,1 `4a r T 3 g,0 T g,2 ´Er,0 ˘´4 3 ∇ ¨pE r,0 Ă u 0 q ´2σ a ´σs 3σ s Ă u 0 ¨∇E r,0 `2 4 3 σ a c E r,0 Ă u 0 2 .
(1.87) We recover Eq. 43 of Krumholz et al. [2007]. The term c r σ a `6a r T 2 g,0 T 2 g,1 `4a r T 3 g,0 T g,2 ´Er,2 ȋs the development at second order of the term κ 0 p4πB ´cEq. Because we do not neglect any terms O `ũ c ˘, some coefficients are slightly different. Krumholz et al.

[2007] discuss the importance of the term 4 3

Ă σ a c E r,0 Ă u 0 2 .
It is interpreted as a "relativistic work", as opposed to the term of "diffusion work", Ă u 0 ¨∇E r,0 . Therefore, it is important for non-equilibrium non-uniform dynamic diffusion problem, whereas it is not important for static diffusion problems.

This chapter is the adaptation of an article published in Astronomy&Astrophysics, see Bloch et al. 2021.

Even though the M 1 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, [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF] 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, to obtain a simpler solver, 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.

This chapter is organized as follows. We go through our new numerical scheme, well-balanced and asymptotic preserving in the diffusive limit in section 2.1. In section 2.2, we present some numerical test cases to show the importance of the asymptotic preserving and well-balanced properties. We also present a physical application about the stability of the ionization front in an H ii region in a massive pre-stellar dense core. Finally, we reach our conclusion and discuss the limitations of the scheme in section 2.3.

Numerical scheme and algorithm

Radiation transport in a static fluid

Let us first introduce some notations: we note ∆x the step along the x-direction. ∆t is the time interval between the current time t n and t n`1 . We write x i the center of the cell i and x i`1 2 the interface between the cell i and the cell i `1. We use the notation u n i to represent the averaged quantity associated with the field u at time t n in the cell i (finite volume). Finally, we note u n i`1 2 to represent the quantity associated with the field u at time t n 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.

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 of the fluid for the hydrodynamics. Because we are interested in radiation hydrodynamics, we will consider a long timescale for the radiative transfer. Therefore, we use a time-implicit integration for the radiative transfer. A similar development can be done with a semi-implicit solver: source terms remain implicit, but the hyperbolic part is time-explicit.

Hyperbolic system

Following González et al. [2007], we discretize the hyperbolic part of eq. 1.15 using a first-order Godunov type solver [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. From Berthon and Turpault [2011], we also introduce an extra parameter α which will be specified in section 2.1.1.3:

$ ' ' ' ' ' & ' ' ' ' ' % E n`1 i " E n i ´∆t ∆x ´αi`1 2 F i`1 2 ´αi´1 2 F i´1 2 ¯`cσ a i ∆t ´ar `T n`1 i ˘4 ´En`1 i F n`1 i " F n i ´∆t ∆x ´Pi `1 2 ´Pi ´1 2 ¯´c∆ttσ s F r u n`1 i ρc v T n`1 i " ρc v T n i ´cσ a i ∆t ´ar `T n`1 i ˘4 ´En`1 i ¯, (2.1)
where F i`1 2 and P i`1 2 are the numerical fluxes given by

F i`1 2 " λ ì`1 2 F n`1 i ´λí `1 2 F n`1 i`1 λ ì`1 2 ´λí `1 2 `λì `1 2 λ í`1 2 λ ì`1 2 ´λí `1 2 `En`1 i`1 ´En`1 i Pi `1 2 "c 2 λ ì`1 2 P n`1 i ´λí `1 2 P n`1 i`1 λ ì`1 2 ´λí `1 2 `λì `1 2 λ í`1 2 λ ì`1 2 ´λí `1 2 `F n`1 i`1 ´F n`1 i ˘, (2.2)
with λ ì`1 2 " max p0, λ max q and λ í`1 2 " min p0, λ min q, where λ max and λ min are the eigenvalues of eq. 1.15. From [START_REF] Berthon | An hllc scheme to solve the m1 model of radiative transfer in two space dimensions[END_REF], we have

λ max,min " c ¨fx ξ ˘?2 b pξ ´1qpξ `2qp2pξ ´1qpξ `2q `3f 2 y q ? 3ξpξ `2q ' , (2.3)
with ξ " a 4 ´3f 2 . See Fig. 1 of González et al. [2007] for more details about the structure of the eigenvalues. tσ s F r u n`1 i is a well-chosen discretization of the term σ s F r in the cell i and at time t n`1 , and is specified in the next section.

Well-balanced modification of the source term

From [START_REF] Berthon | A Well-Balanced Finite Volume Scheme for a Mixed Hyperbolic/Parabolic System to Model Chemotaxis[END_REF], a well-balanced scheme catches the correct steady regime. The steady state, if it exists, is given by

$ ' & ' % E r " a r T 4 g ∇ ¨Fr " 0 c∇ ¨Pr " ´σs F r . (2.4a) (2.4b) (2.4c) Equation 2.4c is discretized by c 2 ´p∇ ¨P q n`1 i`1 2 ´p∇ ¨P q n`1 i´1 2 ¯" ´tσ s F u n`1 i , with p∇ ¨P q n`1 i`1 2 " P n`1 i`1 ´P n`1 i ∆x . An obvious choice for tσ s F u n i is tσ s F u n`1 i " σ s i F n`1 i .
(2.5) However, using this formulation, eq. 2.4c is discretized as

´c 2 ´p∇ ¨P q n`1 i`1 2 `p∇ ¨P q n`1 i´1 2 ¯" σ s i F n`1 i . (2.6)
The radiative flux remains cell-centered and is equal to the divergence of radiative pressure, 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 section 2.2.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 of this balance with the balance between radiative pressure and radiative flux source term in eq. 2.4c), we choose to use an average of a face discretization of the radiative flux source term:

tσ s F u n`1 i " 1 2 ´σs i`1 2 F n`1 i`1 2 `σs i´1 2 F n`1 i´1 2 ¯, (2.7) with $ ' & ' % σ s i`1 2 " 1 2
`σs i `σs i`1 F n`1 i`1 2 " 1 2 `F n`1 i `F n`1 i`1 ˘.
(2.8)

One way to interpret this equation is to remember that

tσ s F u n`1 i " 1 ∆x ż x i`1 2 x i´1 2 σ s pxqF `tn`1 , x ˘dx.
(2.9) Equation 2.5 is obtained with the rectangle rule for numerical integration of eq. 2.9:

tσ s F u n`1 i " x i`1 2 ´xi´1 2 ∆x ˆσs ˆxi´1 2 `xi`1 2 2 ˙F ˆtn`1 , x i´1 2 `xi`1 2 2 ˙" σ s px i qF pt n`1 , x i q " σ s i F n`1 i , (2.10)
whereas eq. 2.7 is given by the trapezoidal rule:

tσ s F u n`1 i " x i`1 2 ´xi´1 2 2∆x ´σs ´xi´1 2 ¯F ´tn`1 , x i´1 2 ¯`σ s ´xi`1 2 ¯F ´tn`1 , x i`1 2 ¯" 1 2 ˆσs i´1 2 F n`1 i´1 2 `σs i `1 2 F n`1 i`1 2 ˙.
(2.11) To have

σ s i`1 2 F n`1 i`1 2 " ´cP n`1 i`1 ´P n`1 i ∆x (2.12)
in the whole domain, we also impose it as boundary condition:

σ s 1 2 F n`1 1 2 " ´c P n`1 1 ´P n`1 0 ∆x , (2.13)
where P n`1 0 is the radiative pressure given by the boundary condition. In that way, the radiative flux is centered at the interfaces of the cells, as well as the divergence of radiative pressure. A von Neumann stability analysis of the modified scheme is done in section 2.A, which shows that this discretization for the source terms is unconditionally stable.

Asymptotic preserving scheme

Now that the choice for tσ s F u n`1 i is determined, we still have to specify our choice for α i`1 2 in eq. 2.1. α i`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. 1.24 as soon as the asymptotic regime is reached, i.e., large opacity and long timescale. Unfortunately, a standard HLL scheme does not have this property (see section 2.2.1). To tackle this issue and get an asymptotic preserving scheme, we choose

α i`1 2 " 1 1 ´3σ s i`1 2 ∆x ´1 ´f 2 i`1 2 ¯λì `1 2 λ í`1 2 c ´λì `1 2 ´λí `1 2 ¯, (2.14) with f i`1 2 " 1 2 `f n i `f n i`1 ˘.
The derivation of eq. 2.14 is done in section 2.B. Other choices can be done, which leads to other schemes, with other properties. For example, the choice done by [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF] leads to a solver more difficult to write, but that benefits from the known properties of an approximate Riemann solver. If σ i`1 2 ∆x goes to 0, α i`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 section 2.B. We show that

$ ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' % E n`1 i,0 " a r `T n`1 i,0 ˘4 F n`1 i,0 " 0 σ s i`1 2 F n`1 i`1 2 ,1 " ´c 3∆x `En`1 i`1,0 ´En`1 i,0 Ȇn`1 i,0 `ρc v T n`1 i,0 " E n i,0 `ρc v T n i,0 `c∆t 3∆x 2 ˜En`1 i`1,0 ´En`1 i,0 σ s i`1 2 ´En`1 i,0 ´En`1 i´1,0 σ s i´1 2 ¸, (2.15)
which is a standard discretization of eqs. 1.22 to 1.24. 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, depending on the test case, we can either enforce f ă 1 (section 2.2.6) or come back to a centered discretization of the source term (section 2.2.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 eq. 2.14 is only the first step to have an asymptotic preserving scheme in the case of a moving fluid (see section 2.3.2).

Coupling to hydrodynamics

Following González et al. [2007], the resolution of the whole system 1.25 describing radiation hydrodynamics is split into three steps:

1. update of the hydrodynamics quantities (eqs. 1.25a to 1.25c without the terms of energy and momentum exchange) using the well-balanced and all-regime solver developed in Padioleau et al.

[2019]; 2. update of the radiative quantities and gas temperature (eqs. 1.15 and 1.16) using the solver developed in section 2.1.1. During this step, the hydrodynamics quantities are frozen; 3. addition of source terms S Er puq and S Fr puq. For simplicity, all source terms which depend on the velocity are treated explicitly. The term σ s c F r in eqs. 1.25b and 1.25e is discretized using the well-balanced scheme proposed in section 2.1.1.2. All the other terms remain cell-centered. This splitting allows reducing the number of equations solved implicitly, making the method more efficient.

Algorithm for nonlinear implicit solver

Newton-Raphson method and linear solver

Because of the Eddington tensor, the eigenvalues in the numerical fluxes, and the a r T 4 g factor, the system is nonlinear. It is solved using a Newton-Raphson method. At each iteration, we have to solve a linear system. Because the system is large (p2 `dqN 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].

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 ´1Kx " b via solving AK ´1y " 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 ´1 is lower than A's one. Preconditioners used here are presented in section 5.1.

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

Numerical results

The implementation of the scheme developed in section 2.1 has been done in the code ARK-RT1 , a fork of the code ARK developed in Padioleau et al. [2019]. The hydrodynamics and gravity part of the solver is similar to ARK and is solved with a well-balanced and all-regime solver. More details about ARK-RT can be found in chapter 5.

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 1 model with a beam test and a shadow test, and the coupling to the hydrodynamics with radiative shocks. We also present a physical application about the stability of the ionization front in an H ii region in a massive pre-stellar dense core. To ease notations, we define the radiative temperature as T r " ´Er ar ¯1 4 .

Marshak wave

From Mihalas and 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, to test the asymptotic preserving scheme developed in section 2.1.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 0 " T r " 300 K, the initial radiative flux is F r " 0. We consider a perfect gas with γ " 5 3 . The hydrodynamics is frozen. The density is constant, such that ρc v " 1 J K ´1 cm ´3, the opacity is also x (cm) Gas temperature (K) This figure shows a snapshot of the gas temperature at time t f " 2 ˆ10 ´4 s, with and without the asymptotic correction and the reference solution. Spatial resolution is 400 points and the opacity is σ " 10 000 cm ´1. constant, with σ a " σ s " σ " 10 000 cm ´1, therefore σ∆x " 25. At time t " 0, a source is lit at the left boundary with T r " 1 000 K.

Reference solution With asymptotic correction Without asymptotic correction

The results are shown in fig. 2.1 at time t f " 2 ˆ10 ´4 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. 1.24. The relative L 2 error between the reference solution and the solution with α i`1 2 given by eq. 2.14 is 1.1%, whereas with the standard HLL scheme the relative L 2 error is 84%. Using the asymptotic correction, we recover the correct behavior in the asymptotic limit.

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 the 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 0 " T r " 300 K, the initial radiative flux is F r " 0. The opacity σ a " σ s " σ is now a function of space: σpxq " # 10 000 cm ´1 if x ă 0.5, 0 if x ą 0.5.

(2.16) At time t " 0, a source is lit at the left boundary with T r " 1 000 K. Figure 2.2 shows the radiative flux at time t f " 10 ´3 s. From eq. 2.4b, 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. 2.5, 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. 2.7 (blue curve), the spurious peak does not appear anymore and the constant steady state is reached.

Using the standard discretization of the source term eq. 2.5, one can show that the numerical scheme is unconditionally stable, in that the error between the numerical solution and the exact solution goes to 0 as ∆x and ∆t go to 0. The spurious peak seems to be due to a lack of precision in the integration of the source term. Using eq. 2.7, the source term is defined at the interfaces of the cells and balances the divergence of radiative pressure, also defined at the interfaces.

Beam

We now perform the same two-dimensional test as in González et al. [2007]; Richling et al. [2001]. The domain r´1, 1sˆr´1, 1s is discretized with 128ˆ128 cells. The initial temperature is T 0 " T r " 300 K, the initial radiative flux is F r " 0, the opacity is σ a " σ s " 0. At time t " 0, a beam with T g " T r " 1 000 K is introduced with an angle of 45 ˝. The beam is located at x " ´1 and y P r´0.875, ´0.75s. Because we are in the free-streaming regime, the propagation of the photons has to be followed, 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 2.3 shows the radiative energy at steady state. The eigenvalues in eq. 2.2 can either be fixed to ˘c or computed using eq. 2.3. Because there is no opacity, 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 2.4 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 stay 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.

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. 2.3 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.

Shadow

Let us now consider a two-dimensional test with source terms. Following González et al. [2007]; Hayes and Norman [2003], we consider a shadow test. The computational domain is a cylinder of length L " 1 cm and radius R " 0.12 cm. It is discretized with 280 ˆ80 cells. A spheroid clump is located at the center of the box, on the symmetric axis: pz c , r c q " p0.5, 0q. The extension of the clump is pz 0 , r 0 q " p0.1, 0.06q. Initially, the medium is at equilibrium with the radiation, with T 0 " T r " 290 K. We consider a homogeneous gas, with ρ 0 " 1 g cm ´3, except for the clump with density ρ 1 " 100ρ 0 . The boundary of the clump is smoothed: ρpr, zq " ρ 0 `ρ1 ´ρ0 1`exp ∆ with ∆ " 10 ˆ´z´zc

z 0 ¯2 `´r´rc r 0 ¯2 ´1˙.
The opacity in the medium is σ a " σ s " σ " σ 0 ´T T 0 ¯´3.5 ´ρ ρ 0 ¯2 with σ 0 " 0.1 cm ´1. At time t " 0, a source is lit at the left boundary with T r " 1 740 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 non-well-balanced scheme: the radiative flux source term is discretized using eq. 2.5. 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 λ ì`1 2 " maxp0.1 ˆc, λ max q and λ í`1 2 " minp´0.1 ˆc, λ min q, where λ max and λ min are given by eq. 2.3. M 1 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.

As in González et al. [2007]; Hayes and Norman [2003], we plot the radial profile of the radiative temperature at the right boundary (fig. 2.6). Using the P 1 model, the radiative pressure is isotropic, therefore the photons go around the obstacle immediately, heating the whole domain. Using the M 1 closure relation, the shadow is better preserved, the temperature behind the obstacle remains at its initial value, 290 K. Because we are not in the diffusion regime 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 section 2.2.2 and the well-balanced modification of the source term is not necessary.

Radiative shocks

Now that we have verified the properties of our scheme with the hydrodynamics frozen, we perform numerical tests to validate the coupling with hydrodynamics. We consider radiative shocks: the gas and the radiation exchange energy and momentum. Following González et al. [2007]; Hayes and Norman [2003]; Ensman [1994], we consider a one-dimensional homogeneous medium, with ρ " 7.78 ˆ10 ´10 g cm ´3 and σ a " σ s " σ " 3.1 ˆ10 ´10 cm ´1. We consider a perfect gas with an adiabatic coefficient γ " 7 5 and a mean molecular weight µ " 1. The length of the domain is 7 ˆ10 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, the initial velocity of the fluid is set to u 0 . According to the value of u 0 , the shock will be subcritical or supercritical. See González et al. [2007] for more details. To compare our results with González et al. [2007]; Hayes and Norman [2003]; Ensman [1994], we plot the temperature as a function of x i " x ´u0 t.

Subcritical shock

We first consider a subcritical shock, the initial velocity is set to u 0 " ´6 km s ´1. Figure 2.7 shows the gas temperature, the radiative temperature, and the reduced flux at three different times: 1.7 ˆ10 4 s, 2.8 ˆ10 4 s, and 3.8 ˆ10 4 s. As expected, the gas and the radiation are not at equilibrium, before nor after the shock. The gas temperature reaches 1 135 K, as in González et al. [2007], whereas it is only 850 K in Ensman [1994].

Supercritical shock

We consider now a supercritical shock, where the initial velocity is set to u 0 " ´20 km s ´1. Figure 2.8 shows the gas temperature, the radiative temperature, and the reduced flux at three different times: 4 ˆ10 3 s, 7.5 ˆ10 3 s, and 1.3 ˆ10 4 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 the subcritical shock's radiative precursor, as intended, and the temperature reaches 5 000 K, as in Ensman [1994]. We also recover the Zel'dovich spike.

We recover the expected results with both subcritical and supercritical shocks, therefore our code is well-suited to study radiation hydrodynamics problems.

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.

Model

We consider the early stage of the development of an H ii region in a massive prestellar dense core [START_REF] Churchwell | Ultra-compact hii regions and massive star formation[END_REF]. 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, i.e., the early stage of the development of the massive star. The temperature reached at this location by infrared heating is of the order of 1 000 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). Highenergy photons emitted by the YSO when entering the main sequence will start to ionize the surrounding gas. This will trigger 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. 1.25. The only photons able to ionize the gas are emitted by the YSO, i.e., there is no local source of ionizing photons. Following [START_REF] Tremblin | Ionization impact on molecular clouds and star formation : numerical simulations and observations[END_REF], we need to modify this model to take into account photo-chemistry and thermal balance. We define the fraction of ionization X " n H `{n H where n H " n H ``n H 0 , n H `is the number of ionized atoms and n H 0 is the number of cold atoms. The evolution of the number of ionized atoms is just the number of incoming photons that interact with the gas minus the number of ionized atoms that recombine (on the spot approximation, see [START_REF] Lesaffre | Dynamical aspects of the interstellar medium[END_REF]. Therefore,

B t pρXq `∇ ¨pρXuq " σ γ F γ n H p1 ´Xq ´mH βX 2 n 2 H , (2.17)
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 ´10 `T 1 K ˘´0.75 cm 3 s ´1 with T the temperature of thermodynamic equilibrium [START_REF] Black | The physical state of primordial intergalactic clouds[END_REF].

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. It is the only source of heating during the ionization, hence the heating rate is given by p1 ´Xqn H F γ σ γ 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 [START_REF] Lesaffre | Dynamical aspects of the interstellar medium[END_REF] to recover the observed temperature around 1 000 K. From [START_REF] Tremblin | Ionization impact on molecular clouds and star formation : numerical simulations and observations[END_REF], the cooling rate is given by βX 2 n 2 H k b T {pγ ´1q. We also add a term of Newtonian forcing B t T g " Tg´T forcing τ 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 H m H and σ a " σ s " σ " σ γ n H , we finally have to solve the following system: %

B t ρ `∇ ¨pρuq " 0 B t pρuq `∇ ¨pρu b u `pIq " ρg `σp1 ´Xq c F r B t pρEq `∇ ¨ppρE `pquq " ρg ¨u `cσp1 ´XqE r ´β ρ 2 X 2 m 2 H k b T g γ ´1 ´ρc v T g ´Tforcing τ forcing B t E r `∇ ¨Fr " ´cσp1 ´XqE r B t F r `∇ ¨Pr " ´cσp1 ´XqF r B t pρXq `∇ ¨pρXuq " σp1 ´XqcE r m H e γ ´βρ 2 X 2 m H .
(2.18)

In this test, we will use the M 1 solver with the asymptotic correction presented in section 2.1.1.3, but we do not use the well-balanced discretization of the source term because of stability issues that will be discussed in section 2.3.1.

Setup

We consider a square domain with a side 5 AU long. We use a setup similar to Padioleau et al. [2019] for compressible convection simulations. The temperature is set at the top and the bottom of the box at 500 K and 1 000 K, 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 7 s. These parameters are chosen to trigger the initial convective motions. We also set the pressure at the bottom: 10 -3 dyn cm -2 [START_REF] Herpin | S-bearing molecules in massive dense cores ***[END_REF]. The density and the pressure are linked by the ideal gas law: p " ρk b Tg m H µ , where µ is the mean molecular weight. The non-ionized medium is made of hydrogen, with µ 1 " 1. When the medium is fully ionized, it is made of atomic nucleus and electrons, so twice as many particles for the same mass. Because the distribution of nucleus and electrons is homogeneous, the mean molecular weight is µ 2 " 0.5. When the medium is partially ionized, we take µ as the mean of the previous values balanced by the fraction of ionization, i.e., µ " p1 ´Xqµ 1 `Xµ 2 . The density is initialized with the recursive formula p i`1 ´pi " 1 2 pρ i `ρi`1 q g∆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

$ ' ' & ' ' % upx, yq " 2 ¨10 ´4c s sin ˆ2π x ´xmid L x ˙sin ˆπ y ´ymid L y vpx, yq " 2 ¨10 ´4c s cos ˆ2π x ´xmid L x ˙cos ˆπ y ´ymid L y ˙, (2.19)
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 σ " σγ ρ m H with σ γ " 6 ˆ10 ´18 cm 2 [START_REF] Lesaffre | Dynamical aspects of the interstellar medium[END_REF]. The radiative energy and flux are set to 0 and the medium is not ionized (X " 0). We initialize the hydrodynamics variables with the steady state described previously. At time t " 0, the bottom boundary of the region is ionized: the reduced flux is set to 1 and the radiative energy is set to F˚eγ c with F ˚" 3 ˆ10 17 cm ´2 s ´1 in the boundary. The boundary conditions for the hydrodynamics variables remain unchanged.

Results

As the initial condition is such that E r is close to 0, 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. 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. The fraction of ionization, which is always between 0 and 1, reaches values between 10 -12 and 10 -6 . The effect of the preconditioner and the MPI domain decomposition is discussed in section 5.2.3. However, the numerical noise does not affect the position of the ionization front.

The stability of the ionization front is an issue that has been discussed for a long time in the literature [START_REF] Mizuta | Hydrodynamic instability of ionization fronts in HiiRegions[END_REF]. 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 [START_REF] Bisbas | starbench: the D-type expansion of an Hii region[END_REF], fig. A3). The dependence of the instability on the grid cast doubts about a possible physical regime. Our test case shows that even with convective motions of large amplitude, the ionization front remains stable.

Discussion and conclusion

Well-balanced discretization of the source term

In section 2.1.1.2 we have proposed 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 and Despres [2008].

Asymptotic limit for radiation hydrodynamics

In section 2.1.1.3, we have presented 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 section 1.C. Most of the schemes proposed in the literature do not preserve this asymptotic regime (e.g., Berthon and Turpault 2011; González et al. 2007). 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 and Despres [2008]. This Lagrange-projection-like strategy is explored in chapter 3 to treat both the asymptotic preserving and well-balanced issues.

2.A Von Neumann stability analysis for the wellbalanced modification of the source term

For simplicity, we split eq. 1.15b into a pure hyperbolic problem B t F r `c2 ∇ ¨Pr " 0 and a source problem B t F r " ´cσ s F r . We focus on the one-dimensional source problem, with periodic boundary conditions on the domain r0, T s ˆr0, 1s 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 (eq. 2.7), the source problem is discretized as

F n`1 j `rj´1 2 F n`1 j´1 `´r j´1 2 `rj`1 2 ¯F n`1 j `rj`1 2 F n`1 j`1 " F n j , (2.20)
with r j`1 2 "

cσ s j`1 2 ∆t 4
. We define the function F n , piecewise constant, such that

F n pxq " # F n j if x j´1 2 ă x ă x j`1 2 0 otherwise. (2.21)
This function is then extended to R by periodicity. F n can now be expanded in a Fourier series:

F n pxq " ÿ kPZ F n pkqe 2ikπx , (2.22)
with F n pkq "

ż 1 0 F n pxqe ´2ikπx dx. (2.23)
We can define the 2-norm of the function F n :

||F n || 2 " ˆż 1 0 pF n pxqq 2 dx ˙1 2 " ˜ÿ kPZ F n pkq 2 ¸1 2 " ˜J ÿ j"1 ∆x `F n j ˘2¸1 2 .
(2.24)

We apply the Fourier transform to eq. 2.20:

F n`1 pkq ´1 `rj´1 2 e ´2ikπ∆x `rj´1 2 `rj`1 2 `rj`1 2 e 2ikπ∆x ¯" F n pkq. (2.25)
We define the amplification factor Apkq as

Apkq " 1 1 `rj´1 2 e ´2ikπ∆x `rj´1 2 `rj`1 2 `rj`1 2 e 2ikπ∆x , ( 2
1 |Apkq| 2 " ´1 `rj´1 2 `rj`1 2 ¯2 `r2 j´1 2 `r2 j`1 2 `2r j´1 2 r j`1 2 cosp4kπ∆xq `2r j´1 2 ´1 `rj´1 2 `rj`1 2 ¯cosp2kπ∆xq `2r j`1 2 ´1 `rj´1 2 `rj`1 2 ¯cosp2kπ∆xq " ¨1 `rj´1 2 p1 ´cosp2kπ∆xqq looooooooooooomooooooooooooon ě0 `rj`1 2 p1 ´cosp2kπ∆xqq looooooooooooomooooooooooooon ě0 ‹ ' 2 `´r j´1 2 sinp2kπ∆xq ´rj`1 2 sinp2kπ∆xq ¯2 looooooooooooooooooooooooomooooooooooooooooooooooooon ě0 ě 1.
(2.27) Since 1 |Apkq| 2 ě 1, we have |Apkq| ď 1.

2.B Numerical scheme in the diffusive limit

We consider the numerical scheme developed in section 2.1.1.3 in the asymptotic regime, with σ s i`1 2 ∆x Ñ 8. Following section 1.4, we introduce the rescaling parameter ε to write the time (resp. the opacity) as Ă ∆t " ε∆t (resp. r σ a " εσ a and r σ s " εσ s ). Using the P 1 closure relation, we have λ ì`1 2 " ´λí

`1 2 " c ? 3 and
ε 2 E n`1 i " ε 2 E n i ´Ă ∆t ∆x ´εα i`1 2 F i`1 2 ´εα i´1 2 F i´1 2 ¯`c Ă ∆t r σ a i ´ar `T n`1 i ˘4 ´En`1 i (2.28a) ε 2 F n`1 i " ε 2 F n i ´Ă ∆t ∆x ´εP i`1 2 ´εP i´1 2 ¯´c Ă ∆t 2 ´r σ s i`1 2 F n`1 i`1 2 `r σ s i´1 2 F n`1 i´1 2 ¯(2.28b) ε 2 ρc v T n`1 i " ε 2 ρc v T n i ´cĂ ∆t r σ a i ´ar `T n`1 i ˘4 ´En`1 i ¯. (2.28c)
Radiative variables are expanded, e.g., E n i " E n i,0 `εE n i,1 `O pε 2 q. Expanding eqs. 2.28a and 2.28b at order 0 leads to

# E n`1 i,0 " a r `T n`1 i,0 ˘4 F n`1 i,0 " 0.
(2.29)

At first order for eq. 2.28b, we have

r σ s i`1 2 F n`1 i`1 2 ,1 `r σ s i´1 2 F n`1 i´1 2 ,1 " ´c 3 E n`1 i`1,0 ´En`1 i´1,0 ∆x `F n`1 i`1,0 ´2F n`1 i,0 `F n`1 i´1,0 ? 3∆x loooooooooooooomoooooooooooooon "0 . (2.30)
Using the boundary condition given by eq. 2.13, we have

r σ s i`1 2 F n`1 i`1 2 " ´c 3∆x `En`1 i`1,0 ´En`1 i,0 ˘(2.31)
in the whole domain. Now, we consider the sum of eqs. 2.28a and 2.28c expanded at second order. If

α i`1 2
" 1, we have

ε 2 `En`1 i,0 `ρc v T n`1 i,0 ˘" ε 2 `En i,0 `ρc v T n i,0 ε Ă ∆t ∆x c 2 ? 3 ´αi`1 2 `En`1 i`1,0 ´En`1 i,0 ˘´α i´1 2 `En`1 i,0 ´En`1 i´1,0 ˘ε 2 Ă ∆t 2∆x α i`1 2 ˆF n`1 i`1,1 `F n`1 i,1 ´c ? 3 `En`1 i`1,1 ´En`1 i,1 ˘ε 2 Ă ∆t 2∆x α i´1 2 ˆF n`1 i,1 `F n`1 i´1,1 ´c ? 3 `En`1 i,1
´En`1 i´1,1 ˘˙, (2.32) whereas the asymptotic development of a standard discretization of eq. 1.24 would be

ε 2 `En`1 i,0 `ρc v T n`1 i,0 ˘" ε 2 `En i,0 `ρc v T n i,0 ˘`ε c 3 ∆t ∆x 2 ˜En`1 i`1,0 ´En`1 i,0 σ s i`1 2 ´En`1 i,0 ´En`1 i´1,0 σ s i´1 2 ¸.
(2.33) So, we are looking for α i`1 2 such that the term of order 1 in eq. 2.32 becomes a term of order 2 with the expected coefficient of diffusion c 3σ s i`1 2 and the term of order 2 becomes a term of order 3 and therefore negligible. In other words, we want the asymptotic development of α i`1 2 to be

2ε ? 3σ s i`1 2
∆x . One way to achieve this is to take

α i`1 2 " 1 1 `?3σ s i`1 2 ∆x 2
.

(2.34) However, in the general case, we do not have λ ì`1 2 " ´λí

`1 2 " c ? 3 .
We can then replace eq. 2.34 by

α i`1 2 " 1 1 ´3σ s i`1 2 ∆x λ ì`1 2 λ í`1 2 c ´λì `1 2 ´λí `1 2 ¯.
(2.35)

Unfortunately, in numerical tests with σ s ∆x close to 1, the condition f ă 1 is not preserved. Because f is close to 1 in this case, we write

α i`1 2 " 1 1 ´3σ s i`1 2 ∆x ´1 ´f 2 i`1 2 ¯λì `1 2 λ í`1 2 c ´λì `1 2 ´λí `1 2 ¯.
(2.36)

We use

f i`1 2 " 1 2 `f n i `f n i`1
˘because numerical experiments have shown good results using this form. In the diffusion regime, because F n`1 i,0 " 0, we recover eq. 2.35. Now that we have the form of α i`1 2

, we can check that the proposed scheme is asymptotic preserving. We have

α i`1 2 " 2ε ? 3 r σ s i`1 2 ∆x `O `ε2 ˘. (2.37)
Therefore,

α i`1 2 F i`1 2 " ´ε c 3 r σ s i`1 2 E n`1 i`1,0 ´En`1 i,0 ∆x `O `ε2 ˘. (2.38)
We finally have

E n`1 i,0 `ρc v T n`1 i,0 " E n i,0 `ρc v T n i,0 `cĂ ∆t 3∆x 2 ˜En`1 i`1,0 ´En`1 i,0 r σ s i`1 2 ´En`1 i,0 ´En`1 i´1,0 r σ s i´1 2 ¸. (2.39)
Equations 2.29, 2.31 and 2.39 are standard discretization of eqs. 1.22 to 1.24, so this scheme is asymptotic preserving.

Chapter 3

An all-regime-like asymptotic preserving scheme The numerical scheme we have shown in the previous chapter does not preserve the admissible states E r ą 0 and f ď 1. We present in this chapter the development of another numerical scheme that preserves these admissible states, as well as the asymptotic preserving property. This scheme being inspired by all-regime schemes for compressible hydrodynamics, we also call it "all-regime scheme". We derive it using a time-explicit integration. A possible extension for a time-implicit solver is presented in chapter 4. For the sake of simplicity, we only consider the case where the hydrodynamics is frozen. The coupling to hydrodynamics could be done as in section 2.1.2.

This chapter is organized as follows. We present our numerical scheme in section 3.1. Buet and Despres [2008] have derived a Lagrange-projection scheme for the M 1 model. Inspired by Chalons et al. [2016], we adapt this scheme with a "diffusion-transport" splitting. Because this leads to a scheme of stencil 2, we also present an all-regime scheme of stencil 1. In section 3.2, we present some numerical test cases to show the good properties of our solvers. Finally, we reach our conclusion in section 3.3 and we discuss the limits of the schemes.

Numerical scheme

We use the same notations as in chapter 2: we note ∆x the step along the xdirection. ∆t is the time interval between the current time t n and t n`1 . We write x i the center of the cell i and x i`1 2 the interface between the cell i and the cell i `1. We use the notation u n i to represent the averaged quantity associated with the field u at time t n in the cell i (finite volume). Finally, we note u n i`1 2 to represent the quantity associated with the field u at time t n and at the interface between cells i and i `1.

Because we use a time-explicit scheme, the time step ∆t is restricted by a CFL condition. For the sake of simplicity, we always use ∆t ă ∆x c in the numerical tests. For physical reasons, the fastest phenomena occur at the speed of light, therefore we use this velocity to restrict the time step, even though slower phenomena could allow us to use a larger time step.

Rewriting the M 1 model

In this chapter, we adapt a scheme developed for hydrodynamics to the M 1 model. We first rephrase our model to have the same form as the Euler equation. Following Buet and Despres [2008], we introduce a variable P R that obeys the transport equation B t `∇ ¨p βq " 0,

with β already defined in section 1.B.2, given by

β " 1 f 3χ ´1 2 F r E r . (3.2)
Despite the resemblance with the mass conservation equation of Euler equations, is not the density of the fluid ρ. However, it will play a similar role in the derivation of the scheme.

In chapter 1, we have shown that

$ & % F r " pq `Er q β P r " β b F r c 2 `qI, (3.3) ' ' ' ' ' % r τ i " τ n i `∆t ∆m i ´βi `1 2 ´βi ´1 2 r v i " v n i ´c2 ∆t ∆m i ´qi `1 2 ´qi ´1 2 r e i " e n i ´∆t ∆m i ´qi `1 2 β i`1 2 ´qi ´1 2 β i´1 2 ¯.
(3.11)

Numerical fluxes β ˚and q ˚will be specified later.

The update of the conservative variables p , v, eq can be written as

$ ' ' ' ' & ' ' ' ' % r L i r i " n i r L i Ą p vq i " p vq n i ´c2 ∆t ∆x ´qi `1 2 ´qi ´1 2 r L i Ą p eq i " p eq n i ´∆t ∆x ´qi `1 2 β i`1 2 ´qi ´1 2 β i´1 2 ¯, (3.12)
where r

L i " 1 `∆t ∆x ´βi `1 2 ´βi ´1 2 ¯.
The last point needed for the resolution of the diffusion step to be complete is the choice of the numerical fluxes β i`1 2 and q i`1 2 . We have chosen to use the fluxes given by Buet and Despres [2008]:

$ ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' % β i`1 2 " n i e n i 3c 2 ``β n i ˘2 β n i ` n i`1 e n i`1 3c 2 `´β n i`1 ¯2 β n i`1 n i e n i 3c 2 ``β n i ˘2 ` n i`1 e n i`1 3c 2 `´β n i`1 ¯2 `?3 4c
q n i ´qn i`1 n i e n i 3c 2 ``β n i ˘2 ` n i`1 e n i`1 3c 2 `´β n i`1 ¯2 q i`1 2 " 3c 2 ``β n i ˘2 n i e n i q n i `3c 2 `´β n i`1 ¯2 n i`1 e n i`1 q n i`1 3c 2 ``β n i ˘2 n i e n i `3c 2 `´β n i`1 ¯2 n i`1 e n i`1 `4c ? 3
β n i ´βn i`1 3c 2 ``β n i ˘2 n i e n i `3c 2 `´β n i`1 ¯2 n i`1 e n i`1 . (3.13)
The derivation of eq. 3.13 is done in section 3.C. Buet and Despres [2008] show that these numerical fluxes preserve the admissible states, i.e., r e i ą 0 and |r v i | cr e i ď 1 (see section 3.D).

Other choices for β ˚and q ˚can be done. Riemann invariants for the diffusion step (eq. 3.10) could be used. It has been explored with S.Bulteau:

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % β i`1 2 " pq n i q ? 3 4 b pc `βn i q `c `βn i`1 ˘´`q n i`1 ˘?3 4 b pc ´βn i q `c ´βn i`1 pq n i q ? 3 4 b pc `βn i q `c `βn i`1 ˘``q n i`1 ˘?3 4 b pc ´βn i q `c ´βn i`1 qi `1 2 " a q n i q n i`1 ˜`c ´βn i`1 ˘pc `βn i q `c `βn i`1 ˘pc ´βn i q ¸1 ? 3 . (3.14)
First numerical experiments have shown results similar to those obtained with eq. 3.13, but with a higher computational cost. Finally, Buet and Despres [2008] show that the solution ´r E i , r F i ¯is independent of the initial value of the density n i .

Transport step

The scheme proposed by Buet and Despres [2008] uses a Lagrange-remap solver. As our all-regime solver, it is split into two steps. The first one, the Lagrangian step, has the same form as our diffusion step. The fields evolve using the Lagrangian formalism. The mesh is distorted to follow the matter. This new mesh is then projected onto the original mesh during the remap step. This approach, although robust, might be difficult to implement in the multi-dimensional case. Inspired by Chalons et al. [2016],

we propose here to replace the remap step by the resolution of the transport step eq. 3.7. The mesh is not changed during the resolution of the diffusion step, therefore solving the transport step is easier than computing the remap step of the Lagrange-remap solver.

Following Padioleau et al.

[2019], the transport system (eq. 3.7) can be written

B t ϑ `Bx pϑβq ´ϑB x β " 0, (3.15)
for ϑ P t , v, eu. We discretize it as

ϑ n`1 i " r ϑ i ´∆t ∆x ´r ϑ i`1 2 β i`1 2 ´r ϑ i´1 2 β i´1 2 ¯`r ϑ i ∆t ∆x ´βi `1 2 ´βi ´1 2 ¯. (3.16)
The term defined at the interface r ϑ i`1 2 is defined by the upwind choice with respect to the velocity β i`1 2 :

r ϑ i`1 2 " $ & % r ϑ i if β i`1 2 ě 0 r ϑ i`1 if β i`1 2 ď 0.
(3.17)

Using eq. 3.12, the update of the conservative quantities is then

w n`1 i " w n i ´∆t ∆x pF ˚p r w i , r w i`1 q ´F˚p r w i´1 , r w i qq , (3.18)
where w " p , F, Eq T and

F ˚¨¨ L F L E L ' , ¨ R F R E R '' " $ ' ' ' ' ' ' & ' ' ' ' ' ' % ¨ L β FL β ˚`c 2 q EL β ˚`q ˚β˚' if β ˚ě 0, ¨ R β FR β ˚`c 2 q ER β ˚`q ˚β˚' if β ˚ď 0, (3.19)
This choice of discretization allows the total update to be conservative. Furthermore, because the transport step is simply the transport of the quantities updated by the diffusion step, all properties preserved by the diffusion step are also preserved by the transport step. In particular, the transport step preserves the admissible states e n`1 i ą 0 and v n`1 i ce n`1 i ď 1, as well as the asymptotic preserving property. For the same reason, the solution at time t n`1 is also independent of the value of the density at time t n . Therefore, in numerical tests, we set to 1 in each cell at the beginning of each time step.

Scheme of stencil 1

When looking at the total update using this diffusion-transport splitting, from time t n to t n`1 , the resulting scheme is of stencil 2. Looking at eq. 3.18, the update of ϑ n`1 i involves, among others, r ϑ i˘1 , which is computed thanks to ϑ n i , ϑ n i˘1 and ϑ n i˘2 . To tackle this issue, eq. 3.17 can be replaced by

r ϑ i`1 2 " # ϑ n i if β i`1 2 ě 0 ϑ n i`1 if β i`1 2 ď 0.
(3.20) Equation 3.18 becomes

w n`1 i " w n i ´∆t ∆x `F˚`wn i , w n i`1 ˘´F ˚`w n i´1 , w n i ˘˘, (3.21)
which is indeed a scheme of stencil 1. This scheme can also be derived as an approximate Riemann solver. It is done in section 3.E. We show in section 3.D that the admissible states E r ą 0 and |Fr| cEr ď 1 are preserved. Looking at eqs. 3.19 and 3.21, the update of E n`1 i and F n`1 i does not depend on . In the numerical tests, we set it to 1 in the whole domain, at the beginning of each time step.

Source terms

Now that we have presented two schemes for the hyperbolic part of the M 1 model, let us consider the addition of the source terms. We aim at designing an asymptotic preserving and well-balanced scheme, i.e., it should capture the diffusive limit eq. 1.24 and the correct steady state, as explained in section 2.1.1.2.

As in section 2.1.1, we use a splitting strategy to treat the source terms. Once we have obtained an intermediate solution E n`1 i , F n`1 i with the scheme described in sections 3.1.3 and 3.1.4 or with the one presented in section 3.1.5, we still have to solve

$ ' & ' % B t E " cσ a `ar T 4 ´EB t F " ´cσ s F B t pρc v T q " ´cσ a `ar T 4 ´E˘, (3.22)
with initial condition E n`1 i , F n`1 i , T n i . To avoid restricting the time step, we use an implicit integration to solve this system, as done in section 2.1.1. This leads to a semi-implicit scheme.

As Buet and Despres [2008], we focus on the radiative flux source term. In section 2.1.1.2, we use an average of a face discretization of this source term.

σ s F is discretized in space as 1 2 ´σs i`1 2 F i`1 2 `σs i´1 2 F i´1 2 ¯, with $ ' & ' % σ s i`1 2 " 1 2 `σs i `σs i`1 F n`1 i`1 2 " 1 2 `F n`1 i `F n`1 i`1 ˘. (3.23)
This leads to solve a linear system and might not preserve the admissible states. Because β and F are colinear, Buet and Despres [2008] suggest using a discretization proportional to

β i`1 2 `βi ´1 2
. However, this requires the resolution of a nonlinear system. An algorithm to solve the nonlinear system arising from the discretization of the M 1 model while preserving the admissible states is presented in chapter 4. Because of the issues raised here, we will use a cell-centered discretization for the radiative flux source term, while being aware of its limitations.

To write eq. 3.13, we first wrote in section 3.C some jump relations:

$ ' ' ' & ' ' ' % ´qi `1 2 ´qn i ¯`4c ? 3
E n i 3c 2 `pβ n i q 2 ´βi `1 2 ´βn i ¯" 0 ´qi `1 2 ´qn i`1 ¯´4c ? 3
E n i`1 3c 2 ``β n i`1 ˘2 ´βi `1 2 ´βn i`1 ¯" 0. (3.24) n i`1 e n i`1 3c 2 `´β n i`1 ¯2 β n i`1 n i e n i 3c 2 ``β n i ˘2 ` n i`1 e n i`1 3c 2 `´β n i`1 ¯2 `?3 4c
ˆqn i ´σs i`1 2 ∆x 2c n i v n i ˙´ˆq n i`1 `σs i`1 2 ∆x 2c n i`1 v n i`1 ˙ n i e n i 3c 2 ``β n i ˘2 ` n i`1 e n i`1 3c 2 `´β n i`1 ¯2 q i`1 2 " 3c 2 ``β n i ˘2 n i e n i ˆqn i ´σs i`1 2 ∆x 2c n i v n i ˙`3c 2 `´β n i`1 ¯2 n i`1 e n i`1 ˆqn i`1 `σs i`1 2 ∆x 2c n i`1 v n i`1 3c 2 ``β n i ˘2 n i e n i `3c 2 `´β n i`1 ¯2 n i`1 e n i`1 `4c ? 3
β n i ´βn i`1 3c 2 ``β n i ˘2 n i e n i `3c 2 `´β n i`1 ¯2 n i`1 e n i`1
.

(3.26) We show in section 3.F that, using eq. 3.26 instead of eq. 3.13, the two schemes developed previously are asymptotic preserving. We also verify this property numerically in section 3.2.2.

Numerical results

We perform a series of verification tests to validate different properties of the two schemes developed in section 3.1. We compare them with a standard HLL solver. Furthermore, we refer to the scheme presented in sections 3.1.3 and 3.1.4 as the "split scheme" and to the solver presented in section 3.1.5 as the "unsplit scheme".

We first verify that both schemes are conservative, especially the split scheme, with a two-dimensional Riemann problem. We then look at numerical diffusion with a beam test. We also show some performance results. Afterward, we explore the asymptotic preserving property with a Marshak wave. Next, we study the well-balanced behavior by reaching a steady state with a jump of opacity. Finally, we highlight all these properties with a shadow test.

Numerical results without source terms

2D Riemann problem

Let us first consider the test described in Blachère and Turpault [2016]. It is a 2D Riemann problem with four states. The domain r0, 1s ˆr0, 1s is discretized with 256 ˆ256 cells. There is no opacity: σ a " σ s " 0. The initial temperature is T 0 " T r " 1 000 K and the radiative flux is set to p1 ´10 ´8q cE r F , where F is a unit 0 0.5 1 0 0.5 1 2), the split scheme and the unsplit scheme.

F " ˆ1 0 ˙Ñ F " ˆ0 ´1˙Ó F " ˆ0 1 ˙Ò F " ˆ´1 0 ˙Ð Figure 3.
vector, piecewise constant. The domain is cut into four states. In each region, the direction of the radiative flux is constant, see fig. Figure 3.3 shows the evolution of the relative error between the expected total radiative energy and the one actually computed in the box at each step. For this test, we use periodic boundary conditions, so no energy should enter or leave the box, therefore, the error should remain at machine precision. Here, it oscillates around the value 10 ´12 . This shows that even with the split scheme, the total energy is conserved. Indeed, the choice of discretization for the transport step allows us to have a conservative scheme. The split scheme being conservative relies on, among others, the fact that we use the same β ˚for the discretization of both diffusion and transport steps.

Beam

We perform the same test as in section 2.2.3, but we do not discuss the properties of the HLL solver, we use it here as a reference. It is the same test as in González et al. [2007]; Richling et al. [2001]. It is the propagation of a beam in the free-streaming regime, with σ a " σ s " 0. The domain r´1, 1s ˆr´1, 1s is discretized with 128 ˆ128 cells. The initial temperature is T 0 " T r " 300 K and the initial radiative flux is F r " 0. At time t " 0, a beam with T g " T r " 1 000 K is introduced with an angle of 45 ˝at x " ´1 and y P r´0.875, ´0.75s.

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. We compare both split and unsplit solvers with a standard HLL solver. The eigenvalues used in the HLL solver can be fixed to ˘c or computed (see eq. 2.2). The impact of this choice has already been discussed in section 2.2.3.

Figure 3.4 shows a horizontal cut at the middle height once the steady state is reached. Using the HLL solver, the full width at half maximum is 22 (resp. 29) cells with computed (resp. fixed) eigenvalues. With both split and unsplit schemes, it is 26 cells.

Table 3.1 shows the computational time needed to reach the steady state with the different solvers: HLL with fixed and computed eigenvalues, and both split and unsplit schemes. Using fixed eigenvalues for the HLL solver allows us to save 40% of the computational time compared to using computed eigenvalues. To reduce computational time with computed eigenvalues, González et al. [2007] suggest tabulating them. They are computed once at the beginning of the simulation and the value needed is then interpolated. The split scheme requires solving two steps, the diffusion step followed by the transport step, hence the increase of 80% of computational cost compared to the HLL solver with computed eigenvalues. The unsplit scheme, which does not require solving these two steps, allows a decrease of 20% of the computational time compared to the HLL solver with computed eigenvalues and an increase of 40% compared to the HLL solver with fixed eigenvalues.

Overall, the unsplit scheme seems to be a good compromise between numerical diffusion and computational cost.

Numerical results with source terms

Marshak wave

Now that we have confirmed the good behavior of our schemes when there is no opacity, let us consider a one-dimensional Marshak wave, in the diffusive limit, to test their asymptotic behavior. We perform a test similar to section 2.2.1. Because we use time-explicit integration here, we cannot reach easily the same final time, therefore, we change some numerical values.

The interval r0, 1s is discretized with 100 points. Initially, the medium is at equilibrium with the radiation: T 0 " T r " 5.87 ˆ10 8 K, the radiative flux is F r " 0, the opacity is constant, with σ a " σ s " 9 900 m ´1, such that σ s ∆x " 100 and ρc v " 1.53 J K ´1 m ´3. At time t " 0, a source is lit at the left boundary with T r " 5.87 ˆ10 9 K.

The results are shown in fig. 3.5 at time t f " 3 ˆ10 ´10 s. We compare different solutions: a reference solution, given by a standard discretization of the diffusion equation 1.24, a standard HLL scheme, the split and unsplit schemes. We also discuss the impact of the source terms in the numerical fluxes β ˚and q

˚.
As shown in chapter 2, the solution obtained with the HLL solver is dominated by the numerical diffusion and does not capture the asymptotic behavior, the relative L 2 error between this solution and the reference solution is around 192%. Using eq. 3.13, without source terms in the numerical fluxes, leads to a relative L 2 error of 194% between the solutions obtained with both split and unsplit solvers and the reference x (m) Gas temperature (K) Reference solution HLL Split, β ˚and q ˚with source terms (eq. 3.26) Unsplit, β ˚and q ˚with source terms (eq. 3.26) Split, β ˚and q ˚without source terms (eq. 3.13) Unsplit, β ˚and q ˚without source terms (eq. 3.13) Figure 3.5 -Marshak wave simulation. This figure shows a snapshot of the gas temperature at time t f " 3 ˆ10 ´10 s, with a reference solution and different solvers: a standard HLL scheme, the split and unsplit schemes, with and without the modification to take into account the source terms in β ˚and q ˚. solution. These solutions are also dominated by numerical diffusion. To tackle this issue, pressure terms in the numerical fluxes are modified to take into account the source terms, as shown by eq. 3.26. With this choice, the relative L 2 error drops to 12.8% and 14.4% for the split and unsplit solvers, respectively. We recover numerically the result shown in section 3.F.

Steady state with a jump of opacity

In the previous cases, the opacity was constant; let us now consider a test with a jump of opacity, in the one-dimensional case, as done in section 2.2.2. We use this test to discuss the well-balanced behavior of the schemes.

The length of the computational domain is 1 cm, it is discretized with 100 cells. Initially, the medium is at equilibrium with the radiation: T 0 " T r " 300 K, the initial radiative flux is F r " 0 and ρc v " 1 J K ´1 cm ´3. The opacity is now a function of space: σ a pxq " σ s pxq " σ pxq, with σ pxq "

# 50 cm ´1 if x ă 0.5, 0 if x ą 0.5. (3.27)
At time t " 0, a source is lit at the left boundary, with T r " 5 000 K. Figure 3.6 shows the radiative flux at time t f " 10 ´6 s, once the steady state is reached. From section 2.1.1.2, when the steady state is reached, we expect the radiative flux to be constant in the box. Using a standard HLL scheme, the value reached in the second half of the box is 50% more than the expected value. As discussed in section 2.2.2, this is due to a lack of precision in the numerical integration of the source terms and can be tackled by using higher order integration of the source terms, which leads to an interface-centered discretization. x (m) Radiative flux (normalized units) HLL Split, β ˚and q ˚with source terms (eq. 3.26) Unsplit, β ˚and q ˚with source terms (eq. 3.26) Split, β ˚and q ˚without source terms (eq. 3.13) Unsplit, β ˚and q ˚without source terms (eq. 3.13) Figure 3.6 -Simulation of a steady state with a jump of opacity. The opacity is piecewise constant, a jump is located at x " 0.5 m (gray line). This figure shows a snapshot of the radiative flux at time t f " 10 ´6 s, with a standard HLL solver, and both the split and unsplit solvers.

Using both split and unsplit solvers, when the numerical fluxes do not take into account the source terms (eq. 3.13), the radiative flux in the second half of the box is also 50% more than the expected value. However, using eq. 3.26, the error is much smaller. We only observe the same spurious peak located at the discontinuity of opacity as in section 2.2.2. Because we still use a splitting strategy to treat implicitly the source terms, this can be improved, either by using β ˚to compute the source terms, as suggested by Buet and Despres [2008], or by using a discretization at the interfaces of the cells, as done in section 2.1.1.2. Both choices would require solving a (non)linear system while preserving the admissible states. An algorithm able to do that is explored in chapter 4.

Shadow test

Finally, let us consider the same two-dimensional test as in section 2.2.4. Again, we do not discuss the properties of the HLL solver, nor the M 1 model, but we compare the results obtained with our all-regime scheme to those obtained with the HLL solver. Following González et al. [2007]; Hayes and Norman [2003], we consider a shadow test. The computational domain is a cylinder of length L " 1 cm and radius R " 0.12 cm. It is discretized with 280 ˆ80 cells. A spheroid clump is located at the center of the box, on the symmetric axis: pz c , r c q " p0.5, 0q. The extension of the clump is pz 0 , r 0 q " p0.1, 0.06q. We consider a homogeneous gas, with ρ 0 " 1 g cm ´3, expect for the clump with density ρ 1 " 100ρ 0 . The boundary of the clump is smoothed:

ρ pz, rq " ρ 0 `ρ1 ´ρ0 1`expp∆q , with ∆ " 10 ˆ´z´zc z 0 ¯2 `´r´rc r 0 ¯2 ´1˙.
The opacity of the medium is σ a " σ s " σ 0 ´T T 0 ¯´3.5 ´ρ ρ 0 ¯2, with σ 0 " 0.1 cm ´1. Initially, the medium is at equilibrium with the radiation, with T 0 " T r " 290 K. At time t " 0, a source is lit at the left boundary with T r " 1 740 K and the reduced flux, aligned with the cylinder, is set to f " 1. Figure 3.7 shows the radiative temperature at the final time t f " 6 ˆ10 ´11 s with different solvers: a classical HLL scheme, the split and unsplit schemes. Because of the high opacity in the clump, the light does not cross it, and we expect the shadow behind it to remain stable.

To compare our results with González et al. [2007]; Hayes and Norman [2003], fig. 3.8 shows the radial profile of radiative temperature at the right boundary. Using both split and unsplit solvers, the shadow is better preserved than using the HLL solver. The solution obtained with these all-regime solvers does not overheat, unlike the one obtained with the HLL solver (see the discussion about the modification of the eigenvalues in section 2.2.4). As discussed in the previous tests, both split and unsplit solvers perform similarly. This all-regime strategy allows keeping numerical diffusion under control, at least when the radiation propagates along the mesh axis.

Discussion and conclusion

In this chapter, we have presented two numerical schemes. The first one is based on a diffusion-transport splitting, similar to the acoustic-transport splitting for Euler equations. However, this leads to a scheme of stencil two, with a high computational cost. The second one is an approximate Riemann solver and relies on a relaxation. Both solvers perform well on classical test cases. Nevertheless, they present some limitations.

Well-balanced discretization of the source terms

Numerical fluxes proposed by Buet and Despres [2008] take into account the source terms and are designed to capture correctly the asymptotic regime. For performance reasons, the source terms are still taken into the center of the cells. As shown in chapter 2, this does not lead to a well-balanced modification of the source terms. It would require computing the source terms at the interfaces of the cells, which calls for the resolution of a linear or nonlinear system, depending on the choice of discretization. However, with source terms taken into account in the numerical fluxes, the steady state with a constant radiative flux in the presence of a discontinuity of opacity is better captured than without source terms in the fluxes.

Time implicit integration

In order to study astrophysical problems, coupling to hydrodynamics is essential. Because of the difference between the timescale of the propagation of the fluid and the photons, solving the radiative transfer requires a time implicit solver. Buet and Despres [2008] show that the schemes proposed in this chapter preserve the admissible states E r ą 0 and f ď 1. Using a Newton-Raphson algorithm to solve the nonlinear system coming from a time-implicit discretization with large time steps might not preserve these admissible states [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF]]. In chapter 4, we explore another strategy to solve this nonlinear system while still preserving the admissible states, based on the work of Pichard [

3.A Equation of evolution for pressure and velocity in the diffusion step

We derive the equations of evolution for the velocity and the pressure. Let us begin with the pressure q. From chapter 1, we have q " c 2 ´β2 3c 2 `β2 E r " 1´χ 2 E r , where χ " 3`4f 2 5`2 a 4´3f 2 is the Eddington factor. Because f " Fr cEr , one has

B t q " Bq BE r B t E r `∇Fr q ¨Bt F r , (3.28)
where ∇ Fr q " ´Bq

BF x r , Bq BF y r , Bq BF z r
¯T and F r " pF x r , F y r , F z r q T . Let us now compute each term individually:

Bq BE r " ´Er 2
Bχ BE r `1 ´χ 2 . (3.29) Because Bχ BE r " Bf 2 BE r dχ df 2 " ´p3c 2 `β2 q f 2 p3c 2 ´β2 q E r , (3.30) one has Bq BE r " c 2 `β2 3c 2 ´β2 . (3.31)
Likewise, one can show that

Bq BF x r " ´F x r 2c 2 E r 3c 2 `β2 3c 2 ´β2 " ´2β x 3c 2 ´β2 , (3.32)
with β " pβ x , β y , β z q T . Similar computations can be done for the x and y coordinates. Therefore,

∇ Fr q " ´2 3c 2 ´β2 β. (3.33)
B t E r and B t F r are given by the diffusion system (eq. 3.6). Therefore,

B t q `4c 2 q 3c 2 ´β2 ∇ ¨β ´c2 ´β2 3c 2 ´β2 β ¨∇q " 0. (3.34)
Using the fact that F r " pE r `qq β, one has

B t β " B t ˆFr E r `q ˙. ( 3

3.B Eigenstructure of the diffusion system

For the sake of simplicity, we focus on the one-dimensional case. We recall that the eigenvalues of the M 1 model are given by [González et al., 2007]:

$ ' ' ' ' ' ' & ' ' ' ' ' ' % λ ´" c n x χ 1 ´aχ 12 `4χ ´4χ 1 f 2 λ 0 " c 2 3χ ´1 2 E r F r λ `" c n x χ 1 `aχ 12 `4χ ´4χ 1 f 2 , (3.37)
where n x is the x-coordinate of n " f f and χ 1 " dχ df " 2f a 4´3f 2 . Let us write eqs. 3.6 and 3.8 as

B t U `A pU q B x U " 0, (3.38)
with U " p , v, e, β, qq T and

A pU q " ¨0 0 0 0 0 0 0 v c 2 0 0 0 4c 2 q c 2 ´β2 β 0 0 0 ´c2 ´β2 3c 2 ´β2 β 3 4q `c2 ´β2 ˘2 3c 2 ´β2 0 0 0 4c 2 q 3c 2 ´β2 ´c2 ´β2 3c 2 ´β2 β ‹ ‹ ‹ ‹ ‹ ‹ ' . (3.39)
0 is eigenvalue with multiplicity 3. Let us write Λ ˘the two remaining eigenvalues. They are the eigenvalues of the 2 ˆ2 matrix

B pU q " ˜´c 2 ´β2 3c 2 ´β2 β 3 4q `c2 ´β2 ˘2 3c 2 ´β2 4c 2 q 3c 2 ´β2 ´c2 ´β2 3c 2 ´β2 β ¸. (3.40) Its characteristic polynomial is Λ 2 ´2β `c2 ´β2 3c 2 ´β2 Λ ´c2 ´β 3c 2 ´β2 " 0. Its roots are Λ ˘" c 2 ´β2
β˘?3c . Using, for example, the symbolic computation Python library SymPy [START_REF] Meurer | Sympy: symbolic computing in python[END_REF], one can check that Λ ˘" λ ˘´λ 0 .

The three eigenvectors associated to the eigenvalue 0 are p1, 0, 0, 0, 0q T , p0, 1, 0, 0, 0q T and p0, 0, 1, 0, 0q T . The field associated to the stationary wave is linearly degenerated. One can check that R `(resp. R ´) is an eigenvector associated to the eigenvalue Λ (resp.

Λ ´), with R

`" ¨?3 `?3c´β 4cq ´c`c ´?3β

˘`? 3c´β c2 ´β2 `?3c´β ˘`3c 2 ´β2 ?3c`β ˘pc 2 ´β2 q ´?3 `c2 ´β2 4cq 1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , R ´" ¨?3 `?3c`β 4cq c `c`?3β ˘`? 3c`β c2 ´β2 `?3c`β ˘`3c 2 ´β2 ?3c´β ˘pc 2 ´β2 q ? 3 `c2 ´β2 4cq 1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . (3.41)
The fields associated to Λ ˘are genuinely nonlinear under the condition |β| ă c. If all the waves were linearly degenerated, the Riemann problem could be solved exactly. Here, we use the result of Buet and Despres [2008].

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . ( 3.47)
After some algebra, the eigenvalues of M are 0 with multiplicity 6 and l ˘with multiplicity 1. Let us notice that the eigenvalues of eqs. 3.44 and 3.45 are equal to the eigenvalues of the diffusion step whenever l ˘" Λ ˘. The eigenvectors associated with the eigenvalue 0 are

¨1 0 0 0 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 1
0 0 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 0
1 0 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 0 0 0 0 1 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 0
0 0 0 0 1 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 0
0 0 0 0 0 1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . (3.48) Ě R `" ¨ al 1 l `` v a `c2 1 l `` e`Π a `b1 a 1 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' and Ě R ´" ¨ al 1 l ´` v a ´c2 1 l ´` e`Π a ´b1 a ´1 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' (3.49)
are the eigenvectors associated with the eigenvalues l ˘.

x x t " 0 pβ L , q L , a L q pβ R , q R , a R q Λ ´Λp β ˚, q ˚, a Lq pβ ˚, q ˚, a Rq

3.C.3 Riemann invariants for a hyperbolic system

Let us recall some notions defined by Bouchut [2004]. We consider a hyperbolic system of the form

B t U `A pU q B x U " 0, (3.50)
where U P R d and A is a d ˆd matrix. Because the system is hyperbolic, A has d real eigenvalues, λ j , and d independent eigenvectors, R j "

´Rp1q j , ¨¨¨, R pdq j ¯T , for j " 1, ¨¨¨, d.
A weak j-Riemann invariant is a scalar function w pU q constant along the characteristic curves. In particular, from eq. 2.123 in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF],

dU pkq R pkq j (3.51)
is constant across the wave λ j , for all k " 1, ¨¨¨, d. For example, if d " 2, U "

`U p1q , U p2q ˘T , R 1 " ´Rp1q 1 , R p2q
1 ¯T and R 2 " ´Rp1q 2 , R p2q
2 ¯T . Then,

dU p1q R p1q 1 " dU p2q R p2q 1 (3.52)
across the wave λ 1 , and

dU p1q R p1q 2 " dU p2q R p2q 2 (3.53)
across the wave λ 2 .

A strong j-Riemann invariant is a scalar function w pU q characterized by B t w pU q `λj pU q B x w pU q " 0.

(3.54)

Bouchut [2004] shows that w pU q is a strong j-Riemann invariant if and only if w pU q is a weak k-Riemann invariant, for k " 1, ¨¨¨, j ´1, j `1, ¨¨¨, d.

3.C.4 Riemann invariants for the relaxed system

Figure 3.9 shows the wave structure for the system made of eqs.

4c ? 3 Π R c 2 ´b2 R pb R
´b˚q ´pΠ R ´Π˚q " 0.

(3.57) Equations 3.56 and 3.57 result in an exact solver for the relaxed system made of eqs. 3.44 and 3.45. We obtain an approximate solver for the original system, eq. 3.42. Using q " c 2 ´β2 3c 2 `β2 E, we recover the same jump relations as Buet and Despres [2008] (see also eq. 3.24). One can show that β ˚and q ˚given by eq. 3.13 are the solution of the system made of eqs. 3.56 and 3.57.

Let us also recall eq. 1.82: (3.67) From the jump relations eq. 3.24, one can show that q i`1 2 ptq ´qi ptq and β i ptq ´βi `1 2 ptq have the same sign, while q i´1 2 ptq ´qi ptq and β i ptq ´βi

B t S r `Bx pβS r q ě 0, ( 3
´1 2
ptq have the opposite sign. One can deduce that ∆m i θ i ptq ds i dt ptq is non negative. Because ∆m i and θ i ptq are non negative, one has ds i dt ptq ě 0 and the semi-discrete scheme for the diffusion step is entropic.

Let us notice that, because F r " pE r `qq β and q " c 2 ´β2 3c 2 `β2 E r , then

1 4 r ˆ3 3c 2 `β2 ˙3 4 `c2 ´β2 ˘1 4 E 3 4
r .

(3.68)

We assume that |β i | ă c at initial time t " 0. Then, s i pt " 0q ą 0. Because s is non-decreasing, s i ptq ą 0 for all t. This leads to

E 3 4 i `c2 ´β2 i ˘1 4 ą 0. (3.69) If E i ą 0, then |β i | ă c.
Buet and Despres [2008] show that

β i`1 2 ď c, (3.70)
under the Courant-Friedrichs-Lewy (CFL) condition ∆t ď ∆x 2c for the diffusion step.

3.D.2 Transport step

We consider now the transport step. We define rβ ˚s˘" β ˚˘|β ˚|

2

. If β ˚ě 0, then rβ ˚s`" β ˚and rβ ˚s´" 0. On the contrary, if β ˚ď 0, then rβ ˚s`" 0 and rβ ˚s´" β ˚.

With this notation, ϑ

i`1 2 β i`1 2 " ϑ i " β i`1 2 ı ``ϑ i`1 " β i`1 2 ı
´, for ϑ P t , F, Eu. After some algebra, eq. 3.16 becomes

ϑ n`1 i " ˆ1 `∆t ∆x " β i`1 2 ı ´´∆t ∆x " β i´1 2 ı `˙ϑ i ptq´∆ t ∆x " β i`1 2 ı ´ϑi`1 ptq`∆ t ∆x " β i´1 2 ı `ϑi´1 ptq.
(3.71)

3.D.3 Split scheme

Let us consider the complete split scheme, the diffusion step followed by the transport step. Equation 3.71 becomes

ϑ n`1 i " ˆ1 `∆t ∆x " β i`1 2 ı ´´∆t ∆x " β i´1 2 ı `˙r ϑ i ´∆t ∆x " β i`1 2 ı ´r ϑ i`1 `∆t ∆x " β i´1 2 ı `r ϑ i´1 .
(3.72) Under the CFL condition eq. 3.70, ϑ n`1 i is a convex combination of r ϑ i´1 , r ϑ i and r ϑ i`1 . Similar computations can be done for the entropy, with inequality in eq. 3.72. By doing the same computations as for the diffusion step, we show that E n`1 i ą 0 and β n`1 i ď c.

3.D.4 Unsplit scheme

The unsplit scheme can be seen as the average of two steps:

ϑ n`1 i " 1 2
`ϑD i `ϑT i ˘, (3.73)
with ϑ D i given by the diffusion step and

ϑ T i " ϑ n i ´∆t ∆x ´βi `1 2 ϑ n i`1 2 ´βi ´1 2 ϑ n i´1 2 ¯, (3.74)
for ϑ P t , F, Eu. Using eq. 3.71, we have

T i " ˆ1 `∆t ∆x " β i`1 2 ı ´´∆t ∆x " β i´1 2 ı `˙ n i ´∆t ∆x " β i`1 2 ı ´ n i`1 `∆t ∆x " β i´1 2 ı ` n i´1 .
(3.75) Under the CFL condition eq. 3.70, n`1 i is a convex combination of n i´1 , n i and n i`1 . Therefore, T i ě 0. For θ P tv, eu, we have

θT i " ` θ˘T i T i " n i T i ˆ1 `∆t ∆x " β i`1 2 ı `´∆t ∆x " β i´1 2 ı ´˙θ n i ´ n i`1 T i ∆t ∆x " β i`1 2 ı ´θ n i`1 ` n i´1 T i ∆t ∆x " β i´1 2 ı `θ n i´1 " α 0 i T i θn i `αì T i θn i`1 `αí T i θn i´1 .
(3.76)

3.E.2 Eigenstructure

Let us now study the eigenstructure of eq. 3.79. We write Ă W " p , v, e, r, µ, ϕ, b, Πq T , then eq. 3.79 can be written in the matrix form

B t Ă W `r M ´Ă W ¯Bx Ă W " 0, (3.80) where r M ´Ă W ¯" ¨0 0 0 b 0 0 r 0 0 0 0 0 b 0 µ c 2 0 0 0 0 0 b ϕ `Π b 0 0 0 b 0 0 0 0 0 0 0 0 b 0 0 0 0 0 0 0 0 b 0 0 0 0 0 0 0 0 ´c2 ´b2 3c 2 ´b2 b 3 4Π `c2 ´b2 ˘2 3c 2 ´b2
0 0 0 0 0 0

4c 2 Π 3c 2 ´b2 ´c2 ´b2 3c 2 ´b2 b ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . (3.81)
Using similar computations as in section 3.B, the eigenvalues of r M are 0 with multiplicity 3, b with multiplicity 3, r Λ `and r Λ ´, with r Λ ˘" c 2 ´b2 b˘?3c . The three eigenvectors associated with the eigenvalue 0 are ¨1 0

0 0 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 1
0 0 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 0
1 0 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . (3.82)
The field associated with the stationary wave is linearly degenerated. The three eigenvectors associated with the eigenvalue b are ¨1 0

0 1 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 1 0 0 1 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , ¨0 0
1 0 0 1 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . (3.83)
The field is also linearly degenerated. Finally, r R `and r R ´are the eigenvectors associated with the eigenvalues r Λ `and r Λ ´respectively, with r R `" ¨?3r `?3c´b 4cΠ ´c`c ´?3b x

˘`? 3c´b c2 ´b `?3c´b ˘`3c 2 ´b2 ?3c`b ˘pc 2 ´b2 q 0 0 0 ´?3 `c2 ´b2 4cΠ 1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , r R ´" ¨?3r `?3c`b 4cΠ c `c`?3b ˘`? 3c`b c2 ´b2 `?3c`b ˘`3c 2 ´b2 ?3c´b ˘pc 2 ´b2 q 0 0 0 ? 3 `c2 ´b2 4cΠ 1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . ( 3
P R Contact wave D " 0, b ‰ 0 " 0 - D ‰ 0, rbs " 0, D " b P R Contact wave D ‰ 0, rbs " 0, D ‰ b " 0 - D ‰ 0, rbs ‰ 0, D ‰ b " 0 - D ‰ 0, rbs ‰ 0, D " b P R Shock wave
x x t " 0 pb L , Π L q pb R , Π R q D " r Λ ´D " r Λ pb ˚, Π ˚q
x t " 0 pb L , Π L q pr L , µ L , ϕ L q pb R , Π R q pr R , µ R , ϕ R q D " r Λ ´D " r Λ pb ˚, Π ˚q pb L , Π L q D " b L pr R , µ R , ϕ R q pr R , µ R , ϕ R q
x x t " 0 pb L , Π L q pr L , µ L , ϕ L q pb R , Π R q pr R , µ R , ϕ R q D " r Λ ´D " r Λ pb ˚, Π ˚q pb ˚, Π ˚q D " b pr L , µ L , ϕ L q pr R , µ R , ϕ R q Figure 3
.12 -Solution of the Riemann problem for b, Π, r, µ, ϕ with D " b ˚ă 0.

x x t " 0 pb L , Π L q pr L , µ L , ϕ L q pb R , Π R q pr R , µ R , ϕ R q D " r Λ ´D " r Λ pb ˚, Π ˚q pb ˚, Π ˚q D " b pr L , µ L , ϕ L q pr R , µ R , ϕ R q Figure 3
.13 -Solution of the Riemann problem for b, Π, r, µ, ϕ with D " b ˚ą 0. We now build an approximate Riemann solver for eq. 3.79. We are looking for a function composed of three states separated by discontinuities as follows (see fig. 3.10):

x x t " 0 pb L , Π L q pr L , µ L , ϕ L q pb R , Π R q pr R , µ R , ϕ R q D " r Λ ´D " r Λ pb ˚, Π ˚q pb R , Π R q D " b R pr L , µ L , ϕ L q pr L , µ L , ϕ L q
$ ' ' ' ' & ' ' ' ' % pr L , µ L , ϕ L , b L , Π L q T if x t ă r Λ ´, pr ˚, µ ˚, ϕ ˚, b ˚, Π ˚qT if r Λ ´ă x t ă r Λ `, pr R , µ R , ϕ R , b R , Π R q T if x t ą r Λ `.
(3.93)

The computation of b ˚and Π ˚is independent of the other variables. We take b ˚and Π ˚given by eq. 3.26, but other choices can be done. Let us now show the different possible wave patterns. Then, we distinguish the different cases corresponding to the localization of the wave D " b:

-D " b L : we take pr ˚, µ ˚, ϕ ˚qT " pr R , µ R , ϕ R q T , see fig. 3.11. -D " b ˚: we take pr ˚, µ ˚, ϕ ˚qT " $ & % pr L , µ L , ϕ L q T if x t ă b ˚, pr R , µ R , ϕ R q T if x t ą b ˚,
(3.94) see figs. 3.12 and 3.13. -D " b R : we take pr ˚, µ ˚, ϕ ˚qT " pr L , µ L , ϕ L q T , see fig. 3.14.

3.E.5 Riemann flux computation

From the wave pattern shown previously, we can deduce the following algorithm to compute the numerical fluxes F ˚on the density and radiative flux variables, and v given by our relaxed system eq. 3.79. To do so, we evaluate the flux function from eq. 3.5 on the state obtained on x t " 0:

1. Compute pb ˚, Π ˚q with eq. 3.26; 2. Compute the flux according to the sign of b ˚:

F ˚¨¨ L F L E L ' , ¨ R F R E R '' " $ ' ' ' ' ' ' & ' ' ' ' ' ' % ¨ L b FL b ˚`c 2 Π EL b ˚`Π ˚b˚' if b ˚ě 0 ¨ R b FR b ˚`c 2 Π ER b ˚`Π ˚b˚' if b ˚ď 0.
(3.95)

3.F Numerical schemes in the diffusive limit

As in section 2.B, we consider both numerical schemes developed in section 3.1 in the asymptotic regime, with σ s i`1 2 ∆x Ñ 8. Following section 1.4, we introduce the rescaling parameter ε to write the time (resp. the opacity) as Ă ∆t " ε∆t (resp. r σ a " εσ a and r σ s " εσ s ). This leads to

ε n`1 i " ε n i ´Ă ∆t ∆x ´ m i`1 2 β i`1 2 ´ m i´1 2 β i´1 2 ¯(3.96a) ε 2 F n`1 i " ε 2 F n i ´ε Ă ∆t ∆x ´F m i`1 2 β i`1 2 `c2 q i`1 2 ´F m i´1 2 β i´1 2 ´c2 q i´1 2 c r σ s i Ă ∆tF n`1 i (3.96b) ε 2 E n`1 i " ε 2 E n i ´ε Ă ∆t ∆x ´Em i`1 2 β i`1 2 `qi `1 2 β i`1 2 ´Em i´1 2 β i´1 2 ´qi ´1 2 β i´1 2 c r σ a i Ă ∆t ´ar `T n`1 i ˘4 ´En`1 i ¯(3.96c) ε 2 ρc v T n`1 i " ε 2 ρc v T n i ´c r σ a i Ă ∆t ´ar `T n`1 i ˘4 ´En`1 i ¯, (3.96d)
where ϑ m i`1 2 is given by eq. 3.17 for the split scheme (ϑ m i`1 2 " r ϑ i`1 2 ) and by eq. 3.20 for the unsplit scheme (ϑ

m i`1 2 " ϑ n i`1 2
), for ϑ P t , F, Eu.

Variables are expanded, e.g., E n i " E n i,0 `εE n i,1 `O pε 2 q. Expanding eq. 3.96c or eq. 3.96d at order 0 leads to

E n`1 i,0 " a r `T n`1 i,0 ˘4 . ( 3
.97) Equation 3.96b at order 0 leads to

F n`1 i,0 " 0. (3.98) Because F n i,0 " `En i,0
`qn i,0 ˘βn i,0 and E `q ą 0, one has

β n i,0 " 0. (3.99)
At order 0 for eq. 3.96a, we have β i`1 2 ,0 ´βi ´1 2 ,0 " 0 for all cells i. Therefore, β i`1 2 ,0 " C, where C is a constant that does not depend on i. With well-chosen boundary conditions, we have β 0`1 2 ,0 " 0, which leads to β i`1 2 ,0 " 0 for all cells i. Moreover,

0 " β i`1 2 ,0 " ? 3c 4 ˆqn i,0 ´Ă σ s i`1 2 2c F n i,1 ˙´ˆq n i`1,0 `Ă σ s i`1 2 2c F n i`1,0 Ėn i,0 `En i`1,0 ô F n i`1,1 `F n i,1 2 " ´c r σ s i`1 2 q n i`1,0 ´qn i,0
∆x .

(3.100)

Using eq. 3.99 in the definition of q, we have q n i,0 " 1 3 E n i,0 and

F n i`1,1 `F n i,1 2 " ´c 3 r σ s i`1 2 E n i`1,0 ´En i,0 ∆x . ( 3

.101)

Let us now consider the sum of eq. 3.96c and eq. 3.96d. The term of order 1 is 0 because β i`1 2 ,0 " 0. We now look at the term of order 2 in the numerical fluxes. For the same reason, it is

E m i`1 2 ,0 β i`1 2 ,1 `qi `1 2 ,0 β i`1 2 ,1 ´Em i´1 2 ,0 β i´1 2 ,1 ´qi ´1 2 ,0 β i´1 2 ,1 . (3.102)
From now on, we assume that all functions are sufficiently smooth to write, for example,

E n i`1,0 " E n i,0 `∆x E n i`1,0 ´En i,0 ∆x " E n i,0 `O p∆xq . (3.103)
We can now look at the different terms involved in eq. 3.102. First,

β i`1 2 ,1 " E n i,0 β n i,1 `En i`1,0 β n i`1,1 E n i,0 `En i`1,0 ´3? 3c 4 ¨qn i`1,1 ´qn i,1 E n i,0 `En i`1,0 loooooomoooooon Op∆xq ´r σ s i`1 2 ∆x 2c F n i`1,2 ´F n i,2 E n i,0 `En i`1,0 loooooomoooooon Op∆xq ‹ ‹ ‹ ‹ ' " α i β n i,1 `p1 ´αi q β n i`1,1 `O p∆xq , (3.104) with α i " E n i,0 E n i,0 `En i`1,0 . Second, q i`1 2 ,0 " q n i,0 E n i,0 `qn i`1,0 E n i`1,0 1 E n i,0 `1 E n i`1,0 `r σ s i`1 2 ∆x 2c F n i`1,1 E n i`1,0 ´Fn i,1 E n i,0 1 E n i,0 `1 E n i`1,0 " E n i`1,0 q n i,0 `En i,0 q n i`1,0 E n i,0 `En i`1,0
`O p∆xq " p1 ´αi q q n i,0 `αi q n i`1,0 `O p∆xq .

(3.105)

Using eqs. 3.104 and 3.105, one has

β n i`1 2 ,1 q n i`1 2 ,0 " α i β n i,1 q n i,0 `p1 ´αi q β n i`1,1 q n i`1,0 `O p∆xq . (3.106)
Furthermore,

α i " E n i,0 E n i,0 `En i`1,0 " 1 2 `O p∆xq , (3.107) which leads to β n i`1 2 ,1 q n i`1 2 ,0 " 1 2 β n i,1 q n i,0 `1 2 β n i`1,1 q n i`1,0 . (3.108)
Likewise,

E m i`1 2 ,0 β n i`1 2 ,1 " α i E m i`1 2 ,0 β n i,1 `p1 ´αi q E m i`1 2 ,0 β n i`1,1 `O p∆xq " 1 2 E m i,0 β n i,1 `1 2 E m i`1,0 β n i`1,1 `O p∆xq . (3.109) Finally, ´Em i`1 2 ,0 `qn i`1 2 ,0 ¯βn i`1 2 ,1 " 1 2 `Em i,0 `qn i,0 ˘βn i,1 `1 2 `Em i`1,0 `qn i`1,0 ˘βn i`1,1 `O p∆xq .
(3.110) Using the unsplit scheme, we have E m

i " E n i and eq. 3.110 becomes

´Em i`1 2 ,0 `qn i`1 2 ,0 ¯βn i`1 2 ,1 " F n i,1 `F n i`1,1 2 `O p∆xq , (3.111)
which leads to

E n`1 i,0 `ρc v T n`1 i,0 " E n i,0 `ρc v T n i,0 `cĂ ∆t 3∆x 2 ˜En i`1,0 ´En i,0 r σ s i`1 2 ´En i,0 ´En i´1,0 r σ s i´1 2 ¸`O p∆xq .
(3.112) Using the split scheme, we can write E m i " r E i " E n i `O p∆tq and eq. 3.110 becomes

E n`1 i,0 `ρc v T n`1 i,0 " E n i,0 `ρc v T n i,0 `cĂ ∆t 3∆x 2 ˜En i`1,0 ´En i,0 r σ s i`1 2 ´En i,0 ´En i´1,0 r σ s i´1 2 ¸`O p∆xq`O p∆tq .
(3.113) Equations 3.112 and 3.113 are consistent with the diffusion equation eq. 1.24 because O p∆xq and O p∆tq go to 0 as ∆x and ∆t go to 0. To summarize, eqs. 3.97, 3.98, 3.101, 3.112 and 3.113 are standard discretization of eqs. 1.22 to 1.24, so these schemes are asymptotic preserving.

Unlike the scheme developed in chapter 2 or the one presented by Berthon and Turpault [2011], we do not need to choose a parameter to recover the asymptotic behavior.

In chapter 2, we presented an asymptotic preserving scheme, with a well-balanced modification of the source terms, that uses a time-implicit integration. Nevertheless, it does not preserve the admissible states E r ą 0 and ||F r || ď cE r . To tackle this issue, we presented in chapter 3 an asymptotic preserving scheme that preserves the admissible states. For the sake of simplicity, we derived it using a time-explicit integration. This chapter is an additional step towards the derivation of a time-implicit and asymptotic preserving scheme that preserves the admissible states: we present an implicit solver that preserves the admissible states. In this chapter, we consider a time-implicit HLL solver, without source terms. Therefore, the asymptotic preserving and well-balanced properties are not considered.

This chapter is organized as follows. We first present in section 4.1 a Jacobi method to solve the nonlinear system arising from the discretization of the M 1 model. Because this method is iterative, its convergence rate can be improved thanks to multigrid acceleration. We explore this technique in section 4.2. In section 4.3, we perform some tests to validate both algorithms. We also show some performance results. Finally, we reach our conclusion in section 4.4.

Nonlinear Jacobi method

We use notations similar to the previous chapters. We note h the step along the x-direction. As previously, ∆t is the time interval between the current time t n and t n`1 . We write x i the center of the cell i. We use the notation u n i to represent the averaged quantity associated with the field u at time t n in the cell i (finite volume).

Time-implicit HLL solver

Let us consider a time-implicit HLL solver for the M 1 model without source terms. For the sake of simplicity, we use fixed eigenvalues in the numerical fluxes (eq. 2.2). This leads to solving the following nonlinear system:

$ ' ' & ' ' % E n`1 i ˆ1 `c∆t h ˙´∆t 2h `cE n`1 i`1 ´F n`1 i`1 ˘´∆t 2h `cE n`1 i´1 `F n`1 i´1 ˘" E n i F n`1 i ˆ1 `c∆t h ˙´c∆t 2h `F n`1 i`1 ´cP n`1 i`1 ˘´c∆t 2h `F n`1 i´1 `cP n`1 i´1 ˘" F n i . (4.1)
As done in chapter 2, eq. 4.1 can be solved using a Newton-Raphson method. However, numerical tests, especially in the free-streaming regime, have shown that the admissible states are not preserved when large time steps are used. A solution is therefore to reduce the time step, but it leads to poor performances when the radiative transfer is coupled to hydrodynamics. [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] show how the resolution of eq. 4.1 can be reduced to the resolution of a linear system (by fixing the nonlinearity at time t n ) and why this approach does not preserve the admissible states either.

In the next section, we present another method to solve eq. 4.1 while preserving the admissible states.

Jacobi method

Let us follow the work of [START_REF] Pichard | Martin Mathematiques appliquees et calcul scientifique Bordeaux[END_REF]. We first define the set of admissible states test already described in sections 2.2.3 and 3.2.1.2. It is the propagation of a beam in the vacuum, therefore it is well-suited to test our solver, as there is no opacity. A square two-dimensional domain is discretized with the same number of cells in both x and y-directions. Because the solver is time-implicit, we set a time step ∆t such that the steady state is reached with one iteration. This test is discussed more precisely in section 4.3.1. Figure 4.1 shows the evolution of the norm of the residual as a function of the number of iterations of the Jacobi method, with different resolutions. As the resolution increases, the number of iterations needed to reach the same residual also increases, from 4 000 iterations with 129 ˆ129 cells up to 15 000 iterations with 513 ˆ513 cells. As shown by table 4.1, performances decrease as the resolution increases, from around 150 cell-updates/s with 129 ˆ129 cells down to 40 cell-updates/s with 513 ˆ513 cells. Indeed, Saad [2003] explains that the convergence rate decreases as the size of the problem increases, resulting in a decrease in performance. For linear problems, classical iterative methods are efficient to compute the high frequencies of the solution, but lack efficiency to compute its low frequencies. However, the computation is easier on a coarser grid with fewer unknowns. This observation has lead to the development of the geometric multigrid technique. The initial guess for the Jacobi algorithm (or any iterative method) is an interpolation of a solution computed on a coarser grid. This method requires information about the geometry of the problem, unlike preconditioners based on Krylov subspace. This additional information allows the multigrid method to be very efficient, but it lacks generality. To tackle this issue, algebraic multigrid methods have been developed, see section 5.1.5.

S " tpE r , F r q , E r ą 0, ||F r || ď cE r u . ( 4
We first present the geometric multigrid algorithm in the linear case (section 4.2.2). Then, in section 4.2.3, we present its extension to the general nonlinear case, leading to the Full Approximation Scheme (FAS). Finally, in section 4.2.4 we apply it to the M 1 model, and we highlight the distinctive features of radiative transfer, i.e., the preservation of the admissible states E r ą 0 and f ď 1.

Linear case

In this section, we present the main ideas of the geometric multigrid method. It does not intend to be a full review of existing work, see e.g., [START_REF] Briggs | A Multigrid Tutorial[END_REF]; Brandt and Livne [2011].

We are interested in solving the linear system

Av " f , (4.6)
with an iterative method, such as Jacobi [Saad, 2003], where A is an invertible matrix. We assume that this system is obtained by the discretization of a PDE and therefore depends on a step h. For example, let us consider the Poisson equation:

# ´∆u pxq " f pxq for x P Ω, u pxq " 0 for x P BΩ, (4.7)
where BΩ is the boundary of Ω. In the one-dimensional case, a finite difference discretization of eq. 4.7 leads to solving a linear system of the form given by eq. 4.6, with

A " 1 h 2 ¨2 ´1 ´1 2 ´1 . . . . . . . . . ´1 2 ´1 ´1 2 ‹ ‹ ‹ ‹ ‹ ' and f " ¨f px 1 q . . . f px n q ‹ ' . (4.8)
As mentioned above, the initial guess used in the Jacobi algorithm is the interpolation of a solution obtained on a coarser grid. This solution computed on the coarse grid can be itself computed thanks to a solution obtained on a third grid, even coarser. This process can be applied recursively until the last grid has only a few unknowns and the problem could be solved with a direct method, which leads to the so-called "nested iterations". From Briggs et al. [2000], this requires the definition of two elements:

-a hierarchy of grids, completed with restriction and prolongation operators; -a smoother to solve the system on a given mesh. Our goal is to solve nonlinear systems coming from the discretization of the M 1 model, then we will use the Jacobi method presented in section 4.1 as smoother. Therefore, we also use a Jacobi method as smoother in the linear case, and we do not investigate the impact of the smoother anymore. In the linear case, one can study the eigenvalues and eigenvectors of A, as well as the iterative method to understand how the smoother allows the quick decrease of the high frequencies of the error.

The hierarchy of grids is handled through recursiveness, we only need to describe the method with two grids. The domain Ω is discretized with two Cartesian meshes: the first one with step h, written Ω h and the second one with step 2h, written Ω 2h . Ω h contains 2 d times more elements than Ω 2h , where d is the dimension of the problem.

h

Ω h 2h Ω 2h P h 2h R 2h h Figure 4.2 -Two-grid cycle. 1 2 3 4 5 6 7 Ω h 1 2 3 Ω 2h 1 1 2 1 2 Figure 4.
3 -Prolongation operator in the one-dimensional case. The blue arrow represents the operator if i h is even and the red arrows represent the operator if i h is odd.

i h , j h i 2h , j 2h 1 i h , j h i 2h , j 2h 1 2 i 2h `1, j 2h 1 2 i h , j h i 2h , j 2h 1 2 i 2h , j 2h `1 1 2 i h , j h i 2h , j 2h 1 4 i 2h `1, j 2h 1 4 i 2h , j 2h `1 1 4 i 2h `1, j 2h `1 1 4
Figure 4.4 -Prolongation operator in the two-dimensional case. The coarse mesh Ω 2h is in blue, the fine mesh Ω h is in red, the weight of the cells in the coarse mesh is in green. The computed cell is colored.

Let us now define the restriction and prolongation operators used for inter-grid operations. We write the prolongation operator P h 2h : Ω 2h Ñ Ω h and the restriction operator R 2h h : Ω h Ñ Ω 2h , see fig. 4.2. We begin with the prolongation operator. We use full weighting operators [START_REF] Strang | 086 mathematical methods for engineers ii[END_REF]. Let us consider the one-dimensional case, with an odd number of cells in the fine mesh. Indices i h even in the fine mesh are the same as i 2h in the coarse mesh (see fig. 4.3). Values corresponding to these cells in the coarse mesh are just moved in the fine mesh. The other values in the fine mesh are obtained by linear interpolation. The prolongation operator writes

`Ph 2h `v2h ˘˘i h " $ & % v 2h i 2h if i h is even, 1 2 `v2h i 2h `v2h i 2h `1˘i f i h is odd. (4.9)
We present now the prolongation operator in the two-dimensional case. A cell i h , j h in the fine mesh Ω h can be mapped directly into the cell i 2h , j 2h in the coarse mesh Ω 2h , with i h " 2i 2h and j h " 2j 2h . The value in a cell of the fine mesh is obtained by the interpolation in one direction, followed by the interpolation in the other direction. The interpolation in the x-direction gives

`Ph 2h `v2h ˘˘2i 2h ,2j 2h " v 2h i 2h ,j 2h and `Ph 2h `v2h ˘˘2i 2h `1,2j 2h " 1 2 `v2h i 2h ,j 2h `v2h i 2h `1,j 2h
˘.

(4.10) The interpolation in the y-direction preserves these values and gives

`Ph 2h `v2h ˘˘2i 2h ,2j 2h `1 " 1 2 `v2h i 2h ,j 2h `v2h i 2h ,j 2h `1˘. (4.11)
Finally, the last value is given by

`Ph 2h `v2h ˘˘2i 2h `1,2j 2h `1 " 1 4 `v2h i 2h ,j 2h `v2h i`1 2h ,j 2h `v2h i 2h ,j 2h `1 `v2h i 2h `1,j 2h `1˘( 4.12) i 2h , j 2h i h , j h 1 16 1 8 1 16 1 8 1 4 1 8 1 16 1 8 1 16
Figure 4.5 -Restriction operator. The coarse mesh Ω 2h is in blue, the fine mesh Ω h is in red, the weight of the cells in the fine mesh is in green. The computed cell is colored.

Using eqs. 4.10 to 4.12 leads to writing the prolongation operator as

`Ph 2h `v2h ˘˘i h ,j h " $ ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' % v 2h i 2h ,j 2h if i h is even and j h is even, 1 2 `v2h i 2h ,j 2h `v2h i 2h `1,j 2h ˘if i h is odd and j h is even, 1 2 `v2h i 2h ,j 2h `v2h i 2h ,j`1 2h ˘if i h is even and j h is odd, 1 4 `v2h i 2h ,j 2h `v2h i`1 2h ,j 2h `v2h i 2h ,j 2h `1 `v2h i 2h `1,j 2h `1ȋ
f i 2h is odd and j 2h is odd. Let us now consider the restriction operator. We first consider the case where the mesh is given by fig. 4.3. The prolongation operator can be made explicit, with appropriate boundary conditions:

v h " `Ph 2h `v2h ˘˘" ¨vh 1 v h 2 v h 3 v h 4 v h 5 v h 6 v h 7 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " ¨vh 2ˆ0`1 v h 2ˆ1 v h 2ˆ1`1 v h 2ˆ2 v h 2ˆ2`1 v h 2ˆ3 v h 2ˆ3`1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " 1 2
¨1 0 0 2 0 0 1 1 0 0 2 0 0 1 1 0 0 2 0 0 1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' ¨v2h 1 v 2h 2 v 2h 3 ' . (4.14)
Even though the prolongation operator is written as an operator, it is linear and can therefore be seen as a matrix, P h 2h . Apart the first and the last rows that handle of boundary conditions, the sum of the coefficients in a row is 1. Values on the fine mesh are convex combinations of values on the coarse mesh.

We take the restriction operator as the transpose of the prolongation operator:

R 2h h " 1 writes v 2h " `R2h h `v2h ˘˘" ¨v2h 1 v 2h 2 v 2h 3 ' " 1 4 ¨1 2 1 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 1 2 1 ' ¨vh 1 v h 2 v h 3 v h 4 v h 5 v h 6 v h 7 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . ( 4.15)
In the general one-dimensional case, we have

`R2h h `vh ˘˘i 2h " 1 4 v h i h ´1 `1 2 v h i h `1 4 v h i h `1. (4.16)
Numerical tests in section 4.3 are done in the two-dimensional case, we also make R 2h h explicit with d " 2 (see also fig. 4.5):

`R2h h `vh ˘˘i 2h ,j 2h " 1 16 v h i h ´1,j h ´1 `1 8 v h i h ´1,j h `1 16 v h i h ´1,j h `1 `1 8 v h i h ,j h ´1 `1 4 v h i h ,j h `1 8 v h i h ,j h `1 `1 16 v h i h `1,j h ´1 `1 8 v h i h `1,j h `1 16 v h i h `1,j h `1.
(4.17)

Let us recall that we are interested in solving eq. 4.6 on the mesh Ω h . We then rewrite it as

A h v h " f h , (4.18)
to highlight the dependence on the mesh. To fix the notations, the discretization of eq. 4.7 on the coarse mesh Ω 2h leads to solving the following linear system:

A 2h v 2h " f 2h . (4.19)
Now that we have made explicit the restriction and prolongation operators used with two grids, we write the two-grid algorithm (algorithm 2).

Algorithm 2 Two-grid algorithm

Pre-smoother: relax ν 0 times A h u h " f h with initial guess v h on the fine mesh

Ω h Restrict the current solution v 2h " R 2h h `uh ˘and the right-hand side f 2h " R 2h h `f h Solve A 2h u 2h " f 2h with initial guess v 2h on the coarse mesh Ω 2h Prolong the solution: vh " P h 2h `u2h
Post-smoother: relax ν 0 times A h u h " f h with initial guess vh on the fine mesh Ω h This process can now be applied recursively until a mesh coarse enough is reached. This leads to the so-called V-cycle algorithm. Figure 4.6 shows the mesh size with different levels. Solving the linear system at a given level is less costly than the previous level because the number of unknowns is reduced. However, numerical experiments have shown that the best performances are obtained by performing a series of V-cycles. The algorithm is made explicit by algorithm 3. At the coarsest level, a direct solver could be used. In the case of radiative transfer, the development of a direct solver is out of reach. Therefore, we will also use the Jacobi method as a coarse grid solver.

Algorithm 3 V-cycle

Choose a maximum number of levels

L max κ " 0 while A h v h pκq ´f h ą ε do v h pκ`1q " V-cycle(0, v h pκq , f h , L max ) κ Ð κ `1 end while function V-cycle(l, v 2 l h , f 2 l h , L max ) Pre-smoother: relax ν l times A 2 l h u 2 l h " f 2 l h with initial guess v 2 l h Restrict the current solution v 2 l`1 h " R 2 l`1 h 2 l h
´u2 l h ¯and the right-hand side

f 2 l`1 h " R 2 l`1 h 2 l h ´f 2 l h īf l `1 " L max then Solve A 2 l`1 h u 2 l`1 h " f 2 l`1 h with initial guess v 2 l`1 h else Recursion: u 2 l`1 h "V-cycle( l `1, v 2 l`1 h , f 2 l`1 h , L max ) end if Prolong the solution: v2 l h " P 2 l h 2 l`1 h ´u2 l`1 h
Post-smoother: Because we mimic this in the linear case, we do not investigate the effect of the coarse grid solver. Now that we have presented the algorithm, let us consider a simple test. We solve eq. 4.7 in the two-dimensional case, with Ω " r0, 1s ˆr0, 1s and f px, yq " 8π 2 sinp2πxq sinp2πyq. Equation 4.7 has an exact solution: ūpx, yq " sinp2πxq sinp2πyq.

relax ν l times A 2 l h u 2 l h " f 2 l h with initial guess v2 l h return u 2 l h end function 200 400 600 800 10 ´3 10 ´2 10 ´1 Number of V-cycles (κ) Residual r h pκq r h p0q L max " 1 L max " 2 L max " 3 L max " 4 L max " 5
We discretize Ω with 129 ˆ129 cells and apply algorithm 3. Figure 4.7 shows the evolution of the residual A h v h ´f h as a function of the number of V-cycles, with different numbers of levels, L max . With L max " 1, we recover the classical Jacobi method. Apart from L max " 1, the residual drops quickly, as L max increases, but it reaches a plateau, around 10 -1 , depending on the value of L max . The same phenomena have been observed with different resolutions.

Algorithm 4 Two-grid correction scheme

Pre-smoother: relax ν 0 times A h u h " f h with initial guess v h on the fine mesh Ω h Restrict the residual r 2h " R 2h h `f h ´Ah u h Solve A 2h e 2h " r 2h on the coarse grid with initial guess 0 on the coarse mesh Ω 2h Correct the fine grid approximation vh " v h `Ph 2h `e2h

Post-smoother: relax ν 0 times A h u h " f h with initial guess vh on the fine mesh Ω h One can show that the presence of smooth modes in the error between the approximate and exact solution ū can cause the relaxation to stall. The proof relies on the study of the eigenvalues and eigenvectors of A, see e.g., Brandt and Livne 2011. As shown by fig. 4.8, a smooth function on the fine mesh will be less smooth on the coarse mesh. An iterative method as Jacobi will converge faster and will smooth the error on the coarse mesh. This gives the idea of restrict and prolong the error between the different meshes. If v is an approximate solution of eq. 4.6 and u is the exact solution, let us define the error as e " u ´v. It verifies Ae " Au ´Av " f ´Av. By defining the residual r " f ´Av, the error verifies the following equation: Residual We will solve the residual equation eq. 4.20 on the coarse mesh. Once the error is computed, a new approximate solution can be computed: v`e. We replace algorithm 2 with algorithm 4, which leads to the so-called "correction scheme".

Ae " r. ( 4
r h pκq r h p0q L max " 1 L max " 2 L max " 3 L max " 4 L max " 5
Figure 4.9 shows the evolution of the norm of residual as a function of the number of V-cycles, using the correction scheme. Unlike fig. 4.7, the residual drops quickly to 0. It decreases faster as L max increases, as expected. When the residual is small enough, around 10 -4 , all curves converge with the same slope. The low resolution at the coarsest meshes is responsible for that.

We have first presented the algorithm without the error and the residual because it would have been easier to adapt for radiative transfer. As shown in section 4.1.3, if the initial guess and the right-hand side are admissible, then the solution obtained with the Jacobi method is also admissible. The restriction and prolongation operators used here (eqs. 4.13 and 4.17) are convex combinations of admissible states, therefore the restricted and prolonged vectors are also admissible. Unfortunately, as shown here, this method does not converge. Because the error and the residual are not admissible in the general case, the algorithm has to be modified for the M 1 model (section 4.2.4).

General nonlinear case

Now that we have presented the multigrid method in the linear case, let us focus on the nonlinear case. Instead of solving eq. 4.6, we solve

A pvq " f , (4.21)
where A is a nonlinear operator. We suppose that we have an iterative method to solve eq. 4.21, such as nonlinear Jacobi or Gauss-Seidel.

To mimic the linear case, we want to use the residual equation on the coarse mesh to approximate the error on the fine mesh e h . Let us recall that the error is defined as e " u ´v, where u is the exact solution of eq. 4.21 and v is an approximate solution. We also define the residual as r " f ´A pvq. The residual equation eq. 4.20 is no longer verified. The definition of the residual leads to r " f ´A pvq " A puq ´A pvq " A pv `eq ´A pvq . We can now solve eq. 4.22 on the coarse mesh:

A 2h `v2h `e2h ˘´A 2h `v2h ˘" r 2h , (4.23)
where r 2h " R 2h h `rh ˘" R 2h h `f h ´Ah `vh ˘˘is an approximation of the residual on the coarse mesh. Likewise, v 2h " R 2h h `vh ˘. By writing u 2h " v 2h `e2h , the error can be computed by solving A 2h `u2h ˘" A 2h `v2h ˘`r 2h , (4.24)

whose right-hand side is known. Once u 2h is computed by solving eq. 4.24, the new approximate solution on the fine mesh can be computed:

v h `Ph 2h ¨u2h ´v2h loooomoooon e 2h ' . (4.25)
This leads to the full approximation scheme (FAS) algorithm 5 [START_REF] Briggs | A Multigrid Tutorial[END_REF].

Algorithm 5 Two-grid full approximation scheme Pre-smoother: relax ν 0 times A h `uh ˘" f h with initial guess v h on the fine mesh

Ω h Restrict the residual r 2h " R 2h h `f h ´Ah `uh ˘˘and the current solution v 2h " R 2h h `uh
Solve A 2h `u2h ˘" r 2h `A2h `v2h ˘with initial guess v 2h on the coarse mesh Ω 2h Correct the approximation on the fine mesh: vh " v h `Ph 2h `u2h ´v2h

Post-smoother: relax ν 0 times A h `uh ˘" f h with initial guess vh on the fine mesh Ω h

The solution of the problem solved on the coarse mesh is the full approximation u 2h " v 2h `e2h , and not the error e 2h . Algorithm 5 is presented for only two meshes, but it can easily be extended to a higher number of meshes using a V-cycle, as done in section 4.2.2. If A is linear, one can show that algorithm 5 is equivalent to algorithm 4.

Algorithm 6 Resolution at coarse level using pseudo-time Choose ∆τ im m " 0 while the steady state is not reached do for each cell i 2h in the coarse grid do

τ ex i 2h " 0 K " 0 ´Ą u 2h ¯p0q i 2h " `u2h ˘m i 2h while τ ex i 2h ă ∆τ im do ´Ą u 2h ¯pK`1q i 2h " ´Ą u 2h ¯pKq i 2h `∆τ ex i 2h ``A 2h `v2h ˘˘i 2h `f 2h i 2h ´`R 2h h `Ah `vh ˘˘˘i 2h with ∆τ ex i 2h such that ´Ą u 2h ¯pK`1q i 2h is admissible K Ð K `1 τ ex i 2h Ð τ ex i 2h `∆τ ex i 2h end while Ą u 2h i 2h " ´Ą u 2h ¯pKq i 2h
end for Solve `u2h ˘m`1 `∆τ im A 2h ´`u 2h ˘m`1 ¯" Ą u 2h with nonlinear Jacobi method (algorithm 1) m Ð m `1 end while

The update of the solution with the correction on the fine grid

v h " v h `Ph 2h `u2h ´v2h ˘(4.29)
can also result in non-admissible states. Let us introduce another pseudo-time τ to compute the steady state of

dv h dτ `vh " v h `Ph 2h `u2h ´v2h ˘, (4.30)
instead of using eq. 4.29. We write ∆τ the interval between the current pseudo-time τ m and τ m`1 . We can now discretize eq. 4.30 with a time-implicit solver: time step than the one required by the explicit scheme to preserve the admissible states. Despite that, the solution obtained with the Newton-Raphson method is more diffusive than the one obtained with the explicit scheme. It reaches only 89% of the maximum value of radiative energy. On the contrary, solutions obtained with the Jacobi method (L max " 1) and with the geometric multigrid algorithm (L max " 2) reach 96% of the maximum of radiative energy, with a much larger time step (CFL " 2 000).

´vh ¯m`1 ´´v h ¯m ∆τ `´v h ¯m`1 " v h `Ph 2h `u2h ´v2h ˘. ( 4
Using the geometric multigrid method should reduce the computational cost. If we reach a low residual such as r h pκq r h p0q " 10 ´5, numerical experiments have shown that there is no gain in computational cost. Therefore, we set it to 10 -2 and we check some properties of the scheme in section 4.3.2. Table 4.2 shows the computational time needed to reach the steady state with different methods: the explicit HLL solver, the implicit HLL solver using the Newton-Raphson method to solve nonlinear solver and the geometric multigrid algorithm with different values for L max . The resolution used is now 257 ˆ257 cells. Using the explicit solver and the Newton-Raphson method, the time step is restricted by the CFL condition, whereas the steady state is reached with only one iteration using the geometric multigrid method. With L max " 1, we recover the Jacobi method and the steady state is reached three times faster than using the explicit solver. Increasing the value of L max leads again to a decrease in computational cost. The time needed to reach the steady state with L max " 1 is more than twice the time needed with L max " 4.

As shown by table 4.3, when L max increases, the computational time per V-cycle also increases, from 0.48 s with L max " 1 to 1.75 s with L max " 4. With L max ą 1, the computational time per V-cycle is more than three times the computational time per cycle with L max " 1. However, the total computational time (table 4.2) decreases when L max increases. As shown by fig. 4.11, the number of V-cycles to reach the same residual decreases as L max increases, from 1 267 with L max " 1 to 400 with L max " 4.

Table 4.4 shows the memory consumption with the different methods. As described in section 2.1.3, the Newton-Raphson method requires solving large sparse linear systems. For performance reasons, we store the matrix and its preconditioner. This leads 10

Scheme

Memory consumption (MB) Explicit 25 Newton-Raphson 127

L max " 1 63 L max " 2 74 L max " 3 77 L max " 4 78
Table 4.4 -Memory consumption for the explicit solver, the implicit solver using Newton-Raphson method and the geometric multigrid algorithm different values for L max . 0 0.5 1 0 0.5 1

F " ˆ1 0 ˙Ñ F " ˆ0 ´1˙Ó F " ˆ0 1 ˙Ò F " ˆ´1 0 ˙Ð Figure 4.
12 -Initial condition for 2D Riemann problem.

to a higher memory footprint than the other methods, 127 MB instead of 25 MB for the explicit solver. Furthermore, the choice of the preconditioner has an impact on the memory consumption of the method. This is discussed in section 5.2. We use a "matrix-free" approach for the Jacobi method, i.e., we do not store the operator A, but we access it by computing A pvq. However, using the geometric multigrid algorithm requires storing temporary values at coarse levels, hence the increase of memory consumption with L max , from 63 MB with L max " 1 to 78 MB with L max " 4. As L max increases, more values have to be stored, but each level is coarser than the previous one, therefore fewer variables per additional level are needed. For example, with 257 ˆ257 cells at the fine level and L max " 4, there is only 32 ˆ32 additional cells compared to L max " 3.

Performance results obtained here are obtained on small configurations: 257 ˆ257 cells for a two-dimensional problem. One can expect a better speed-up in the threedimensional case.

2D Riemann problem

Let us now consider the same two-dimensional Riemann problem as in section 3.2.1.1 and Blachère and Turpault [2016]. The domain r0, 1sˆr0, 1s is discretized with 257ˆ257 cells. The initial temperature is T 0 " T r " 1 000 K. The domain is cut into four states, in each of them the initial radiative flux is constant. It is set to p1 ´10 ´8q cE r F , with F given by fig. 4.12.

Using periodic boundary conditions, no energy should enter nor leave the box, therefore, the total radiative energy should be conserved at machine precision. Figure 4.13 shows the evolution of the relative error between the expected total radiative energy and the one actually computed in the box at each time step. Using both the explicit solver and the Jacobi method, with the same time step ∆t, the relative error oscillates around the value 10 -12 . For the Jacobi method, the residual r h pκq r h p0q is set to 5 ˆ10 -3 . Even for a quite high value of residual, the scheme is conservative.

We only checked that the Jacobi method (L max " 1) is conservative for performance reasons. As shown by fig. 4.14, with a large time step (CFL " 2 000), increasing the value of L max does not help to reduce the number of V-cycle to reach the same residual. When increasing the value of L max , the computational time per V-cycle increases, therefore the total computational time also increases, unlike the beam test. Because we are in the free-streaming regime, the propagation of the photons has to be followed. In sections 2.2.3 and 2.2.4, we did it by reducing the time step. Here, we still use a large time step, but the Jacobi method needs lots of iterations to converge. ). Using the Jacobi method, we only need one time step to reach the same final time as the explicit solver with more than 150 time steps. Because we use a larger time step, the solution obtained with the Jacobi method is more diffusive than the one obtained with the explicit solver. However, in this test, numerical diffusion increases as the residual decreases.

As discussed in section 3.2.1.1, the reduced flux stays close to 1 during the simulation. Overall, the methods presented in this chapter are stable, even for very stiff problems. But it should be used carefully in physical problems that require following the propagation of photons. In the previous chapters, we presented numerical schemes to solve the M 1 model of radiative transfer. However, studying three-dimensional astrophysical problems requires the use of High Performance Computing (HPC) tools, to reach large resolutions.

In this chapter, we focus on the scheme developed in chapter 2. Its implementation has been done in the code ARK-RT 1 , a fork of the code ARK developed in Padioleau et al. [2019]. The hydrodynamics and gravity part of the solver is similar to ARK and is solved explicitly with a well-balanced and all-regime solver. For the radiation part, the scheme is implicit and a Newton-Raphson method is used to solve the nonlinear system arising from the discretization of the M 1 model. At each iteration, a large sparse ill-conditioned linear system has to be solved. In order to do this, we use the library Trilinos [Trilinos Project Team], especially its second generation of packages.

We first present this library in section 5.1. In section 5.2, we provide details of some implementation features of ARK-RT. We also show some performance results and we discuss the impact of preconditioners and domain decomposition on a physical problem. Finally, we reach our conclusion in section 5.3.

Trilinos

The packages we used are Kokkos (section 5.1.1) for shared memory computation, Tpetra (section 5.1.2) for distributed vectors and matrices, Belos (section 5.1.3) for linear solvers, Ifpack2 (section 5.1.4) for classical preconditioners, and MueLu (section 5.1.5) for the algebraic multigrid (AMG) preconditioner. We now detail each of them. We only present the way we used Trilinos, it can be used for many others applications. This work does not intend to be a full review of the capabilities of Trilinos.

Kokkos

As new architectures have more and more cores, the distributed memory model is not enough 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 more and more heterogeneous, for example, multi-cores, many-cores, or accelerators such as GP-GPUs. Each architecture requires its own interface, such as OpenMP or C++11 threads for multi-cores and many-cores 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, so running the code on a different architecture will result in bad performance.

The package Kokkos [Edwards et al., 2014] tackles this issue. The user has a unique code that can be compiled with different shared memory models such as OpenMP or CUDA. We call these programming models backends.

As CUDA programming model, Kokkos relies on a host/device abstraction. The host dispatches the parallel work to the device. The device can be a GPU or a multicore CPU, for example.

The memory space abstraction defines where data are stored. It can be on the host (Kokkos::HostSpace) or on the device (e.g., Kokkos::CudaSpace for a GPU). Kokkos provides multidimensional arrays, whose storage is optimized according to the architecture. These arrays are called Kokkos::View.

The execution space specifies where a function is executed. The different possibilities are -on the host with a sequential execution (Kokkos::Serial), -on the host using POSIX threads (Kokkos::Threads), -on the host using OpenMP (Kokkos::OpenMP), -on the device, if a GPU is used (Kokkos::Cuda). Finally, the execution policy dictates how the function is executed. Kokkos provides three of them. Kokkos::parallel_for is a loop pattern. Each iteration executes the body of the function. Kokkos gives access to an index, but there is no guarantee about the order of the iterations. The second pattern is Kokkos::parallel_reduce that implements a reduction. Each iteration produces a result, and they are accumulated into a single value. The last pattern is Kokkos::parallel_scan, which is a cumulative reduction. For example, if the reduction operation is a sum, all partial sums are stored.

Tpetra

The package Tpetra [Tpetra Project Team] handles distributed linear algebra objects, such as sparse matrices and vectors. There are two levels of parallelism: the Message Passing Interface (MPI) for distributed memory parallelism and a shared memory programming model, handled through Kokkos (see section 5.1.1).

Tpetra relies on another package, Teuchos [Teuchos Project Team], that provides some tools common to all packages of Trilinos. The only class we mention here is Teuchos::RCP<T>. This class provides smart pointers, with an automatic garbage collector. The memory is managed automatically, the user does not deallocate it manually.

Before presenting the different objects we used, let us briefly present the template parameters used by Tpetra.

Template parameters

Trilinos is a general-purpose framework, the user has to specify several data types. Therefore, most of the objects use template parameters, such as local_ordinal_type, global_ordinal_type, scalar_type, and node_type.

The scalar_type is the type of values in the matrix or vectors. Common cases are float or double. The choice of float or double has an impact on the performances and the precision of the code.

Because Tpetra handles distributed objects, the rows of the matrix or vector are distributed across the MPI processes. Each of them is associated with a unique global index, whose type is global_ordinal_type. This type indicates the maximum size of the object, for example, the maximal number of rows of the matrix. If we use int as global_ordinal_type, the maximum number of rows in the matrix would 2 31 ´1, about two billion. In the three-dimensional case, there is five variables per cell (two scalars, the radiative energy and the gas temperature, and three components for the radiative flux). We are then restricted to 754 3 cells, so the total number of unknowns is less than 2 31 ´1. Therefore, Trilinos suggest using more than 64 bits integers. The default type is long long.

Inside a MPI process, the row of the matrix or vector is also associated with a local index. Its type is local_ordinal_type and is usually int.

The last type shared by Tpetra's objects is the node_type. It specifies the shared memory programming model, such as OpenMP or CUDA.

Most of the classes defined by Tpetra use template parameters. Let us consider a matrix class, with the four template parameters presented previously. A class with template parameters does not generate any code until all the parameters are specified. Therefore, part of the code of Trilinos is compiled only when the application code is compiled. Each file of the application code which uses a matrix has to include the header file where the matrix class is declared. So, the matrix class will be compiled several times, even if the same template parameters are used. This leads to a large compilation time when the application code is built against Trilinos.

To tackle this issue, explicit template instantiation (ETI) is used. When Trilinos is built, the template parameters are explicitly specified, and the classes are compiled. This increases the compilation time of Trilinos, but significantly reduces the compilation time of the application code. Listing 5.2 shows a way to simplify future code by defining template parameters used by Tpetra. Let us notice that these types can also be defined using default values when Trilinos is compiled with ETI (e.g. using scalar_type = Tpetra::MultiVector <>::scalar_type;).

Map

Using Tpetra vocabulary, a map describes how data are distributed across the MPI processes. Let us notice that, despite the name, it is not the same concept as C++ Standard Template Library (STL) std::map, although there are some common ideas.

Let us consider a matrix with N rows, distributed across p processes. For the sake of simplicity, let us assume that N is a multiple of p. Using Tpetra's default behavior, each process handles N p contiguous rows. This case is optimized by Tpetra, it should provide the best performances. From Tpetra documentation the constructor of a map must be called collectively, i.e., all processes have to call it. Listing 5.3 shows how to build a map.

Multivector

A multivector is a set of distributed vectors. Their distribution across the MPI processes is the same, i.e., they all have the same map. A multivector can contain only one vector. Having a set of vectors can be useful, for example, to solve a linear system with several right-hand sides, but the same matrix. Local to each MPI process, the multivector can be seen as a view, in the sense given by Kokkos (section 5.1.1).

Let us briefly recall the algorithm to solve radiation hydrodynamics (section 2.1.2):

1. update of the hydrodynamics quantities;

2. update of the radiative quantities and gas temperature using the implicit solver developed in section 2.1.1; 3. addition of source terms.

The first and the last steps do not require the use of whole library Trilinos, only Kokkos and MPI are used. The code ARK implements only explicit solvers, it is not coupled with Trilinos. Views used to store the hydrodynamics quantities and views coming from a multivector have to interact. For performance reasons, it should be done through a Kokkos::parallel_for, to avoid data transfers between the device and the host. Listing 5.4 shows how it can be done. In our application, the values of the matrix have to be updated at each iteration of the Newton-Raphson method, but its graph (the localization of the nonzero elements in the matrix) does not change. To obtain the best performances, Trilinos recommends building first the graph and then the matrix. The structure of the matrix is not allowed to change. The graph is distributed as the matrix, but values are not stored. Listings 5.5 to 5.7 show how to build the matrix obtained with the discretization of the Laplacian operator in the one-dimensional case. The matrix is

¨2 ´1 ´1 2 ´1 . . . . . . . . . ´1 2 ´1 ´1 2 ‹ ‹ ‹ ‹ ‹ ' .
(5.1)

In particular, listing 5.5 shows how to build this tridiagonal matrix with such a graph. Tpetra provides several methods to update the coefficients of the matrix. One of them uses only global indices, which is the way recommended by Tpetra. Listing 5.6 shows how this can be done. However, this function can only be called from 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 Tpetra. This will increase the computational cost. Another way is to use local indices. The package Kokkos Kernel, 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 (see listing 5.7). For performance reasons, we use local indices.

// Get t h e p a r t o f t h e matrix l o c a l t o t h e MPI

Belos

As discussed in section 2.1.3.1, the matrix of the linear system to be solved using the scheme developed in chapter 2 is large, sparse, and not symmetric. We solve the linear system with a biconjugate gradient stabilized (BiCGSTAB) method [Van der Vorst, 1992]. The package Belos [Belos Project Team] provides an implementation of this algorithm. Belos gives also access to other linear solvers we did not explore. Listing 5.8 shows how to create the linear problem associated with AX=B and the associated solver manager. It also shows how to solve the problem.

Ifpack2

In theory, the BiCGSTAB algorithm converges with a finite number of iterations. However, the convergence can be slow for physical applications. A similar phenomenon has already been observed in section 4.2.1. Instead of using a geometric multigrid method as in section 4.2, we use a preconditioning technique (see e.g., Saad 2003). The linear system to be solved is modified into another system, with the same solution, but easier to solve. Instead of solving Ax " b, we solve the right preconditioned system AK ´1Kx " b via solving AK ´1y " 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 ´1 is lower than A's one.

The package Ifpack2 [Ifpack2 Project Team] provides classical preconditioners. One of them, the Jacobi or diagonal preconditioner is given by

K ij " # A ij if i " j 0 otherwise. (5.2)
Ifpack2 also provides a relaxed incomplete LU factorization with level k fill (Ifpack2 ::RILUK). Instead of solving Ax " b, one can write A " LU, where L is a lower triangular matrix and U is an upper triangular matrix and then solve Ly " b followed by Ux " y. However, even if A is sparse, L and U are not sparse in general. This leads to a phenomenon called fill-in. Memory consumption can become a bottleneck. Therefore, other matrices L and U can be used, such that A « LU. The product LU is used as a preconditioner. For example, the nonzero elements of the new matrices L and U can be the same as A, leading to the method ILU(0). A more accurate method, but with a higher memory footprint is ILU(k), where the nonzero elements of L and U are the same as A k`1 . All tests in section 5.2 are performed with the ILU(0) method.

Ifpack2 implements another incomplete factorization (Ifpack2::ILUT) where the maximum number of entries to keep in each row of L and U is set according to the number of nonzero elements of A in the same row. Entries smaller than a threshold are dropped. See [START_REF] Saad | Ilut: A dual threshold incomplete lu factorization[END_REF] for more details.

Finally, Ifpack2 gives access to an addition Schwarz domain decomposition method (Ifpack2::AdditiveSchwarz). The problem is divided into smaller domains and the local results are added to each other. Each domain corresponds exactly to one MPI process, to avoid extra communications. More details can be found, for example in Prokopenko et al. [2016]. Listing 5.9 -Demonstration of the construction of an Ifpack2 preconditioner Listing 5.9 shows how to create a preconditioner using Ifpack2 and how it interacts with Belos linear problem.

MueLu

The last package we present here is MueLu [Muelu Project Team]. We used it as algebraic multigrid (AMG) 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 efficient to compute the high frequencies of the solution, but lack efficiency to compute its low frequencies. However, the computation is easier on a coarser grid with fewer unknowns. We recover some ideas already developed in section 4.2. The idea of the multigrid solver is to build a coarser grid, then solve the problem on this coarse grid and finally interpolate the solution on the fine grid. We can then 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 non squared matrices. From Saad [2003], here are the main steps of the method:

1. pre-smoothing: a few iterations of a simple method such as Jacobi or an incomplete factorization are performed, to get the value x;

2. the residual r " b ´Ax is projected over the coarse grid with the restriction operator R, to get the residual equation RAPy " Rr;

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 x " x `Py; 5. post-smoothing: a few iterations of a simple method are again performed to get the solution x.

The solution 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 (see section 4.2), whereas, with an algebraic multigrid solver, they are automatically generated, using data from the matrix. Algorithms provided by Ifpack2 (see section 5.1.4) are used as smoothers. Listing 5.10 shows how to create a preconditioner using MueLu and how it interacts with Belos linear problem. In the next section, we present some performance results obtained with the different packages of Trilinos. Especially, we explore the impact of the different preconditioners.

Application to radiative transfer and performance results

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 [START_REF] Kestener | Implementing high-resolution fluid dynamics solver in a performance portable way with kokkos[END_REF], 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, 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. See section 5.1.2 for the different ways to update the matrix and the right-hand side.

Performances

Thanks to Trilinos, we can use many preconditioners. Unfortunately, they do not behave 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 @ 2.6 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 is getting a piece of the whole domain of 1 500 2 cells, therefore the size of the system increases with the number of MPI tasks. The resolution is close to the one we are aiming for three-dimensional simulations. Figure 5.1 shows the mean number of iterations for the linear solver to converge as a function of the number of cells. For all preconditioners, the number of iterations remains constant, 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. Figure 5.2 shows the speed-up as a function of the number of MPI processes. The speed-up reaches a plateau of 80% to 90% of maximum performance, depending on the preconditioner. is constant, one can expect the number of iterations to also remain constant when the number of MPI processes increases. However, using the algebraic multigrid (orange curve) and the incomplete factorizations (green curve), 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 using two or four MPI processes. Furthermore, all tested preconditioners and the linear solver requires several communications per iteration, which likely become the main cost when the local resolution decreases.

Thanks to Kokkos, we can use exactly the same code on different architectures like Sandy Bridge processors and NVIDIA GP-GPUs (e.g., K80). Unfortunately, the memory required by the AMG preconditioner with a 1 500 2 simulation is larger than the memory available on a K80 GPU. For the next tests, we use a lower resolution of 1 000 2 cells. Table 5.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 needs several thousands of 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 CPU, the implicit solver is around 160 times faster than the explicit solver, whereas, on GPU, it is only 11 times faster. Figure 5.5 compares the computational time with different preconditioners, on both CPU and GPU. Except for the implicit solver using the AMG preconditioner, all solvers are faster on GPU than CPU, up to three times faster for the relaxation preconditioner. Part of the AMG algorithm probably remains sequential. On CPU, the AMG preconditioner is faster than the relaxation preconditioner. The other preconditioners are slower, up to a factor 8 between the relaxation and the additive Schwarz domain decomposition on GPU.

Figure 5.6 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 on the host, between 0.125 GB and 0.208 GB, 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).

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 CPU but is less efficient on GPU. Both incomplete factorizations and the additive Schwarz domain decomposition are slightly less efficient than the AMG preconditioner. Overall we have found the AMG preconditioner or relaxation method are a good compromise between stability and performances.

The performances we obtained thanks to Kokkos and Trilinos are encouraging for the study of astrophysical problems. The time step given by the hydrodynamics can be written as CF L ∆x c . Using a relaxation as a preconditioner, we need CF L ě 50 on CPU and CF L ě 100 on GPU to save computational time, whereas, using an incomplete factorization, we need CF L ě 250 on CPU and CF L ě 1 000 on GPU. We need a larger CFL number on GPU because the explicit solver is more efficient on GPU than CPU.

Effect of preconditioner and MPI domain decomposition on the expansion of H ii region

In the test case described in section 2.2.6, some numerical noise appears, as a consequence of the long timescales. Let us recall that a time-implicit scheme is used, with large time steps for the radiative transfer. At each time step, the Newton-Raphson method is used and, at each iteration of this algorithm, an ill-conditioned linear system is solved, using an iterative process. This results in the appearance of some numerical noise.

We have performed the simulation described previously, with different numbers of MPI processes and different preconditioners. The physical domain is either distributed over 4ˆ4 MPI processes (figs. 5.7a and 5.7c) or 2ˆ2 MPI processes (figs. 5.7b and 5.7d). We have also tried two preconditioners which allowed us to reach the final time with reasonable computational time: a standard ILU(k) factorization (figs. 5.7a and 5.7b) and an additive Schwarz domain decomposition (figs. 5.7c and 5.7d).

Figure 5.7 shows snapshots of the fraction of ionization and the velocity field at the final time t f " 10 10 s. The shape of the small structures produced by the numerical noise varies with the number of MPI processes and the preconditioner. Furthermore, the propagation of the ionization front creates some velocity that also depends on the number of MPI processes and preconditioners. However, the position of the ionization front is not affected by these parameters.

Discussion and conclusion

In this chapter, we explored some features of the library Trilinos for the resolution of linear systems. We also studied their performances on different architectures, such The goal of this thesis is to develop an asymptotic preserving implicit scheme for the M 1 model for radiative transfer. This scheme should also preserve the admissible states, E r ą 0 and ||Fr|| cEr ď 1. In chapter 1, we have presented the radiative transfer equation. This model is costly to solve numerically, we derive a two-moment model by averaging the specific intensity over the direction of propagation and the frequency of the photons. We have to specify a closure relation, to express the radiative pressure as a function of the radiative energy and the radiative flux. We have chosen to use the M 1 model for its good properties in both optically thick and thin media. We compute the radiative pressure by maximizing the radiative entropy. We also have investigated the behavior of this model in the diffusive limit, with long timescale and large opacity. Finally, we have presented the model used for radiation hydrodynamics, where the hydrodynamics is no longer frozen.

In chapter 2, we have developed a first numerical scheme. This scheme is asymptotic preserving and time-implicit. It is based on a standard HLL scheme, numerical fluxes are modified with an asymptotic correction that allows us to capture the correct behavior in the diffusive limit. We also propose a well-balanced modification of the source term, to capture the correct steady state with a discontinuity of opacity. The nonlinear system arising from this discretization is solved using the Newton-Raphson method. Using this numerical scheme, we study the stability of the ionization front in a massive pre-stellar dense core in the presence of convection. However, numerical experiments have shown that this scheme does not preserve the admissible states.

In chapter 3, we have derived a second numerical scheme, asymptotic preserving and entropic. Because this scheme is entropic, we have shown that it preserves the admissible states. It is based on a splitting of operator, first developed for compressible hydrodynamics. We have rewritten the M 1 model in a form similar to the Euler equations. The first scheme we presented in this chapter is an adaptation of the work of Chalons et al. [2016]. We have then modified this scheme to obtain an unsplit scheme, with stencil 1. For both schemes, we have used the numerical fluxes proposed by Buet and Despres [2008], developed for a Lagrange-remap method. These fluxes take into account source terms and allow the scheme to be asymptotic preserving. For the sake of simplicity, these two schemes are presented with a time-explicit integration.

In chapter 4, we have presented an algorithm to solve the nonlinear system arising from a time-implicit HLL solver. We used the Jacobi method developed by [START_REF] Pichard | Martin Mathematiques appliquees et calcul scientifique Bordeaux[END_REF]. It is an iterative process and each iteration results in an admissible state. Nevertheless, it should be used carefully in the free-streaming regime. Increasing the number of iterations can increase the numerical diffusion (see section 4.3.2). Furthermore, the Jacobi method can be slow to converge when the resolution increases. To tackle this issue, we have used a nonlinear geometric multigrid algorithm [START_REF] Briggs | A Multigrid Tutorial[END_REF]. To help the convergence on a fine grid, a coarser grid is built, and the problem is solved on the coarse grid. This process can be applied recursively to solve the problem on the coarse grid. The Jacobi method is used as a smoother and coarse grid solver. The quantities moved between the different grids are the residual of the equation and the error. These quantities are not admissible in most cases. Following [START_REF] Kifonidis | On multigrid solution of the implicit equations of hydrodynamics -experiments for the compressible euler equations in general coordinates[END_REF], we introduce a pseudo-time at the coarse level and we use it to preserve the admissible states.

Finally, in chapter 5, we have presented the library Trilinos, used for linear algebra in a high performance computing context. We used it to solve the linear systems involved in the scheme developed in chapter 2. The package Kokkos allows us to Title : Development and implementation of numerical schemes for radiation hydrodynamics Keywords : radiation hydrodynamics, asymptotic preserving scheme, all-regime scheme, geometric multigrid, high performance computing, linear algebra Abstract : The problem of radiative transfer describes the interaction between light and matte, therefore it appears in many astrophysical systems. This work aims at developing numerical methods for radiation hydrodynamics, when the radiation interacts with a fluid. We use the M1 model that is able to accurately capture the two main regimes arising in problems involving radiative transfer: the optically thin medium in which photons are freestreaming and the optically thick medium in which photons are constantly interacting and obey a diffusion equation in the asymptotic limit. We present several numerical schemes, with different properties: the correct behavior in the asymptotic limit, a time-implicit integration, and the preservation of the admissible states. Their implementation aims at being ready for highperformance computing, on exascale architectures. We apply one of them to the study of the propagation of an ionization front in a massive prestellar dense core.

4. 5

 5 Figure 1.1 -Nébuleuse d'Orion et différentes zones optiques. Source : https://esahubble. org/news/heic0601

1. 4

 4 Description of the work . . . . . . . . . . . . . . . . . . . . . 1.4.1 Radiative transfer . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 A first asymptotic preserving solver . . . . . . . . . . . . . . 1.4.3 An all-regime-like asymptotic preserving scheme . . . . . . . 1.4.4 Nonlinear Jacobi method and geometric multigrid . . . . . . 1.4.5 High performance computing and linear algebra . . . . . . . . 1.5 Publications and communications . . . . . . . . . . . . . . . 1.6 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 2 -

 12 Figure 1.2 -Orion nebula and different optical zones. Source: https://esahubble.org/ news/heic0601

Figure 1 .Figure 1 . 1 -

 111 Figure 1.1 summarizes the notations.Let us now derive the radiative transfer equation by computing the energy balance in an elementary volume dV with curvilinear length dl (see fig.1.2). We recall that the radiation can be absorbed by the surrounding gas. Photons could also be emitted

Figure 1 . 2 -

 12 Figure 1.2 -Elementary volume in which the energy balance is computed.

Figure 1 . 3 -

 13 Figure 1.3 -Spherical coordinates.

  cos θq 4 dθ " 0.

r 2 r 2

 22 pm log m ´pm `1q log pm `1qq dr, pM α log M α ´pM α `1q log pM α `1qq dr, (1.73)

Figure 2 .

 2 Figure2.1 -Marshak wave simulation. This figure shows a snapshot of the gas temperature at time t f " 2 ˆ10 ´4 s, with and without the asymptotic correction and the reference solution. Spatial resolution is 400 points and the opacity is σ " 10 000 cm ´1.

x

  Figure2.2 -Simulation of a steady state with a jump of opacity. The opacity is piecewise constant, a jump is located at x " 0.5 cm (gray line). This figure shows a snapshot of the radiative flux at time t f " 10 ´3 s.

Figure 2 . 3 -Figure 2 . 4 -

 2324 Figure 2.3 -Beam simulation. The figure shows the radiative energy. The eigenvalues are fixed to ˘c (left panel) or calculated with eq. 2.3 (right panel).

Figure 2

 2 Figure 2.5 -Shadow simulation. Snapshots of the radiative temperature at time t f " 10 ´10 s with different closure relations: P 1 model (upper panel), M 1 model with fixed eigenvalues (middle panel) and M 1 model with computed eigenvalues (lower panel).

Figure 2 .Figure 2

 22 Figure 2.6 -Shadow simulation. The figure shows the radial profiles of the radiative temperature at time t f " 10 ´10 s with different closure relations: P 1 model, M 1 model with fixed eigenvalues and M 1 model with computed eigenvalues.

Figure 2

 2 Figure 2.7 -Subcritical shock simulation. The figure shows snapshots of gas temperature, radiative temperature, and reduced flux at different times: 1.7 ˆ10 4 s, 2.8 ˆ10 4 s, and 3.8 ˆ10 4 s.

Figure 2

 2 Figure 2.8 -Supercritical shock simulation. The figure shows snapshots of gas temperature, radiative temperature, and reduced flux at different times: 4 ˆ10 3 s, 7.5 ˆ10 3 s, and 1.3 ˆ10 4 s.

  Time (s) Position of the ionization front (AU)Without the initial velocity perturbation With the initial velocity perturbation

Figure 2 .

 2 Figure 2.9 -Evolution of the position of the ionization front as a function of time, with and without the initial velocity perturbation.

Figure 2 .

 2 9 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[START_REF] Tremblin | Three-dimensional simulations of globule and pillar formation around hii regions: turbulence and shock curvature[END_REF].

)Figure 2 . 10 -

 210 Figure 2.10 -Snapshots of the fraction of ionization and the velocity field at the final time t f " 10 10 s without the initial velocity perturbation (left panel) and with it (right panel).

Figure 2 .

 2 Figure2.10 shows the ionization front at the final time t f " 10 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. The fraction of ionization, which is always between 0 and 1, reaches values between 10 -12 and 10 -6 . The effect of the preconditioner and the MPI domain decomposition is discussed in section 5.2.3. However, the numerical noise does not affect the position of the ionization front.

Contents 3 . 1

 31 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Rewriting the M 1 model . . . . . . . . . . . . . . . . . . . . . 3.1.2 Diffusion-transport splitting . . . . . . . . . . . . . . . . . . . 3.1.3 Diffusion step . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Transport step . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.5 Scheme of stencil 1 . . . . . . . . . . . . . . . . . . . . . . . . 3.1.6 Source terms . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Numerical results without source terms . . . . . . . . . . . . 3.2.2 Numerical results with source terms . . . . . . . . . . . . . . 3.3 Discussion and conclusion . . . . . . . . . . . . . . . . . . . 3.3.1 Well-balanced discretization of the source terms . . . . . . . . 3.3.2 Time implicit integration . . . . . . . . . . . . . . . . . . . . 3.4 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.A Equation of evolution for pressure and velocity in the diffusion step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.B Eigenstructure of the diffusion system . . . . . . . . . . . . 3.C Derivation of β ˚and q ˚. . . . . . . . . . . . . . . . . . . . . 3.C.1 Relaxed set of equations . . . . . . . . . . . . . . . . . . . . . 3.C.2 Eigenstructure . . . . . . . . . . . . . . . . . . . . . . . . . . 3.C.3 Riemann invariants for a hyperbolic system . . . . . . . . . . 3.C.4 Riemann invariants for the relaxed system . . . . . . . . . . . 3.D Admissible states . . . . . . . . . . . . . . . . . . . . . . . . . 3.D.1 Diffusion step . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.D.2 Transport step . . . . . . . . . . . . . . . . . . . . . . . . . . 3.D.3 Split scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.D.4 Unsplit scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 3.E Derivation of the unsplit scheme . . . . . . . . . . . . . . . 3.E.1 Relaxed set of equations . . . . . . . . . . . . . . . . . . . . . 3.E.2 Eigenstructure . . . . . . . . . . . . . . . . . . . . . . . . . . 3.E.3 Jump conditions . . . . . . . . . . . . . . . . . . . . . . . . .

  3.1. Figure 3.2 shows the radiative energy (figs. 3.2a, 3.2c and 3.2e) and the reduced flux (figs. 3.2b, 3.2d and 3.2f) at the final time 2 ˆ10 -11 s. Solutions are obtained with a standard HLL solver (figs. 3.2a and 3.2b), the split scheme (figs. 3.2c and 3.2d) and the unsplit scheme (figs. 3.2e and 3.2f). As shown by figs. 3.2b, 3.2d and 3.2f, the reduced flux stays close to 1 during the simulation. Even for a very stiff problem, our all-regime schemes can preserve the admissible states E r ą 0 and f ď 1, as shown in section 3.D.

Figure 3 . 2 -Figure 3 . 3 -Figure 3 . 4 -

 323334 Figure 3.2 -Snapshots of radiative energy and reduced flux at final time t f " 2 ˆ10 ´11 s with the HLL solver, the split scheme and the unsplit scheme.

Figure 3 . 9 -

 39 Figure 3.9 -Solution of the Riemann problem for b, Π, and a.

  .84) Value of D and b rrs , rµs , rϕs Wave D " 0, b " 0

Figure 3 . 10 -

 310 Figure 3.10 -Solution of the Riemann problem for b and Π.

Figure 3 . 11 -

 311 Figure 3.11 -Solution of the Riemann problem for b, Π, r, µ, ϕ with D " b L .

Figure 3 .

 3 Figure 3.14 -Solution of the Riemann problem for b, Π, r, µ, ϕ with D " b R .

Figure 4 .

 4 Figure 4.6 -V-cycle with two to four levels. Circles represent the different meshes, black lines represent the restriction and prolongation operators.

Figure 4 .

 4 Figure 4.7 -Evolution of the residual as a function of the number of V-cycle (algorithm 2), with different values of L max .

Figure 4 . 8 - 10 ´1

 4810 Figure 4.8 -Restriction of a smooth function.

Figure 4 . 9 -

 49 Figure 4.9 -Evolution of the residual as a function of the number of V-cycle (algorithm 4), with different values of L max .

L

  max " 1 L max " 2 L max " 3 L max " 4

Figure 4 . 11 -

 411 Figure 4.11 -Evolution of the residual as a function of the number of V-cycle, with different values for L max to reach the steady state for the beam problem.

Figure 4 .

 4 Figure 4.13 -Evolution of the relative error on radiative energy as a function of time. With periodic boundary conditions, the radiative energy is conserved at machine precision.

L

  max " 1 L max " 2 L max " 3 L max " 4

Figure 4 .

 4 Figure 4.14 -Evolution of the residual as a function of the number of V-cycle, with different values for L max for the two-dimensional Riemann problem.

´2Figure 4 .

 4 Figure 4.15 -Snapshots of radiative energy at final time t f " 10 ´11 s with the explicit solver and the Jacobi method with different residuals.

Figure 4 .

 4 Figure 4.15 shows snapshots of radiative energy at the final time 10 -11 s with the explicit solver (fig. 4.15a) and Jacobi method with different residuals (figs. 4.15b and 4.15c). Using the Jacobi method, we only need one time step to reach the same final time as the explicit solver with more than 150 time steps. Because we use a larger time step, the solution obtained with the Jacobi method is more diffusive than the one obtained with the explicit solver. However, in this test, numerical diffusion increases as the residual decreases.

  Sections 5.1 and 5.2 are the adaptation of an article published in Astronomy&Astrophysics, see[START_REF] Bibliography | A high-performance and portable asymptotic preserving radiation hydrodynamics code with the m1 model[END_REF]

  i n t main ( i n t argc , c h a r * argv [ ] ) { // I n i t i a l i z e Kokkos Kokkos : : i n i t i a l i z e ( argc , argv ) ; // D e c l a r e a view with 3 * 15 e l e m e n t s Kokkos : : View <d o u b l e * [3] > view ( " view_name " , 1 5 ) ; // I n i t i a l i z e t h e view through a p a r a l l e l l o o p Kokkos : : p a r a l l e l _ f o r ( " loop_name " , 1 5 , [ = ] ( c o n s t i n t i ) { view ( i , 0 ) = i ; view ( i , 1 ) = i * i ; view ( i , 2 ) = i * i * i ; } ) ; // Shut down Kokkos Kokkos : : f i n a l i z e ( ) ; } Listing 5.1 -Example of code using Kokkos: declaration and initialization of a Kokkos::View Listing 5.1 is an example of code showing how to declare a Kokkos::View and how to initialize it with a parallel loop.

  u s i n g s c a l a r _ t y p e = d o u b l e ; u s i n g g l o b a l _ o r d i n a l _ t y p e = l o n g l o n g ; u s i n g l o c a l _ o r d i n a l _ t y p e = i n t ; u s i n g node_type = Kokkos : : Compat : : KokkosCudaWrapperNode ; Listing 5.2 -Definition of template parameters used by Tpetra.

  u s i n g map_type = Tpetra : : Map<l o c a l _ o r d i n a l _ t y p e , g l o b a l _ o r d i n a l _ t y p e , node_type >; // C r e a t e a communicator , such a s MPI_COMM_WORLD Teuchos : : RCP<c o n s t Teuchos : : Comm<i n t >> comm = Tpetra : : getDefaultComm ( ) ; // L o c a l number o f e l e m e n t s c o n s t l o c a l _ o r d i n a l _ t y p e numLocalElements = . . . ; // T o t a l number o f e l e m e n t s // L o c a l number o f e l e m e n t s * number o f MPI p r o c e s s e s c o n s t g l o b a l _ o r d i n a l _ t y p e numGlobalElements = numLocalElements * comm´>g e t S i z e ( ) ; // Gl ob a l i n d i c e s s t a r t a t 0 c o n s t g l o b a l _ o r d i n a l _ t y p e indexBase = 0 ; // C r e a t e t h e map Teuchos : : RCP<c o n s t map_type> map = Teuchos : : r c p ( new map_type ( numGlobalElements , indexBase , comm) ) ; Listing 5.3 -Demonstration of the construction of a Tpetra::Map

u s i n

  g vec_type = Tpetra : : MultiVector <s c a l a r _ t y p e , l o c a l _ o r d i n a l _ t y p e , g l o b a l _ o r d i n a l _ t y p e , node_type >; // View c o n t a i n i n g hydrodynamics q u a n t i t i e s Kokkos : : View<d o u b l e * , Kokkos : : D e f a u l t E x e c u t i o n S p a c e > hydroData = . . . ; // C r e a t e m u l t i v e c t o r s c o n t a i n i n g one v e c t o r Teuchos : : RCP<vec_type> X = Teuchos : : r c p ( new vec_type (map , 1 ) ) ; Teuchos : : RCP<vec_type> B = Teuchos : : r c p ( new vec_type (map , 1 ) ) ; // S e t i n i t i a l g u e s s t o 0 X´>p u t S c a l a r ( 0 . ) ; // I n i t i a l i z e RHS G´>f i l l C o m p l e t e ( ) ; // Use t h e graph t o b u i l d t h e matrix Teuchos : : RCP<matrix_type> A = Teuchos : : r c p ( new matrix_type (G) ) ;// T e l l t h e matrix we a r e done c h a n g i n g i t s s t r u c t u r e A´>f i l l C o m p l e t e ( ) ; Listing 5.5 -Demonstration of the construction of a tridiagonal matrix with a static graph

/

  / T e l l t h e matrix we a r e about t o change i t s v a l u e s A´>r e s u m e F i l l ( ) ; // Update t h e c o e f f i c i e n t s f o r ( l o c a l _ o r d i n a l _ t y p e l o c a l I n d e x = 0 ; l o c a l I n d e x < numLocalElements ; l o c a l I n d e x ++) { // Ask t h e map t h e g l o b a l i n d e x c o r r e s p o n d i n g t o t h e l o c a l i n d e x c o n s t g l o b a l _ o r d i n a l _ t y p e g l o b a l I n d e x = map´>g e t G l o b a l E l e m e n t ( l o c a l I n d e x ) ; // Use g l o b a l i n d e x t o a c c e s s t h e e l e m e n t s o f t h e matrix i f ( g l o b a l I n d e x == 0 ) A´>r e p l a c e G l o b a l V a l u e s ( g l o b a l I n d e x , Teuchos : : t u p l e ( g l o b a l I n d e x , g l o b a l I n d e x +1) , Teuchos : : t u p l e ( 2 . , ´1.) ) ; e l s e i f ( g l o b a l I n d e x == numGlobalElements ´1) A´>r e p l a c e G l o b a l V a l u e s ( g l o b a l I n d e x , Teuchos : : t u p l e ( g l o b a l I n d e x ´1, g l o b a l I n d e x ) , Teuchos : : t u p l e ( ´1. , 2 . ) ) ; e l s e A´>r e p l a c e G l o b a l V a l u e s ( g l o b a l I n d e x , Teuchos : : t u p l e ( g l o b a l I n d e x ´1, g l o b a l I n d e x , g l o b a l I n d e x +1) , Teuchos : : t u p l e ( ´1. , 2 . , ´1.) ) ; } // T e l l t h e matrix we a r e done c h a n g i n g i t s v a l u e s A´>f i l l C o m p l e t e ( ) ;Listing 5.6 -Updating the coefficients of the matrix from the host (non-parallel update)

  p r o c e s s auto l o c a l M a t r i x = A´>g e t L o c a l M a t r i x ( ) ; // Get t h e l o c a l map , which can be used from Kokkos k e r n e l auto localMap = map´>getLocalMap ( ) ; // Use a p a r e l l e l l o o p t o update t h e c o e f f i c i e n t s o f t h e matrix Kokkos : : p a r a l l e l _ f o r ( " Update matrix c o e f f i c i e n t s " , numLocalElements , [ = ] ( c o n s t i n t l o c a l I n d e x ) { // Ask t h e map t h e g l o b a l i n d e x c o r r e s p o n d i n g t o t h e l o c a l i n d e x c o n s t g l o b a l _ o r d i n a l _ t y p e g l o b a l I n d e x = localMap . g e t G l o b a l E l e m e n t ( l o c a l I n d e x ) ; i f ( g l o b a l I n d e x == 0 ) { // T e l l Tpetra which column ( s ) t o update l o c a l _ o r d i n a l _ t y p e c o l [ 2 ] ; c o l [ 0 ] = l o c a l I n d e x ; c o l [ 1 ] = l o c a l I n d e x +1; // T e l l Tpetra which v a l u e ( s ) t o s e t s c a l a r _ t y p e v a l [ 2 ] ; v a l [ 0 ] = 2 . ; v a l [ 1 ] = ´1.; // Update t h e matrix l o c a l M a t r i x . r e p l a c e V a l u e s ( l o c a l I n d e x , c o l , 2 , v a l ) ; } e l s e i f ( g l o b a l I n d e x == numGlobalElements ´1) { // T e l l Tpetra which column ( s ) t o update l o c a l _ o r d i n a l _ t y p e c o l [ 2 ] ; c o l [ 0 ] = l o c a l I n d e x ´1; c o l [ 1 ] = l o c a l I n d e x ; // T e l l Tpetra which v a l u e ( s ) t o s e t s c a l a r _ t y p e v a l [ 2 ] ; v a l [ 0 ] = ´1.; v a l [ 1 ] = 2 . ; // Update t h e matrix l o c a l M a t r i x . r e p l a c e V a l u e s ( l o c a l I n d e x , c o l , 2 , v a l ) ; } e l s e { // T e l l Tpetra which column ( s ) t o update l o c a l _ o r d i n a l _ t y p e c o l [ 3 ] ; c o l [ 0 ] = l o c a l I n d e x ´1; c o l [ 1 ] = l o c a l I n d e x ; c o l [ 2 ] = l o c a l I n d e x +1; // T e l l Tpetra which v a l u e ( s ) t o s e t s c a l a r _ t y p e v a l [ 3 ] ; v a l [ 0 ] = ´1.; v a l [ 1 ] = 2 . ; v a l [ 2 ] = ´1.; // Update t h e matrix l o c a l M a t r i x . r e p l a c e V a l u e s ( l o c a l I n d e x , c o l , 3 , v a l ) ; } } ) ; Listing 5.7 -Updating the coefficients of the matrix from the device (parallel update)

  u s i n g op_type = Tpetra : : Operator<s c a l a r _ t y p e , l o c a l _ o r d i n a l _ t y p e , g l o b a l _ o r d i n a l _ t y p e , node_type >; u s i n g problem_type = B e l o s : : LinearProblem<s c a l a r _ t y p e , vec_type , op_type >; // Create , c o n f i g u r e , and r e t u r n t h e s p e c i f i e d s o l v e r B e l o s : : S o l v e r F a c t o r y <s c a l a r _ t y p e , vec_type , op_type> f a c t o r y ; Teuchos : : RCP<B e l o s : : SolverManager<s c a l a r _ t y p e , vec_type , op_type>> s o l v e r = f a c t o r y . c r e a t e ( "BiCGSTAB" ) ; // L i n e a r system AX = B t o be s o l v e d , and i t s a s s o c i a t e d i n f o r m a t i o n s Teuchos : : RCP<problem_type> problem = Teuchos : : r c p ( new problem_type (A, X, B) ) ; // I f needed , do something with a p r e c o n d i t i o n n e r // Prepare t h e problem t o s o l v e t h e l i n e a r system problem´>setProblem ( ) ; // T e l l t h e s o l v e r what problem you want t o s o l v e s o l v e r ´>setProblem ( problem ) ; // S o l v e t h e l i n e a r system // B e l o s : : ReturnType can t a k e t h e v a l u e s B e l o s : : Converged o r B e l o s : : Unconverged B e l o s : : ReturnType r e s u l t = s o l v e r ´>s o l v e ( ) ; // Get t h e number o f i t e r a t i o n s needed t o s o l v e t h e system c o n s t i n t numIters = s o l v e r ´>getNumIters ( ) ; // Get t h e a c h i e v e d t o l e r a n c e r e a c h e d by c a l l i n g s o l v e ( ) c o n s t s c a l a r _ t y p e a c h i e v e d T o l = s o l v e r ´>a c h i e v e d T o l ( ) ; // Do something t o r e a d t h e s o l u t i o n X, t o change t h e v a l u e s o f t h e matrix A and t h e r i g h t ´hand s i d e B // T e l l t h e s o l v e r manager t o p r e p a r e t h e s o l v e r f o r t h e next s o l v e s o l v e r ´>r e s e t ( B e l o s : : Problem ) ; // S o l v e a g a i n t h e l i n e a r system r e s u l t = s o l v e r ´>s o l v e ( ) ; Listing 5.8 -Demonstration of the construction of a Belos solver

  u s i n g i f p a c k 2 _ p r e c _ t y p e = I f p a c k 2 : : P r e c o n d i t i o n e r <s c a l a r _ t y p e , l o c a l _ o r d i n a l _ t y p e , g l o b a l _ o r d i n a l _ t y p e , node_type >; // C r e a t e an i n s t a n c e o f t h e p r e c o n d i t i o n e r , h e r e i t i s a d i a g o n a l p r e c o n d i t i o n e r I f p a c k 2 : : Factory f a c t o r y ; Teuchos : : RCP<ifpack2_prec_type > K = f a c t o r y . c r e a t e <matrix_type >( "DIAGONAL" , A) ; // S e t up t h e graph s t r u c t u r e o f t h e p r e c o n d i t i o n e r K´> i n i t i a l i z e ( ) ; // S e t up t h e n u m e r i c a l v a l u e s i n t h e p r e c o n d i t i o n e r K´>compute ( ) ; // Use K a s r i g h t p r e c o n d i t i o n e r f o r B e l o s l i n e a r problem problem´>s e t R i g h t P r e c (K) ; // S o l v e t h e l i n e a r system s o l v e r ´>s o l v e ( ) ; // Do something t o change t h e v a l u e s o f t h e matrix A and t h e r i g h t ´hand s i d e B // Compute a g a i n n u m e r i c a l v a l u e s i n t h e p r e c o n d i t i o n e r K´>compute ( ) ; // T e l l t h e s o l v e r manager t o p r e p a r e t h e s o l v e r f o r a next s o l v e s o l v e r ´>r e s e t ( B e l o s : : Problem ) ; // S o l v e a g a i n t h e l i n e a r system s o l v e r ´>s o l v e ( ) ;

u s i n g

  muelu_prec_type = MueLu : : TpetraOperator<s c a l a r _ t y p e , l o c a l _ o r d i n a l _ t y p e , g l o b a l _ o r d i n a l _ t y p e , node_type >; // C r e a t e a MueLu p r e c o n d i t i o n e r t h a t can be used by Tpetra Teuchos : : RCP<muelu_prec_type> K = MueLu : : C r e a t e T p e t r a P r e c o n d i t i o n e r <s c a l a r _ t y p e , l o c a l _ o r d i n a l _ t y p e , g l o b a l _ o r d i n a l _ t y p e , node_type >(A) ; // Use K a s r i g h t p r e c o n d i t i o n e r f o r B e l o s l i n e a r problem problem´>s e t R i g h t P r e c (K) ; // S o l v e t h e l i n e a r system s o l v e r ´>s o l v e ( ) ; // Do something t o change t h e v a l u e s o f t h e matrix A and t h e r i g h t ´hand s i d e B // Update t h e p r e c o n d i t i o n e r MueLu : : R e u s e T p e t r a P r e c o n d i t i o n e r (A, * K) ; // T e l l t h e s o l v e r manager t o p r e p a r e t h e s o l v e r f o r a next s o l v e s o l v e r ´>r e s e t ( B e l o s : : Problem ) ; // S o l v e a g a i n t h e l i n e a r system 14 s o l v e r ´>s o l v e ( ) ; Listing 5.10 -Demonstration of the construction of a MueLu preconditioner

  Figures 5.3 and 5.4 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. It is set to 2 048 2 cells. Because the global resolution 4

Figure 5

 5 Figure 5.1 -Number of iterations in a weak scaling test. Each MPI process treats 1 500 2 cells. We have tested different preconditioners: Jacobi with damping (Relaxation), algebraic multigrid (AMG), standard ILU(k) factorization (RILUK), variant of the standard ILU factorization (ILUT) and additive Schwarz domain decomposition (Schwarz).

Figure 5 . 2 -

 52 Figure 5.2 -Speed-up as a function of the number of MPI processes in a weak scaling test. Each MPI process treats 1 500 2 cells. We have tested different preconditioners: Jacobi with damping (Relaxation), algebraic multigrid (AMG), standard ILU(k) factorization (RILUK), variant of the standard ILU factorization (ILUT) and additive Schwarz domain decomposition (Schwarz).

Figure 5 . 3 -

 53 Figure 5.3 -Number of iterations in a strong scaling test. The global resolution is 2 048 2 cells. We have tested different preconditioners: Jacobi with damping (Relaxation), algebraic multigrid (AMG), standard ILU(k) factorization (RILUK), variant of the standard ILU factorization (ILUT) and additive Schwarz domain decomposition (Schwarz).

Figure 5 . 4 -

 54 Figure 5.4 -Speed-up as a function of the number of MPI processes in a strong scaling test. The global resolution is 2 048 2 cells. We have tested different preconditioners: Jacobi with damping (Relaxation), algebraic multigrid (AMG), standard ILU(k) factorization (RILUK), variant of the standard ILU factorization (ILUT) and additive Schwarz domain decomposition (Schwarz).

Figure 5 . 5 -

 55 Figure 5.5 -Computational time for the implicit solver with different preconditioners (Jacobi with damping (Relaxation), algebraic multigrid (AMG), ILU(k) factorization (RILUK), slightly modified variant of ILU factorization (ILUT) and additive Schwarz domain decomposition (Schwarz)) on different architectures (Sandy Bridge CPU and K80 NVIDIA GPU). The resolution is 1000 2 cells.

RelaxationFigure 5

 5 Figure 5.6 -Memory consumption for the implicit solver with different preconditioners (Jacobi with damping (Relaxation), algebraic multigrid (AMG), ILU(k) factorization (RILUK), slightly modified variant of ILU factorization (ILUT) and additive Schwarz domain decomposition (Schwarz)) on different architectures (Sandy Bridge CPU and K80 NVIDIA GPU). The resolution is 1000 2 cells.

4 ˆ4

 4 MPI processes, ILU(k) factorization.

2 ˆ2 2 ˆ2

 22 MPI processes, ILU(k) factorization. MPI processes, Schwarz domain decomposition.

Figure 5

 5 Figure 5.7 -Snapshots of the fraction of ionization and the velocity field at the final time t f " 10 10 s without the initial velocity perturbation (left panel) and with it (right panel). The physical domain is distributed across different numbers of MPI processes and different preconditioners have been used. Figure 5.7a is the same figure as fig. 2.10.

  Maison du doctorat de l'Université Paris-Saclay 2 ème étage aile ouest, Ecole normale supérieure Paris-Saclay 4 avenue des Sciences, 91190 Gif sur Yvette, France Titre : Développement et implémentation de schémas numériques pour l'hydrodynamique radiative Mots clés : hydrodynamique radiative, schémas préservant l'asymptotique, schéma tout-régime, multigrille géométrique, calcul haute performance, algèbre linéaire Résumé : Le transfert radiatif décrit les interactions entre la lumière et la matière, il apparaît donc dans de nombreux systèmes astrophysiques. Ce travail vise à développer des méthodes numériques pour l'hydrodynamique radiative, lorsque le rayonnement interagit avec un fluide. Nous utilisons le modèle M1 qui est capable de capturer précisément les deux régimes rencontrés dans les problèmes impliquant du transfert radiatif : le milieu optiquement fin dans lequel les photons circulent librement et le milieu optiquement épais dans lequel les photons interagissent constamment et obéissent à une équation de diffusion dans la limite asymptotique. Nous présentons plusieurs schémas numériques, avec différentes propriétés : le bon comportement dans la limite de diffusion, une intégration implicite en temps et la préservation des états admissibles. Leur implémentation vise à être prête pour le calcul haute performance, sur des architectures exascales. Nous appliquons l'un d'eux à l'étude de la propagation d'un front d'ionisation dans un coeur dense préstellaire massif.

C Diffusive limit for radiation hydrodynamics . . . . . . . . .

  Gary E. Thomas and Knut Stamnes. Radiative Transfer in the Atmosphere and Ocean. Van Henson, and Steve McCormick. A Multigrid Tutorial, 2nd Edition. 01 2000. ISBN 978-0-89871-462-3. 12, 13 Christophe Chalons, Pierre Kestener, Samuel Kokh, and Maxime Stauffert. A large time-step and well-balanced Lagrange-Projection type scheme for the shallow-water equations. July 2016. URL https://hal.archives-ouvertes.fr/hal-01297043. Mathematical modelling for dose depositon in photontherapy. PhD thesis, 2016. URL http://www.theses.fr/2016BORD0177. Thèse de doctorat dirigée par Dubroca, Bruno et Frank, Martin Mathematiques appliquees et calcul scientifique Bordeaux 2016. 13 Sections 1.4, 1.5 and 1.C are the adaptation of an article published in Astron-omy&Astrophysics, see Bloch et al. 2021.

	working paper or preprint. 13 Christophe Buet and Bruno Despres. A gas dynamics scheme for a two moments model of radiative transfer. 11 2008. 13 Teddy Pichard. Chapter 1 Radiative transfer

2002. 10 Rolf Kuiper, Hubert Klahr, Henrik Beuther, and Thomas Henning. CIRCUM-VENTING THE RADIATION PRESSURE BARRIER IN THE FORMATION OF MASSIVE STARS VIA DISK ACCRETION. The Astrophysical Journal, 722(2): 1556-1576, oct 2010. doi: 10.1088/0004-637x/722/2/1556. URL https://doi.org/ 10.1088/0004-637x/722/2/1556. 10 William Briggs, K. Kifonidis and E. Müller. On multigrid solution of the implicit equations of hydrodynamics -experiments for the compressible euler equations in general coordinates. A&A, 544:A47, 2012. doi: 10.1051/0004-6361/201116979. URL https: //doi.org/10.1051/0004-6361/201116979. 14 Contents 1.1 Radiative transfer equation . . . . . . . . . . . . . . . . . . 1.2 Moment model . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Closure relation . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Diffusive limit . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Radiation hydrodynamics . . . . . . . . . . . . . . . . . . . 1.6 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A Computation of source terms . . . . . . . . . . . . . . . . . 1.B Maximization of radiative entropy . . . . . . . . . . . . . . 1.B.1 Radiative entropy . . . . . . . . . . . . . . . . . . . . . . . . 1.B.2 Maximization of entropy . . . . . . . . . . . . . . . . . . . . . 1.B.3 M 1 closure relation . . . . . . . . . . . . . . . . . . . . . . . . 1.B.4 Back to the radiative entropy . . . . . . . . . . . . . . . . . . 1.

Chapter 2 A first asymptotic preserving solver Contents 2.1 Numerical scheme and algorithm

. . . . . . . . . . . . . . . 2.1.1 Radiation transport in a static fluid . . . . . . . . . . . . . . 2.1.2 Coupling to hydrodynamics . . . . . . . . . . . . . . . . . . . 2.1.3 Algorithm for nonlinear implicit solver . . . . . . . . . . . . . 2.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Marshak wave . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Steady state with a jump of opacity . . . . . . . . . . . . . . 2.2.3 Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.4 Shadow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.5 Radiative shocks . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.6 Expansion of H ii region . . . . . . . . . . . . . . . . . . . . . 2.3 Discussion and conclusion . . . . . . . . . . . . . . . . . . . 2.3.1 Well-balanced discretization of the source term . . . . . . . . 2.3.2 Asymptotic limit for radiation hydrodynamics . . . . . . . . .

2.4 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.A Von Neumann stability analysis for the well-balanced mod- ification of the source term . . . . . . . . . . . . . . . . . . 2.B Numerical scheme in the diffusive limit . . . . . . . . . . .

  .26) and we have then F n`1 pkq " Apkq F n pkq. By induction, we have F n pkq " pApkqq n F 0 pkq. The coefficient F n pkq remains bounded if and only if |Apkq| ď 1. In this case, for all k P Z, F n`1 pkq ď F n pkq . Therefore, ||F n`1 || 2 ď ||F n || 2 ď ||F 0 || 2 and the scheme is unconditionally stable. We now have to prove that |Apkq| ď 1:

  1 -Initial condition for 2D Riemann problem.

	Scheme	Computational time (s)
	HLL with computed eigenvalues	44
	HLL with fixed eigenvalues	26
	Split	79
	Unsplit	36

Table 3

 3

.1 -Computational time with the HLL solver, with computed and fixed eigenvalues (eq. 2.

  2016].

	François Bouchut, Christian Klingenberg, and Knut Waagan. A multiwave approximate
	riemann solver for ideal mhd based on relaxation ii: numerical implementation with
	3 and 5 waves. Numerische Mathematik, 115(4):647-679, 2010. ISSN 0029-599X. 77
	François Bouchut. Nonlinear Stability of Finite Volume Methods for Hyperbolic Con-
	servation Laws and Well-Balanced Schemes for Sources. Frontiers in Mathematics.
	1st ed. 2004.. edition, 2004. ISBN 3-7643-6665-6. 79
	Eleuterio Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics, pages
	163-212. 03 2009. ISBN 978-3-540-25202-3. doi: 10.1007/b79761_5. 79
	Christophe Berthon and Rodolphe Turpault. Asymptotic preserving hll schemes. Nu-
	merical Methods for Partial Differential Equations, 27:1396 -1422, 11 2011. doi:
	10.1002/num.20586. 91

  3.44 and 3.45. We consider only one intermediate state for the velocity b ˚and the pressure Π ˚, and two intermediate states for a, a L and a R. From the last equation of eq. 3.45, a is a strong Riemann invariant associated with the central wave. Therefore, a is a weak Riemann invariant associated with both waves l ´and l `. Then, a L " a L and a R " a R . After integration, eq.

							3.55
	becomes	4c ? 3	Π L c 2 ´b2 L	pb L	´b˚q `pΠ L	´Π˚q " 0.	(3.56)
	Similar computations across the wave l `give

Using eqs. 3.49 and 3.51, we have adb " dΠ (3.55) across the wave l ´.

  .62) with equality for smooth functions. This equation can be split into the diffusion step B t S r `Sr B x β ě 0,

										(3.63)
	followed by the transport step
								B t S r `βB x S r ě 0.	(3.64)
	Using S r " s and the first equation of eq. 3.6, eq. 3.63 rewrites
								B t s ě 0.	(3.65)
	Following Buet and Despres [2008], let us now prove that the entropy s is non
	decreasing in the semi-discrete case for the diffusion step. The semi-discrete diffusion
	step is
			$ ' ' ' ' ' & ' ' ' ' % ' ∆m i ∆m i i	dτ i dt dt de i dv i dt	ptq " β i`1 2 ptq " ´´q i`1 2 ptq ´βi ptqβ i`1 2 ´1 2 ptq ptq ´qi ptq " ´c2 ´qi `1 2 ptq ´qi ´1 2 ptq ´1 2 ∆m ptqβ i´1 2 ptq ¯. (3.66)
	Using eq. 3.66 in eq. 3.61, and after some algebra, we have
	∆m i θ i ptq	ds i dt	ptq " ´´q i`1 2	ptqβ i`1 2	ptq `βi ptqq i`1 2	ptq `qi ptqβ i`1 2	ptq ´qi ptqβ i ptq qi
										"
				´1 2	ptqβ i´1 2	ptq `βi ptqq i´1 2	ptq `qi ptqβ i´1 2	ptq ´qi ptqβ i ptq
			´qi	`1 2	ptq ´qi ptq ¯´β i ptq ´βi	`1 2	ptq ¯´´q	i´1 2	ptq ´qi ptq ¯´β i ptq ´βi	´1 2	ptq ¯.

  |Fr| cEr " 4c 3c 2 `β2 |β|. One can show that |Fr| cEr ď 1 if and only if |β| ď c or |β| ě 3c. For physical reasons, we only consider the case |β| ď c.

	From eq. 3.60, one can show that
	S r "	4 3	ca

Table 3 .

 3 2 -Summary of the different cases, according to the value of D and b.

Table 4 .

 4 Figure 4.1 -Evolution of the residual as a function of the number of iterations of Jacobi method, with different resolutions.

	.2)

1 -Number of cell-updates per second of Jacobi method, with different resolutions.

Table 4 .

 4 This process can be applied locally, only in the cells where it is needed. The whole algorithm, including these modifications, is presented in section 4.A. In the next section, we present numerical results to show the gain in computational time. 2 -Computation time to reach the steady state with the explicit solver, the implicit solver using Newton-Raphson method and the geometric multigrid algorithm different values for L max .

	.31)

Table 4 .

 4 3 -Computational time per V-cycle with different values of L max .

  Table5.1 -Computational time with both explicit and implicit solvers on CPU and GPU. With the implicit solver, the matrix is updated in a parallel or a non-parallel way. The implicit solver uses the AMG preconditioner.

https://gitlab.erc-atmo.eu/erc-atmo/ark-rt/tree/v1.0.0

d`Ph2h ˘T , where d is the dimension of the problem (d P t1, 2,

3u). Thanks to the factor 1 2 d , values on the coarse mesh are also convex combinations of values on the fine mesh. In the one one-dimensional case given by fig.4.3, the restriction operator

5.1 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

Remerciements

Remerciements

Je ne sais pas si j'ai réussi à transcrire ici à quel point j'ai apprécié ces trois années, grâce à toutes les personnes que j'ai pu rencontrées, y compris celles que je n'ai pas nommées ici.

with q " c 2 ´β2 3c 2 `β2 E r , where β " ||β||. Let us notice that β is a velocity and q is a pressure.

Let us introduce an specific energy e and a specific flux v such that # E r " e F r " v. (3.4) Despite the notation, v is not a velocity. Because of the similarly with Euler equations, we use the same vocabulary as hydrodynamics: we call a density. Let us notice that the admissible states E r ą 0 and ||Fr|| cEr ď 1 are equivalent to ą 0, e ą 0 and ||v|| ce ď 1. Let us first consider the case without source terms, with σ a " σ s " 0. Using eqs. 3.1, 3.3 and 3.4, the M 1 model (eq. 1.15) can be written as

B t `∇ ¨p βq " 0 B t p vq `∇ ¨p v b βq `c2 ∇q " 0 B t p eq `∇ ¨p eβ `qβq " 0.

(3.5) Equation 3.5 has the same form as the Euler equations. The main difference is the link between the momentum v and the velocity β: in the hydrodynamics case, it is linear (the momentum is ρu and the velocity is u), whereas it is not linear for the M 1 model (eq. 3.2). Let us now focus on the discretization of eq. 3.5.

Diffusion-transport splitting

Following Chalons et al. [2016], we use a splitting strategy to solve eq. 3.5. In the hydrodynamics case, the Euler system is split into an acoustic step and a transport step. Here, we split eq. 3.5 into a diffusion step and a transport step. The choice of the terms "diffusion" and "transport" will be explained later. Transport phenomena of the form β ¨∇ are separated from the other terms. We first solve the diffusion step $ ' & ' % B t ` ∇ ¨β " 0 B t p vq ` v∇ ¨β `c2 ∇q " 0 B t p eq ` e∇ ¨β `∇ ¨pqβq " 0, (3.6) followed by the transport step $ ' & ' % B t `β ¨∇ " 0 B t p vq `β ¨∇ p vq " 0 B t p eq `β ¨∇ p eq " 0.

(3.7) On one hand, in the diffusive limit, we have F r " 0 at first order. Because F r and β are colinear and E r `q ‰ 0, we also have β " 0 at first order. Therefore, the variables are not impacted by the transport step at first order. Furthermore, one can show that q " Er 3 at first order, which is equivalent to the P 1 closure relation. Let us notice the source terms play an important role in the diffusive limit, as shown in section 1. 4, therefore, the analysis done here is not whole. Recovering the diffusion equation eq. 1.24 brings into play velocity β at order 1. More details about the integration of source terms can be found in section 3.1.6.

On the other hand, in the transport regime, β " cn, where n is a unit vector aligned with the radiative flux, and q " 0. The diffusion step simplifies into B t ϑ `cϑ∇ ¨n " 0 for ϑ P t , v, eu. In particular, if n is constant, all variables remain constant during the diffusion step. They are then transported at velocity c by the transport step.

Using eq. 3.6, one can deduce the equation of evolution of velocity β and pressure q:

´c2 ´β2 3c 2 ´β2 β∇ ¨β `c2 ´β2 4q ∇q ´c2 ´β2 2q p3c 2 ´β2 q β pβ ¨∇qq " 0 B t q `4c 2 q 3c 2 ´β2 ∇ ¨β ´c2 ´β2 3c 2 ´β2 β ¨∇q " 0.

(3.8)

Detailed computations are done in section 3.A.

Let us now study the eigenstructure of eqs. 3.6 and 3.7. For the sake of simplicity, we focus on the one-dimensional case. We write λ ´, λ 0 and λ `the eigenvalues of the M 1 model, with λ ´ď λ 0 ď λ `. From González et al. [2007],

where n x is the x-coordinate of n " f f and χ 1 " dχ df . The eigenvalues for the transport system (eq. 3.7) are β x " λ 0 , with multiplicity 3. Let us now focus on the acoustic step (eq. 3.6). Equation 3.8 is introduced to ease computations. Equations 3.6 and 3.8 involve five eigenvalues: λ ´´λ 0 , 0, λ `´λ 0 , the multiplicity of the eigenvalue 0 is 3. The fields associated with 0 (resp. λ ˘´λ 0 ) are linearly degenerate (resp. genuinely nonlinear). See section 3.B for more details. Both diffusion and transport systems are hyperbolic.

We now focus on the resolution of each system in the one-dimensional case. The extension to higher dimensions is straightforward.

Diffusion step

From now on, because we consider only the one-dimensional case, we slightly change our notations: β (resp. v) is no longer the norm of the vector β (resp. v), but it is the x-coordinate of the vector.

Let us first write the diffusion system (eq. 3.6) using Lagrangian variables pτ, v, eq, where τ " 1 . We introduce a mass variable dm " pt n , xq dx where the time is frozen at instant t n . Equation 3.6 becomes

One can check that eq. 3.13 is indeed the solution of this system.

To develop an asymptotic preserving scheme, Buet and Despres [2008] replaced eq. 3.24 with

whose solution is

3.C Derivation of β ˚and q

For the sake of simplicity, we focus on the one-dimensional case. We do not present the extension to higher dimensions. Equations 3.6 and 3.8 become

Let us recall from section 3.B that the eigenvalues of eq. 3.42 are 0 with multiplicity 3 and Λ ˘" c 2 ´β2 β˘?3c with multiplicity 1. After some algebra, one can show that the two last equations from eq. 3.42 rewrite

(3.43)

3.C.1 Relaxed set of equations

We now follow Bouchut et al. [2010]. Let us first introduce a relaxation parameter λ. We also introduce an approximation of the velocity b « β, pressure Π « q, and eigenvalues l ˘« Λ ˘. Finally, we introduce a coefficient a which is an approximation of 4cq ? 3pc 2 ´β2 q . We use these approximations in the flux terms:

We add the corresponding evolution equations to the system:

We formally recover eq. 3.42 when λ Ñ 8. We mimic this behavior by setting λ to 0 in eq. 3.45. More details about this relaxation method can be found in section 3.E.1. We aim at building an exact solver for eqs. 3.44 and 3.45. This solver will be an approximate solver for eq. 3.42.

3.C.2 Eigenstructure

Let us now study the eigenstructure of eqs. 3.44 and 3.45. We write W " p , v, e, b, Π, l `, l ´, aq T , then eqs. 3.44 and 3.45 can be written in the matrix form

where

3.D Admissible states

Let us focus on the case without source terms. We follow Buet and Despres [2008]. To prove that the schemes developed in the chapter preserve the admissible states, we first prove that they are entropic. As done for the radiative energy and the radiative flux, let us introduce a specific entropy s such that S r " s.

3.D.1 Diffusion step

Using eq. 1. [START_REF] Mizuta | Hydrodynamic instability of ionization fronts in HiiRegions[END_REF] Let us recall some results obtained in section 1.B:

(3.60) Using eq. 3.60 and q " c 2 ´β2 3c 2 `β2 E r , eq. After some algebra, one can show that α 0 i ě 0, α ì ě 0, α í ě 0 and α 0 i `αì `αí " T i . θT i is a convex combination of θn in the neighbor cells. Similar computations can be done for the entropy, therefore s T i ě 0 and S T i " T i s T i ě 0. We can conclude as done in section 3.D.

3.E Derivation of the unsplit scheme

The goal of this section is to build an approximate Riemann solver for eq. 3.5. For the sake of simplicity, we focus on the one-dimensional case, but it can be extended to higher dimensions without any difficulty.

3.E.1 Relaxed set of equations

Let us introduce an approximation of the density r « , velocity b « β, radiative flux µ « v, energy ϕ « e and pressure Π « q and a relaxation parameter λ in the flux terms:

We add the corresponding evolution equations to the system:

(3.78)

The relaxed system is an approximation of the original system, which can be recovered in the limit λ Ñ 8. We solve it using an operator splitting technique. We first solve

, p e, ϕq , pβ, bq , pq, Πqu, followed by the transport system without source terms, with λ " 0. We end up with the following system:

The field is genuinely nonlinear.

3.E.3 Jump conditions

Let us now write the jump conditions to solve the Riemann problem associated with eq. 3.79. From eq. 3.77, we only need the variables r, µ, ϕ, b and Π in an intermediate state to update the variables of interest , v and e. We are not interested in the update of the relaxed variables r, µ, ϕ, b, and Π. The intermediate state for b and Π will be given by Buet and Despres [2008], therefore we focus on the intermediate state for r, ϕ and µ.

We introduce the artificial density variable Θ to rewrite eq. 3.79 in conservative variables:

(3.85)

The corresponding Rankine-Hugoniot jump condition for a discontinuity of speed

´D r s `rrbs " 0 ´D r vs `rbµs `c2 rΠs " 0 ´D r es `rbϕs `rbΠs " 0 ´D rΘrs `rΘbrs " 0 ´D rΘµs `rΘbµs " 0 ´D rΘϕs `rΘbϕs " 0 ´D rΘs `rbΘs " 0,

where r¨s is the difference from either side of the discontinuity. We introduce the following quantify:

to write the jump condition as

´D r s `rrbs " 0 ´D r vs `rbµs `c2 rΠs " 0 ´D r es `rbϕs `rbΠs " 0 m rrs " 0 m rµs " 0 m rϕs " 0 rms " 0.

(3.88)

We now have to consider several cases:

1. If D " 0, then eq. 3.88 becomes

rbrs " 0 rbµs `c2 rΠs " 0 rbϕs `rbΠs " 0 m rrs " 0 m rµs " 0 m rϕs " 0 rms " 0,

(3.90) 1.a. If D " 0 and b " 0, then we have m " 0 and rrs P R, rµs P R, rϕs P R.

1.b. If D " 0 and b ‰ 0, then rrs " 0, rµs " 0, rϕs " 0.

2.a. We first assume that rbs " 0. Then eq. 3.88 becomes

´D r s `b rrs " 0 ´D r vs `b rµs `c2 rΠs " 0 ´D r es `b rϕs `b rΠs " 0 m rrs " 0 m rµs " 0 m rϕs " 0 pb ´Dq rΘs " 0.

(3.91) 2.a.i. If D ‰ 0, rbs " 0 and D " b, then m " 0 and rrs P R, rµs P R, rϕs P R.

2.a.ii. If D ‰ 0, rbs " 0 and D ‰ b, then rrs " 0, rµs " 0, rϕs " 0. 3.2 summarizes the different cases.

3.E.4 Riemann problem solution

Let us consider a piecewise initial data defined by

(3.92)

Chapter 4

Nonlinear Jacobi method and geometric multigrid

Contents

By writing v " p¨¨¨, E i , F i , ¨¨¨q T P S N , where N is the number of cells, eq. 4.1 rewrites

3)

The operators L, D, and R contain the terms depending on `En`1

.1 can also be seen as

where A is a nonlinear operator and f P S N is a known vector.

From [START_REF] Pichard | Martin Mathematiques appliquees et calcul scientifique Bordeaux[END_REF], we solve eq. 4.4 using algorithm 1.

Algorithm 1 Nonlinear Jacobi method

Let us notice that if A is a linear operator, algorithm 1 simplifies into the classical Jacobi method for a tridiagonal matrix (see e.g., Saad 2003). [START_REF] Pichard | Martin Mathematiques appliquees et calcul scientifique Bordeaux[END_REF] also shows that this algorithm converges to the unique solution v n`1 because the operator A is contractant.

Preservation of the admissible states

We assume that v n i P S, and we show that v n`1 i obtained with algorithm 1 is also admissible, i.e., v n`1 i P S. Let us show that if v n`1,pkq P S N , then v n`1,pk`1q P S N . Using Proposition 5.1 in [START_REF] Pichard | Martin Mathematiques appliquees et calcul scientifique Bordeaux[END_REF], one can show that L and R are stable, i.e., L pvq P S and R pvq P S if v P S.

Because f i " v n i , it is admissible and f i `L pv i´1 q `R pv i`1 q P S. Finally, in this particular case, D " `1 `c ∆t h ˘I, where I is the identity matrix. So,

We have shown that D ´1, L, and R are stable. Therefore, v n`1,pk`1q i is admissible. By induction, the proof is complete.

Geometric multigrid (GMG)

Convergence problematic

Even though the algorithm described in section 4.1 does not include the treatment of source terms for the M 1 model, it is the early development of a time-implicit solver suitable for the study of astrophysical problems. Therefore, it is reasonable to investigate its performances with academic problems. The one considered here is the beam

Application to the HLL solver for the M 1 model

Let us now apply algorithm 5 to the HLL solver for the M 1 model, without source terms. The operator A is the one described in section 4.1, and we use algorithm 1 as smoother and coarse grid solver. The system to be solved on the coarse mesh is

As shown in section 4.1.3, if the initial guess for v 2h and the right-hand side A 2h `v2h ˘r 2h are admissible, then the solution obtained with algorithm 1 is also admissible. However, numerical experiments have shown that, in general, A 2h `v2h ˘`r 2h is not admissible, which leads to a non-admissible solution.

To tackle this issue, we follow Kifonidis and Müller [2012] and we introduce a pseudo-time τ . Instead of solving eq. 4.26, we look for the steady state in pseudo-time of the following equation:

When the steady state is reached, du 2h dτ " 0, and we recover eq. 4.26. Let us notice that the pseudo-time τ is completely independent of the physical time step ∆t.

Equation 4.27 is a (nonlinear) system of ordinary differential equations in the variable τ . We use notation similar to the physical time for the discretization in pseudotime. ∆τ is the interval between the current pseudo-time τ m and τ m`1 . We choose ∆τ such that u 2h is admissible. Using the definition of the residual, the right-hand side of eq. 4.27 becomes A 2h `v2h ˘`f 2h ´R2h h `Ah `vh ˘˘. We want to solve eq. 4.27 with a stable scheme, for all ∆τ . Therefore, we use a splitting strategy. For the scheme to be stable, the left-hand side has to be taken implicitly. This leads to solving eq. 4.27 as

´`u 2h ˘m`1 ¯" 0.

(4.28)

The first equation in eq. 4.28 is explicit in pseudo-time. We can always choose a value for ∆τ such that Ą u 2h is admissible. Let us notice that the right-hand side is fixed and the left-hand side is local to a cell.

The second equation in eq. 4.28 is implicit in pseudo-time. With arguments similar to section 4.1.3, `u2h ˘m`1 is admissible as soon as Ą u 2h is admissible.

Unfortunately, choosing ∆τ such that Ą u 2h is admissible can result in small values for ∆τ . To reduce the computational cost, we slightly change the algorithm. We do not use the same pseudo-time step in all cells and for the implicit step. We choose a pseudo-time step ∆τ im for the implicit step and we iterate the explicit step with another pseudo-time step, local to each cell, until ∆τ im is reached. This results in algorithm 6. We use this algorithm as smoother and coarse solver in algorithm 5.

The next question to be solved is how to choose ∆τ im . Numerical experiments have shown that a large value for ∆τ im allows reducing the computational cost per V-cycle, but can result in a very low decrease of the norm of the residual. To avoid such a result, we use an adaptive pseudo-time step. When the norm of the residual decreases fast enough, ∆τ im is increased. On the contrary, when the norm of the residual decreases slowly, ∆τ im is decreased. Explicit, CFL " 0.45 Newton-Raphson, CFL " 0.2 GMG, L max " 1, CFL " 2 000 GMG, L max P t2, 3, 4u, CFL " 2 000

Numerical results

We perform some verification tests to validate the performances of the algorithms presented in sections 4.1 and 4.2. Because we developed them without source terms, we focus on tests without opacity. We consider the same tests as in section 3.2.1. We compare our results with those obtained with a time-explicit HLL solver and a time-implicit HLL solver using a Newton-Raphson method (see section 2.1.3). We will compare the performances with different time steps. We write ∆t " CFL h c . For the time-explicit solver, one should have CFL ď 1 to respect the Courant-Friedrichs-Lewy (CFL) condition. The eigenvalues in the numerical fluxes (eq. 2.2) are fixed to ˘c.

We always use the same parameters for the geometric multigrid method. At the finest level, the number of iterations for the pre-and post-smoothers is ν 0 " 3. When using iterations in pseudo-time, the number of iterations in pseudo-time m is set to 3 and the number of iterations for the smoothers is ν l " 1, for l ‰ 0. The initial value for ∆τ im is 10 -3 . These parameters are chosen because they give reasonable performances in most cases and can easily be used in physical problems.

Beam

We first consider the same test as in sections 2.2.3 and 3.2.1.2. It is the same test as in González et al. [2007]; Richling et al. [2001]. We consider the propagation of a beam in the free-streaming regime. The domain is r´1, 1s ˆr´1, 1s and it is discretized with the same number of cells in the x and y-directions. The initial temperature is T 0 " T r " 300 K and the initial radiative flux is F r " 0. At time t " 0, a beam is introduced with T g " T r " 1 000 K with an angle of 45 ˝at x " ´1 and y P r´0.875, ´0.75s.

The domain is first discretized with 129 ˆ129 cells. Figure 4.10 shows a horizontal cut at the middle height once the steady state is reached. With the Newton-Raphson method, we are not able to reach large time steps. Indeed, we have to use a smaller

Discussion and conclusion

In this chapter, we have first presented a Jacobi method to solve the nonlinear system arising from the discretization of the M 1 model without source terms with a time-implicit HLL solver. The admissible states are preserved, even with large CFL numbers. This method is iterative, the convergence rate decreases when the resolution increases. To tackle this issue, we use a nonlinear geometric multigrid algorithm. The Jacobi method first described is used as smoother and coarse grid solver. However, this algorithm relies on the residual equation, which does not preserve the admissible states. Instead of solving this equation, we introduce a pseudo-time, and we look for a steady state in pseudo-time. Numerical experiments have shown the good performances obtained with this algorithm.

Source terms

In order to study astrophysical problems, the integration of source terms is indispensable. It needs to be done carefully, regarding the asymptotic preserving property already discussed in chapters 2 and 3. The schemes developed in chapter 3 are asymptotic preserving and preserve the admissible states. The HLL solver used here could be replaced by the all-regime schemes presented in chapter 3. This will change the nonlinear operator A. It requires to choose carefully which terms are taken at iteration k and k `1 in algorithm 1. The geometric multigrid algorithm presented in section 4.2 does not rely on the form of the nonlinear operator, as long as the iterations of the Jacobi method preserve the admissible states. Therefore, the geometric multigrid method can be used to study astrophysical problems with opacity as soon as the Jacobi algorithm is extended to the all-regime scheme.

Performances

Although performance results shown in section 4.3 are promising, some choices were made and others were not explored. For example, only a V-cycle is used, but there exist other possibilities: W-cycle, F-cycle,. . . The restriction and prolongation operators can also have an impact on the performances.

Another well-known method to reduce the number of iterations performed by the Jacobi method is to replace it with a Gauss-Seidel algorithm. Even though the tests presented in section 4.3 are obtained with a sequential code, we aim at using a parallel implementation to study astrophysical situations. In the linear case, the Jacobi method is known to be easier to make parallel than the Gauss-Seidel method.

Bibliography

4.A Overall algorithm

Algorithm 7 presents the whole algorithm to compute v n`1 with v n known. When we presented the V-cycle in section 4.2.2, we used a recursive form. We present here a version based on loops.

Algorithm 7 Overall algorithm

Choose a maximum number of levels L max κ " 0 f h " v n v h p0q " v n . From now on, the exponent is the mesh size, it is no longer the time while r pκq ą ε do for l " 0, ¨¨¨, L max ´1 do if l " 0 then Pre-smoother: solve A h ´uh pκq ¯" f h with ν 0 iterations of the nonlinear Jacobi method (algorithm 1), with initial guess v h pκq Compute residual: r h pκq " f h ´Ah ´uh pκq ēlse Pre-smoother: compute steady state in pseudo-time of

A 2 l h ´v2 l h pκq ¯`r 2 l h pκq using algorithm 6, with ν l iterations of the nonlinear Jacobi method , with initial guess v 2 l h pκq Compute residual:

for Solver at coarsest level: compute steady state in pseudo-time of

pκq using algorithm 6, with ν Lmax iterations of the nonlinear Jacobi method , with initial guess v 2 Lmax h pκq for l " L max ´2, ¨¨¨, 0 by step ´1 do Prolongation:

each cell l i in the grid at level l do if ´v2 l h pκq ¯li is not admissible then use eq. 4.32 to update ´v2 l h pκq ¯li end if end for if l " 0 then Post-smoother: solve A h ´uh pκq ¯" f h with ν 0 iterations of the nonlinear Jacobi method (algorithm 1), with initial guess v h pκq else Post-smoother: compute steady state in pseudo-time of

A 2 l h ´v2 l h pκq ¯`r 2 l h pκq using algorithm 6, with ν l iterations of the nonlinear Jacobi method , with initial guess v 2 l h pκq end if end for

CrsMatrix

Using the scheme presented in chapter 2, there are only a few nonzero elements per row. Using a sparse matrix is essential to save both computational time and memory. The storage of the matrix is done with a compressed sparse row (CSR) or compressed row storage (CRS) format (see e.g., [START_REF] Buluç | Parallel sparse matrix-vector and matrix-transpose-vector multiplication using compressed sparse blocks[END_REF]). The class Tpetra::CrsMatrix implements such storage. The rows of the matrix are distributed across the MPI processes. All the nonzero elements of a row are owned by the same process. as multi-cores and GPUs, thanks to the package Kokkos. Finally, we discussed the impact of the choice of the preconditioner and the MPI domain decomposition on a physical application.

Vectorization

We did not explore all the capabilities of Kokkos. In particular, Kokkos provides tools for vectorization through hierarchical parallelism. Some CPUs can perform operations on multiple floating points with the same number of cycles as the scalar version. Enabling vectorization through Kokkos should improve the performances [START_REF] Padioleau | Development of "all-régime[END_REF].

Implementation of other numerical schemes

Numerical schemes developed in chapters 3 and 4 are only implemented in sequential codes. Parallelism in the code ARK, which implements an all-regime solver for hydrodynamics, is handled through Kokkos and MPI. Because our all-regime scheme (chapter 3) has the same structure as the scheme for hydrodynamics, its implementation in ARK should be easily done. However, we have not done it because we aim for an implicit solver.

The Jacobi method, presented in section 4.1, can easily be made parallel. It could also be done with Kokkos and MPI. However, the geometric multigrid algorithm (section 4.2) raises the problem of load balancing when the coarser grids, with only a few unknowns, are reached. The Tpetra Project Team. The Tpetra Project Website. 119

Bibliography

The Teuchos Project Team. The Teuchos Project Website. 119 target different architectures, such as multi-cores, many-cores, and GP-GPUs. We have discussed the performances obtained thanks to Kokkos and the different preconditioners provided by Trilinos. We have also investigated the impact of the MPI domain decomposition and the preconditioner on the stability of the ionization front in a massive pre-stellar dense core. Several perspectives arise from this work. To study astrophysical problems, an implicit and asymptotic preserving scheme is required. It should also preserve the admissible states. We have presented several schemes, each of them with two of these three properties. Developing a scheme with all of these properties, using the schemes presented in chapters 3 and 4, has not been done because of time constraints.

Such a scheme could be altered with a well-balanced modification of the source terms to capture properly the steady state with a discontinuity of opacity. In chapter 2, we have presented such a modification. For performance reasons, we did not take it into account in the schemes presented in chapter 3.

So far, the schemes presented in chapters 3 and 4 are implemented in a sequential C++ code. To study astrophysical problems, high performance computing is mandatory. We can achieve portability and performance portability with the library Kokkos.

Inspired by the resolution of the Poisson equation for self-gravity [Guillet and Teyssier, 2011], the geometric multigrid algorithm could be extended in the context of Adaptive Mesh Refinement (AMR), but is beyond the scope of this thesis.

Bibliography