Dans cette thèse, nous nous intéressons à la résolution des équations SPN du transport de neutrons au sein des coeurs de réacteurs nucléaires à eau pressurisée. Ces équations s'écrivent naturellement sous une forme mixte, les inconnues sont le flux de neutrons et le courant de neutrons, mais peuvent aussi s'écrire sous une forme primale, où le flux est la seule inconnue du problème. Les équations SPN du transport des neutrons forment un problème aux valeurs propres généralisé. Dans notre étude nous commenons par le problème source associé et ensuite nous étudions le problème aux valeurs propres. Un coeur de réacteur est composé de différents milieux: le combustible, le fluide caloporteur, le modérateur... à cause de ces hétérogénéités de la géométrie, le flux solution du problème source peut être peu régulier. Nous montrons que ce problème est bien posé sous ses formes mixte et primale. Nous trouvons aussi une estimation d'erreur a priori pour l'approximation de la solution par la méthode des éléments finis dans les deux formes du problème dans le cas où la solution est peu régulière. Pour le problème aux valeurs propres sous sa forme primale, nous appliquons la théorie déjà existante pour l'approximation des problème aux valeurs propres. Dans le cas mixte, les théories déjà développées ne s'appliquent pas. Nous proposons ici une nouvelle méthode pour étudier la convergence de la méthode des éléments finis mixtes pour les problèmes aux valeurs propres. Pour les solutions peu régulières, la montée en ordre de la méthode des éléments finis n'améliore pas l'approximation du problème, il faut raffiner le maillage aux alentours des singularités de la solution. La géométrie des coeurs de réacteur se prête bien aux maillages cartésiens, mais leur raffinement augmente vite leur nombre de degrés de liberté. Pour palier à cette augmentation, nous proposons ici une méthode de décomposition de domaine qui permet d'utiliser des maillages globalement non-conformes. Finalement, nous appliquons cette méthode pour un cas test industriel de réacteurs à eau pressurisée. vi

Introduction

This thesis, done at the "Commissariat à l'Énergie Atomique et aux Énergies Alternatives" (CEA) and at the "École Nationale Supérieur des Techniques Avancées" (ENSTA), aims at proposing the numerical analysis of a finite element discretization for nuclear reactor core simulation. Moreover, this work proposes a non-conforming domain decomposition method.

In order to simulate nuclear reactor core, one has to solve the neutron transport equation. This equation describes the dependence of the neutron density on 7 variables which are :

-three for the space; -two for the direction; -one for the energy; -one for the time.

Using the multigroup discretization (for the energy dependence) and projecting the equation on the spherical harmonics (for the direction dependence), one obtains the multigroup SP N equations. In the neutronic platform APOLLO3 R the solver MINOS solves numerically these equations with a Finite Element method. Due to the high number of isotopes and the complicated dependence of the cross sections on the direction and the energy, the nuclear data can be massive. In order to reduce the number of these data, the cross sections are pre-processed. One step of this preprocessing is a homogenization, which transforms the cross sections into coefficients, which are piecewise regular. We call a material an area of the reactor where all the homogenized cross sections are regular. In figure 1 we present an example of the pre-processed geometry of a part of a nuclear reactor core, where each color corresponds to a material. In this realistic configuration, i.e. the cross sections are piecewise regular, the points where three or more materials intersect are allowed (crosspoints). Therefore, the solution of the multigroup SP N equations has a low-regularity. The numerical analysis done in this manuscript addresses the case of low regularity solution. In MINOS, the triangulation used for the Finite Element Method is Cartesian and is constructed from the material geometry such that each element is included in exactly one material cell. In the left component of figure 2, we show the coarser triangulation of the geometry given in figure 1. Due to the non-conformity of the geometry, lots of CONTENTS elements with large aspect ratio can appear. For this geometry, the coarser triangulation is composed of 162 elements. In order to reduce the number of elements, one can use a non-conforming triangulation, as shown in the right component of figure 2. This latter triangulation has 114 elements, with better aspect ratio.

Conforming triangulation

Non-conforming triangulation We propose a domain decomposition method in order to treat non-conforming triangulation with finite element method. We also carry out the numerical analysis of this method in the case of low-regularity solution. In this work, we focus on the numerical analysis of the method. Moreover, we describe the algorithmic aspect of this method. The sequential implementation of this method has been done in MINOS. We do not address the high performance computing aspect but the implementation of this domain decomposition method could be parallelized [START_REF] Lathuilière | Méthode de Décomposition de Domaine Pour Les Équations Du Transport Simplifié En Neutronique[END_REF].

After the introduction of the physical model in chapter 1, this manuscript is composed of three parts. First, in part I, we look at the continuous problem of the multigroup SP N neutron transport equations. then in part II, we focus on the discrete version of these equations with some finite element methods. Finally, in part III, we show some numerical applications of the resolution of the multigroup SP N neutron transport equations. More precisely, in the first part, in chapter 2 we study the diffusion approximation under its primal formulation and its mixed formulation. Then, in chapter 3 we generalize the previous results to the multigroup SP N neutron transport equations.

In the second part, in chapter 4, we carry out the numerical analysis of a H 1 -conforming finite element method to solve the multigroup SP N neutron transport equations while in chapter 5 we carry out the numerical analysis of a mixed finite element method. Then, in chapter 6, we introduce the L 2 -jump domain decomposition method.

In the last part, in chapter 7, we describe the algorithm to solve the eigenvalue problem given by the multigroup SP N neutron transport equations. Then, we illustrate the theoretical result and the L 2 -jump domain decomposition method with two testcases. Finally, in chapter 8, we derive an a posteriori error estimate for the mixed neutron diffusion problem.

Chapter 1

Modelling

A nuclear power plant produces electricity thanks to a coolant heated by a nuclear reactor which creates steam to run a turbine so that this one produces electricity. The nuclear reactor core is a delimited area where a fission chain reaction occurs. When a neutron collides a fissile nucleus there is a probability that a fission occurs. The fission emits new neutrons which can induce new fissions, thereby generating a chain of fission events. Moreover, each fission is accompanied by the release of energy which heats the coolant.

The model used to describe the physics of the nuclear reactor core is the neutron transport equation. Classically the resolution of this equation is decomposed in two steps. The first step consists in homogenizing nuclear data on each fuel rod or assembly. The second step is the core calculation which approximates the solution of the neutron transport equation with the homogenized nuclear data.

The "Commissariat à l'Énergie Atomique et aux Énergies Alternatives" (CEA) is a French public government-funded research organisation in the areas of energy, defence and security, information technologies and health technologies. One aim of the CEA energy research is to develop some neutronic simulation tools for nuclear power plants such that for instance APOLLO3 R . More precisely, APOLLO3 R is a shared platform with Électricité de France (EDF) for the nuclear reactor core physics. The purpose of this platform is to simulate the nuclear reactor core behaviour by solving the neutron transport equation.

Our work takes part in the core calculation step inside APOLLO3 R . The solution of the neutron transport equation depends on the neutron position in space, the neutron motion direction and the neutron energy (only the steady state case is studied in this manuscript so that we do not consider the time dependence), each of these dependencies are discretized with different methods. We consider here the multigroup approximation for the energy, the simplified spherical harmonics (SP N ) and a finite element method for the space. When in the steady state case there is no sources, the multigroup SP N equations correspond to an eigenvalue problem. We are looking for the fundamental mode of this eigenvalue problem.

In order to find this fundamental mode the multigroup SP N neutron transport eigenvalue problem is transformed into its variational form. As this equation is naturally written under its mixed form, we approximate the solution with Raviart-Thomas-Nédélec (RTN) finite element method. The approximation of eigenvalue problem with finite element has been studied by Osborn et al. [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF][START_REF] Babuška | Eigenvalue Problem[END_REF] in general cases, and more particularly for eigenvalue problem under mixed form by Boffi et al. [START_REF] Boffi | On the convergence of eigenvalues for mixed formulations[END_REF][START_REF] Boffi | On the problem of spurious eigenvalues in the approximation of linear elliptic problems in mixed form[END_REF][START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF]. In this work, we study the approximation of the multigroup SP N neutron transport equation with finite element. Moreover, we propose here a domain decomposition method for mixed RTN finite element method which can be used on non-conforming triangulations.

But first of all, we establish the neutron transport equation which models the evolution of the neutron population. Let R represents the delimited area of the reactor. The neutron population is described thanks to the neutron density N (x, Ω, E, t) in cm -3 .Mev -1 .sr -1 , in the phase space (x, Ω, E, t) representing:

• the position of the measure x = (x, y, z) T ∈ R 3 ,

• the direction of motions of the neutrons Ω ∈ S 2 ,

• the kinetic energy of the neutrons E ∈ R + ,

• the time of the measure t ∈ R + .

The velocity of the neutrons is given by the vector

v = v(E)Ω,
where the link between the velocity norm and the kinetic energy is

E = 1 2 mv 2 .
The quantity N (x, Ω, E, t)dxdΩdE can be interpreted as the number of neutrons in a infinitesimal volume dxdΩdE around (x, Ω, E) at time t.

The interaction between a neutron and a nucleus leads to different events as fission, absorption and scattering. Neutron cross sections provide a quantitative measure of the probability for various types of neutron-nuclear reactions to occur. For the scattering interaction, the general concept of neutron cross section, is useful as it links the incident neutron velocity to the one of the neutron after the collision. Thus, the microscopic differential scattering cross section represents the probability of an incident neutron with an energy E and a direction Ω to be scattered by a nucleus of the isotope i with an energy E and a direction Ω . This probability is denoted by

σ s,i (Ω → Ω , E → E ) in cm 2 .sr -1 .Mev -1 .
For other reactions, like for the fission, we introduce σ x,i (Ω, E) in cm 2 , to characterize the probability that a neutron with a kinetic energy E and an incident direction Ω incident upon a nucleus of the isotope i will produce an event of type x. These possible events are fission, scattering, radiative capture and (n, 2n) scattering. Another interpretation of this quantity is the cross sectional area presented by the target nucleus to the incident neutron. Particularly, we define the fission microscopic cross section σ f,i (Ω, E).

Finally, we define the total microscopic cross section as

σ t,i (Ω, E) = S 2 ∞ 0 σ s,i (Ω → Ω , E → E )dE dΩ + σ f,i (Ω, E) + σ ,i (Ω, E), (1.1)
where σ (Ω, E) represents all possible reactions but the scattering and the fission, as sterile capture and radiative capture. The total microscopic cross section is the probability that a neutron with kinetic energy E and a direction Ω will interact with a nucleus of the isotope i.

All the microscopic cross sections are positive, for any reaction x and isotope i:

σ x,i ≥ 0.
Moreover, it holds that σ ,i (Ω, E) > 0.

(1.

2)

The dependence of the microscopic cross sections over the direction of the neutron is much weaker than the one over the neutron energy [START_REF] Duderstadt | Nuclear Reactor Analysis[END_REF]. A nucleus has only discrete stable energy levels. Then, the reaction during a neutron-nucleus collision depends of the new energy level of the nucleus which is mainly defined by the neutron energy. In nuclear reactor engineering, the media are modelled to be isotropic:

Hypothesis 1.1 (Isotropic media). Every medium in the reactor is isotropic, so that all the microscopic cross sections except the scattering do not depend on Ω. Moreover, the scattering microscopic cross section σ s (Ω → Ω, E → E) does not depend on Ω and Ω but only on the cosine of the angle between these two directions c = cos ΩΩ .

In a nuclear reactor core, the materials are composed by several isotopes. As the microscopic cross sections consider only the probability of a neutron which collides a nucleus of isotope i to produce an event of type x, we need macroscopic cross section Σ x to described the probability of a neutron which collides a nucleus in the materials to produce an event of type x. The macroscopic cross sections take into account the composition of the material. Thus, we define the probability that a neutron with kinetic energy E and a direction Ω incident upon a nucleus of a material will produce an event of type x by unit length by the formula

Σ x (x, •, t) = i σ x,i (•)C i (x, t), (1.3)
where C i is the isotope concentration (in cm -3 ) in the material. The macroscopic cross sections are expressed in cm -1 (in cm -1 .sr -1 .Mev -1 for the scattering macroscopic cross section). More particularly, we introduce:

• Σ s (x, c, E → E , t) the scattering macroscopic cross section;

• Σ f (x, E, t) the fission macroscopic cross section;

• Σ t (x, E, t) the total macroscopic cross section which can be interpreted as the inverse of the neutron mean free path in the material.

• Σ (x, E, t) the macroscopic cross associated to all possible reactions but the scattering and the fission. From inequality (1.2), it holds:

∀(x, E, t) ∈ R × R * × R + , Σ (x, E, t) > 0.
(1.4)

Equation (1.1), which links the microscopic cross sections, can be extend to the macroscopic cross sections. By multiplying equation (1.1) by the isotope concentration and summing over the isotopes, one obtains:

Σ t (x, E, t) = S 2 ∞ 0
Σ s (x, c, E → E , t)dE dc + Σ f (x, E, t) + Σ (x, E, t).

(1.5)

The product of the velocity and the density, v(E)N (x, Ω, E, t) represents the number of neutrons with a velocity v(E)Ω passing through an infinitesimal surface perpendicular to Ω located in x per energy unit, solid angle unit and time unit and is denoted by: ψ(x, Ω, E, t) = v(E)N (x, Ω, E, t).

(1.6)

This quantity is called the neutron angular flux and is measured in cm -2 .sr -1 .Mev -1 .s -1 . The total reaction rate is defined:

τ t (x, Ω, E, t) = Σ t (x, E, t)ψ(x, Ω, E, t),
and measures the total number of neutron-nucleus interactions. When a nucleus fission occurs, 2 or 3 neutrons are emitted. The fission yield is the number of neutrons emitted by a fission generated by a neutron with an energy E ν p (E).

The fission spectrum described the probability of a neutron generated by a fission to have an energy E and is denoted by:

χ p (E) in Mev -1 ,
where E is the new neutrons energy. Moreover, a nucleus, produced by a fission, can be transformed by radioactive decay, and the result of this decay can emit some neutrons by spontaneous fission decay. Those neutrons are called delayed neutrons. The radioactive decay constant associated to the radioactive decay of an isotope i is denoted by:

λ d,i in s -1 .
The number of neutrons emitted by radioactive decay of a fission product

ν d .
The energy distribution of the neutron emitted by radioactive decay is denoted by:

χ d (E) in Mev -1 ,
where E is the new neutron energy. We define the macroscopic radioactive decay associated to the radioactive decay phenomenon:

Λ d (x, t) = i λ d,i C i (x, t) in cm -3 .s -1 .
Let us note that nuclear reactions not only have an influence on neutron population; they also induce variations in the population of atomic nuclei.The evolution in concentrations 1.1 : The Neutron Transport Equation of the various isotopes is governed by the Bateman equations, taking into account atomic nuclei generations and disappearances through nuclear reactions and radioactive decay process [START_REF][END_REF].

In section 1.1 we obtain the neutron transport equation, then in section 1.2, we obtain the Bateman equations, in section 1.3, we derive its steady-state. Next, in section 1.4, we discretize it in energy. Finally, in section 1.5, we treat the angular dependence.

All the nuclear quantities defined in this chapter, above and below, are summarized in appendix A. For more details about nuclear data and the neutronics, one can refer to [START_REF] Duderstadt | Nuclear Reactor Analysis[END_REF][START_REF] Bussac | Traité de Neutronique[END_REF].

The Neutron Transport Equation

To estimate the power of a nuclear reactor core, one must study the neutron distribution.

Establishing the neutron balance in an infinitesimal volume lets us model the evolution of the angular neutron flux inside the reactor core. In an infinitesimal volume dxdΩdE around (x, Ω, E, t), the neutron density variations are given by

dN dt (x, Ω, E, t) = dx dt • grad x N (x, Ω, E, t) + dv dt • grad v N (x, Ω, E, t) + ∂ t N (x, Ω, E, t)
From the fundamental law of the mechanics, we know that the acceleration of a neutron is proportional to the sum of the forces applied to him. As there is no other forces applied to a free neutron than the gravity which we neglect, its acceleration is null, dv dt = 0. Moreover, from the definition of the speed, it stands v = dx dt . Then, the neutron density evolution can be written as:

dN dt (x, Ω, E, t) = v(E)Ω • grad x N (x, Ω, E, t) + ∂ t N (x, Ω, E, t).
Using the definition of the neutron angular flux, equation (1.6), it holds:

dN dt (x, Ω, E, t) = 1 v(E)
∂ t ψ(x, Ω, E, t) + Ω • grad x ψ(x, Ω, E, t).

(1.7)

The variations of the neutron density comes from different phenomena:

• Neutron loss due to the interactions with some nucleus:

-Σ t (x, E, t)ψ(x, Ω, E, t), (1.8a)

• Neutron transfer after collision:

∞ 0 S 2 Σ s (x, c, E → E, t)ψ(x, Ω , E , t)dΩ dE , (1.8b)
• Fission sources:

χ p (E) 4π ∞ 0 S 2 ν p (E )Σ f (x, E , t)ψ(x, Ω , E , t)dΩ dE , (1.8c)
• Delayed neutron sources:

χ d (E) 4π ν d Λ d (x, t), (1.8d)
• External sources: S ext (x, Ω, E, t).

(1.8e)

The particle derivative of the neutron density (1.7) is equal to the sum of all the neutron variations (1.8a)-(1.8e). This equality gives us the neutron transport equation:

1 v(E)
∂ t ψ(x, Ω, E, t) + Ω • grad x ψ(x, Ω, E, t) = -Σ t (x, E, t)ψ(x, Ω, E, t)

+ ∞ 0 S 2 Σ s (x, c, E → E, t)ψ(x, Ω , E , t)dΩ dE + χ p (E) 4π ∞ 0 S 2 ν p (E )Σ f (x, E , t)ψ(x, Ω , E , t)dΩ dE + χ d (E) 4π Λ d (x, t) + S ext (x, Ω, E, t).
(1.9)

The neutron-nucleus interactions also changes the composition of the materials inside the reactor core. To model this evolution, we introduce to our model some equations which describe the evolution of isotope concentration inside the reactor core.

The Bateman Equations

After a neutron-nucleus collision, an isotope i can change into an isotope i . The probability that a neutron with an energy E colliding an isotope i produces an events x, except a scattering reaction, and change the nucleus to an isotope i is denoted by:

σ * i→i ,x (E).
Moreover, we have the following relation:

σ i,x (E) = i σ * i→i ,x (E).
(1.10)

For the scattering event, the probability that a neutron with an energy E and a direction Ω colliding an isotope i is scattered with an energy E and a direction Ω and changes the nucleus to an isotope i is denoted by:

σ * i→i ,s (c, E → E ).
Moreover, we have the following relation:

σ i,s (c, E → E ) = i σ * i→i ,s (c, E → E ).
(1.11)

The probability of an isotope i nucleus changes into an isotope i after colliding a neutron at a position x at a time t is given by:

ζ i→i (x, t) = x =s ∞ 0 S 2 σ * i→i ,x (E)ψ(x, Ω, E, t)dEdΩ + ∞ 0 S 2 ∞ 0 S 2 σ * i→i ,s (c, E → E )ψ(x, Ω, E, t)dEdΩdE dΩ .
Therefore, the number of isotope i changed into an other isotope after colliding a neutron at position x and at time t is given by

i =i ζ i→i (x, t)C i (x, t),
(1.12) and the number of isotope i which appears by the transformations of an other isotope after colliding a neutron at a position x and a time t is given by

i =i ζ i →i (x, t)C i (x, t).
(1.13) Moreover, each isotope i nucleus can be changed into an isotope i by radioactive decay, the probability of this to occur is given by λ i→i . Then, the number of isotope i disappearing by radioactive decay is

i =i λ i→i C i (x, t), (1.14)
and the number of isotope i appearing by radioactive decay is

i =i λ i →i C i (x, t). (1.15)
Therefore, the evolution of the concentration of isotope i at a position x at a time t is obtained by summing equations (1.12)-(1.15) and reads:

∂ t C i (x, t) = i =i (ζ i →i (x, t) + λ i →i )C i (x, t) - i =i (ζ i→i (x, t) + λ i→i )C i (x, t).
(1.16)

In the case of pressurized water reactor equations (1.16) can be simplified [START_REF] Reuss | Précis de Neutronique[END_REF][START_REF] Chauvet | Méthode Multi-Échelle Pour La Résolution Des Équations de La Cinétique Neutronique[END_REF]. Indeed, we consider only the evolution of precursors. The precursors are radioactive isotopes that emit neutrons with a given delay. The precursors disappear only by radioactive decay and appears by fission. Then, the evolution of the concentration of a precursor i reads:

∂ t C i (x, t) = i =i ζ i →i,f (x, t)C i (x, t) - i =i λ i→i C i (x, t), (1.17)
where

ζ i →i,f (x, t) = ∞ 0 S 2 σ * i →i,f (E)ψ(x, Ω, E, t)dEdΩ.
(1.18)

The system composed of equation (1.16) for all the isotopes is called the generalized Bateman equations. To model the physics of a nuclear reactor, the neutron transport equation and the generalized Bateman equations have to be taken into account together.

The Steady State Equation

Physical Point of View

Here we are interested in the steady state of the reactor core without any external source, S ext = 0, when the chain reaction has already started. Thus, it stands ∂ t ψ(x, Ω, E, t) = 0, and the neutron transport equation reads:

Ω • grad x ψ(x, Ω, E) + Σ t (x, E)ψ(x, Ω, E) = ∞ 0 S 2 Σ s (x, c, E → E)ψ(x, Ω , E )dΩ dE + χ p (E) 4π ∞ 0 S 2 ν p (E )Σ f (x, E )ψ(x, Ω , E )dΩ dE + χ d (E) 4π ν d Λ d (x, t). (1.19)
Moreover, it stands, for all precursors i, ∂ t C i (x, t) = 0, therefore, Bateman equations (1.17) become:

0 = i =i ζ i →i,f (x, t)C i (x, t) - i =i λ i→i C i (x, t).
(1.20)

By summing equation (1.20) over the precursors, one obtains:

i i =i ζ i →i,f (x, t)C i (x, t) = i i =i λ i→i C i (x, t). (1.21)
Using the definition of ζ i →i,f (1.18), equation (1.10), noting that for all isotope i it stands σ * i→i,f = 0 and the definition of Σ f (1.3), we obtain:

∞ 0 S 2 Σ f (x, E)ψ(x, Ω, E, t)dEdΩ = i i =i λ i→i C i (x, t). (1.22)
Moreover, we have, for all precursors i:

i =i λ i→i = λ d,i .
Therefore, using the definition of Λ d , equation (1.22) becomes:

∞ 0 S 2 Σ f (x, E)ψ(x, Ω, E, t)dEdΩ = Λ d . (1.23)
We denote the number of neutrons emitted by fission and radioactive decay by

ν(E) = ν p (E) + ν d .
We also define β d (E), the quantity such that

ν d = β d (E)ν(E) and ν p (E) = (1-β d (E))ν(E).
Using equation (1.23) in equation (1.19), one obtains:

Ω • grad x ψ(x, Ω, E) + Σ t (x, E)ψ(x, Ω, E) = ∞ 0 S 2 Σ s (x, c, E → E)ψ(x, Ω , E )dΩ dE + 1 4π ∞ 0 S 2 χ(E, E )ν(E )Σ f (x, E )ψ(x, Ω , E )dΩ dE , (1.24)
where for all energy E and E

χ(E, E ) = (1 -β d (E ))χ p (E) + β d (E )χ d (E).

: The Steady State Equation

Note that equation (1.24) is a homogeneous problem and that ψ = 0 is a solution to it. For a given reactor, the flux ψ = 0 is in general the unique solution to this problem. This seems to be in contradiction with the real physical situation in which the stationary reactor has a non zero flux. But in practice, the core never fully reaches stable conditions, thus the flux does not exactly satisfy equation (1.24). Therefore, we look for a non zero flux that could be considered as a representation of the system under nearly steady conditions. For this purpose, we relax equation (1.24) by altering the fission yield term in (1.24) by a factor λ -1 , ν(E ) becomes λ -1 ν(E ). With this correction on the fission sources, one obtains the following generalized eigenvalue problem:

Ω • grad x ψ(x, Ω, E) + Σ t (x, E)ψ(x, Ω, E) = ∞ 0 S 2 Σ s (x, c, E → E)ψ(x, Ω , E )dΩ dE + λ -1 1 4π ∞ 0 S 2 χ(E, E )ν(E )Σ f (x, E )ψ(x, Ω , E )dΩ dE .
(1.25)

In equation (1.25), the coefficient λ -1 is the eigenvalue of the problem. The eigenvalue λ represents the multiplication factor of neutrons emitted by fission. More precisely, λ is the ratio of the number of neutrons emitted by fission over the neutrons lost by absorption or leakage. In order to control the fission chain reaction in a reactor core, the control rods, which are composed of an neutron absorbant medium, are inserted or removed from the core. The resolution of the eigenvalue problem (1.25) shows us if the control rods need to be inserted or removed from the core in order to maintain the fission reaction chain.

Remark 1.2. One can choose to alter with a multiplicative coefficient other terms in equation (1.25) instead of the fission source term. One can refer to [START_REF] Veladre | A Comparison of the Eigenvalue Equations in k, α, λ and γ in Reactor Theory[END_REF] for the other alternatives. In those other cases, the eigenvalue λ is still multiplication factor but it does not represent the same ratio. For instance, if one decide to change the scattering source term which the fission source term, λ is the ratio of neutrons emitted by collision over the loss.

We call a physical solution of equation (1.25) any positive solution, as a matter of fact the angular flux can not be negative on some part of the reactor.

Mathematical Point of View

Thanks to Krein and Rutman in [START_REF] Krein | Linear Operators Leaving Invariant a Cone in a Banach Space[END_REF], the following theorem allows us to know which solutions are physical.

Theorem 1.3 (Krein-Rutman reported in [START_REF] Brezis | Analyse Fonctionelle -Théorie et Applications[END_REF]). Let E be a Banach space and C be a convex cone of E with 0 as apex. We assume that C is closed, C = ∅ and C ∩ -C = {0}. Let T c be a compact operator on E such that T c (C\{0}) ⊂ C, and we denote σ(T c ) the spectrum of T c . Thus, there exists u in C and λ > 0 such that T c u = λu; moreover λ is the unique eigenvalue of T c associated to an eigenvector in C and the multiplicity of λ is one. Finally, it stands:

λ = max µ∈σ(Tc) |µ|,
One can apply this theorem to the inverse transport operator, as done in [START_REF] Dautray | Analyse Mathématique et Calcul Numérique Pour Les Sciences et Les Techniques[END_REF]. Indeed, the set of all the square integrable functions over

R × S 2 × R + , L 2 (S 2 × R + , L 2 (R))
, is a Banach space and the set of all its positive functions is a cone which satisfies all the hypothesis of theorem 1.3. The steady state neutron transport operator T from

L 2 (S 2 × R + , H 1 (R)) to L 2 (S 2 × R + , L 2 (R)) is defined as, for all ψ in L 2 (S 2 × R + , H 1 (R)): T (ψ)(x, Ω, E) = Ω • grad x ψ(x, Ω, E) + Σ t (x, E)ψ(x, Ω, E) - ∞ 0 S 2 Σ s (x, c, E → E)ψ(x, Ω , E )dΩ dE . The fission operator F from L 2 (S 2 × R + , H 1 (R)) to L 2 (S 2 × R + , L 2 (R)) is defined as, for all ψ in L 2 (S 2 × R + , H 1 (R)): F (ψ)(x, Ω, E) = 1 4π ∞ 0 S 2 χ(E, E )ν(E )Σ f (x, E )ψ(x, Ω , E )dΩ dE .
Therefore, equation (1.25) can be written as

T -1 F ψ(x, Ω, E) = λψ(x, Ω, E). As we considerate T -1 from L 2 (S 2 × R + , L 2 (R)) to L 2 (S 2 × R + , L 2 (R)) and thanks to Rellich theorem [AdFo03] H 1 (R) is compactly embedded in L 2 (R), T -1 is compact from L 2 (S 2 × R + , L 2 (R)) to L 2 (S 2 × R + , L 2 (R)) and thus T -1 F is compact too.
We apply the theorem of Krein-Rutman, theorem 1.3, to T -1 F , and we find that the only physical solution of equation (1.25) is associated to the greatest eigenvalue of the inverse operator T -1 F . We denote the criticality of the reactor core by k eff such that

k eff = max λ λ.
The criticality of a reactor core characterizes the physical state of the core:

• if k eff = 1: the reactor core is in a steady state and the nuclear chain reaction is self-sustaining. The reactor is said to be critical;

• if k eff > 1: there are more neutrons produced than neutrons lost. The chain reaction races. The reactor is said to be supercritical;

• if k eff < 1: there are less neutrons produced than neutrons lost. The chain reaction vanishes. The reactor is said to be subcritical.

Energy Discretization and Homogenization

The only physical data available are the pointwise microscopic cross sections. Therefore, in order to solve equation (1.25), one has to evaluate the macroscopic cross sections. These macroscopic cross sections are given as energy piecewise constants and spatial piecewise polynomial functions. The use of energy piecewise constant functions in neutronics is called the multi-group method. The spatial treatment is done by some homogenization techniques of the reactor geometry on a smaller scale [START_REF] Sanchez | Assembly Homogenization Techniques for Core Calculations[END_REF][START_REF] Coste-Delclaux | Modélisation Du Phénomène d'autoprotection Dans Le Code de Transport Multigroupe APOLLO2[END_REF]. We will not give details on how to obtain these values. First we assume that the neutron energy can not be smaller than E min > 0 and higher than E max > E min . Let G be an integer which represents the number of energy groups and

(E g ) G+1 g=1 such that E G+1 = E min , E 1 = E max , for all g < g, E g < E g and [E min , E max ] = G g=1 [E g+1 , E g ], see figure 1.1. Each energy subinterval [E g+1 , E g ] is called an energy group.
1.4 : Energy Discretization and Homogenization

0 E G+1 E G E g+1 E g E 2 E 1 E Group g Figure 1
.1: Representation of the energy groups.

Moreover, we approach the scattering macroscopic cross sections Σ t , Σ f , the fission yield ν and spectrum χ by constant value on each energy group.

∀1 ≤ g ≤ G, ∀E ∈ [E g+1 , E g ] :    Σ t (x, E) ≈ Σ g t (x); ν(E)Σ f (x, E) ≈ ν g Σ g f (x); χ(E, E ) ≈ χ g,g .
The macroscopic cross section depends only on the energy groups of the incident and scattered neutron energy:

∀ 1 ≤ g ≤ G, ∀ 1 ≤ g ≤ G, ∀E ∈ [E g +1 , E g ] : E g E g+1 Σ s (x, c, E → E)dE ≈ Σ g →g s (x, c).
Finally, we denote:

∀1 ≤ g ≤ G, ψ g (x, Ω) = E g E g+1 ψ(x, Ω, E)dE.
By integrating equation (1.25) over the enery group [E g+1 , E g ], g = 1, G, and using the trapezoidal rule,

+∞ 0 dE ≈ G g=1 E g E g+1
dE, for the scattering term, one can obtain the multi-group approximation of the steady state neutron transport equation which reads, for each energy group g ∈ {1, • • • , G}:

Ω • grad x ψ g (x, Ω)+Σ g t (x)ψ g (x, Ω) = G g =1 S 2 Σ g →g s (x, c)ψ g (x, Ω )dΩ + λ -1 1 4π G g =1 χ g,g ν g Σ g f (x) S 2 ψ g (x, Ω )dΩ .
(1.26)

This approximation corresponds to a P 0 approximation of the neutron flux over the energy.

Remark 1.4. As the fission spectrum χ corresponds to the normed energy distribution of the neutrons emitted by a fission, in the case of the one-group approximation, it stands χ 1,1 = 1.

One can easily obtain from equation (1.5) that on each energy group 1 ≤ g ≤ G, it stands:

Σ g t (x) = G g =1 S 2 Σ g→g s (x, c)dc + Σ g f (x) + Σ g (x)
.

(1.27)

Remark 1.5. Another way to model the energy dependence is to use the probability table [START_REF] Levitt | The Probability Table Method for Treating Unresolved Neutron Resonances in Monte Carlo Calculations[END_REF][START_REF] Mosca | Conception et Développement d'un Mailleur Énergétique Adaptatif Pour La Génération Des Bibliothèques Multigroupes Des Codes de Transport[END_REF]. This method consists in associating to every energy groups a discrete probability law of the value of the microscopic cross sections. This probability law is used to evaluate and homogenize the macroscopic cross sections.

Angular Discretization

To discretize the solution over the direction Ω, there exists several methods, let us mention the two most common, called the discrete ordinate method (S N ) and the spherical harmonic method (P N ). The first one consists in evaluating the solution over a set of quadrature points of the unit sphere and the second one is a projection of the solution and the macroscopic cross sections over the spherical harmonics. One can refer to [START_REF] Davidson | Neutron Transport Theory[END_REF] for both methods. There exists also the simplified P N method (SP N ) which can be derived from the P N method. In this work, we use the SP N method, which was described for the first time by Gelbard in [START_REF] Gelbard | Application of Spherical Harmonics Method to Reactor Problems[END_REF]. The approximation is done by projecting the solution and the macroscopic cross sections over the spherical harmonics and assuming that the flux is locally planar.

The Multigroup P N Transport Equations

Here, we describe how to obtain the multigroup P N transport equations. We recall from [START_REF] Hochstadt | The Functions of Mathematical Physics[END_REF], that any complex-valued function in L 2 (S 2 ) can be expanded in terms of the spherical harmonics [START_REF] Hochstadt | The Functions of Mathematical Physics[END_REF], Y m n for n ∈ N and m ∈ Z such that |m| ≤ n, therefore, for any energy groups g in {1, • • • , G}, the angular flux can be denoted by:

ψ g (x, Ω) = ∞ n=0 n m=-n ψ g n,m (x)Y m n (Ω),
where

ψ g n,m (x) = S 2 ψ g (x, Ω)Y m n (Ω)dΩ ∈ C.
One can remark that the fission integral in equation (1.26), can be simplified as

1 4π G g =1 χ g,g ν g Σ g f (x) S 2 ψ g (x, Ω )dΩ = 1 4π G g =1 χ g,g ν g Σ g f (x)ψ g 0,0 (x).
The scattering macroscopic cross sections depend on c = cos ΩΩ , where ΩΩ is the angle between the incident direction and the new one. In order to evaluate the scattering source integral, we want to expand the scattering macroscopic cross sections over the spherical harmonics. We recall from [START_REF] Hochstadt | The Functions of Mathematical Physics[END_REF], that any function in L 2 (] -1; 1[) can be expanded on Legendre polynomials, P n for n ∈ N, therefore, for any energy groups g and g in {1, • • • , G}, Σ g→g s can be rewritten as

Σ g→g s (x, c) = ∞ n=0 2n + 1 4π Σ g→g s,n (x)P n (c), (1.28) where Σ g→g s,n (x) = 2π 1 -1 Σ g→g s (x, c)P n (c)dc.
More particularly, the first moment of the scattering macroscopic cross section from the energy group g ∈ {1, • • • , G} to the energy group g ∈ {1, • • • , G} is given by:

Σ g→g s,1 = 2π 1 -1 Σ g→g s (c)P 1 (c)dc.
Using the addition theorem B.4, one obtains the following decomposition of the scattering macroscopic cross sections:

Σ g→g s (x, c) = ∞ n=0 Σ g→g s,n (x) n m=-n Y m n (Ω)Y m n (Ω ).
(1.29) Therefore, the scattering source integral can be rewritten as

S 2 Σ g →g s (x, c)ψ g (x, Ω )dΩ = ∞ n=0 n m=-n Σ g →g s,n (x)ψ g n,m (x)Y m n (Ω).
Using this expression in equation (1.26) and we integrate over the unit sphere against the conjugate of the spherical harmonic Y m n : We recall that Ω can be described with its spherical coordinates (ϑ, θ), see figure 1

∞ n =0 n m =-n grad x ψ g n ,m (x) • S 2 ΩY m n (Ω)Y m n (Ω)dΩ + Σ g t (x)ψ g n,m (x) = G g =1 Σ g →g s,n (x)ψ g n,m (x) + λ -1 G g =1 χ g,g ν g Σ g f (x)δ n,0 δ m,0 ψ g 0,0 (x).
                             Ω x Y m n (ϑ, θ) = 1 2 b x,y (n, m)Y m + n -(ϑ, θ) -b x,y (n + , -m + )Y m + n + (ϑ, θ) +b x,y (n + , m -)Y m - n + (ϑ, θ) -b x,y (n, -m)Y m - n -(ϑ, θ) ; Ω y Y m n (ϑ, θ) = i 2 b x,y (n, m)Y m + n -(ϑ, θ) -b x,y (n + , -m + )Y m + n + (ϑ, θ) -b x,y (n + , m -)Y m - n + (ϑ, θ) + b x,y (n, -m)Y m - n -(ϑ, θ) ; Ω z Y m n (ϑ, θ) = b z (n + , m)Y m n + (ϑ, θ) + b z (n, m)Y m n -(ϑ, θ). Where b z (n, m) = (n + m)(n -m) (2n + 1)(2n -1) and b x,y (n, m) = (n -m -1)(n -m) (2n + 1)(2n -1) .
In the previous equations, some terms may appear without being defined, for instance when n = 0 in the third equation Y m -1 appears, in this case they are null, for more explanation the reader can refer to appendix B. To obtain the P N transport equations, one can use the previous relation to integrate the integral term in (1.26) and obtain:

b x,y (n + , -m + ) 2 (∂ x + i∂ y ) ψ g n + ,m + (x) + b x,y (n + , m -) 2 (∂ x -i∂ y ) ψ g n + ,m -(x) + b x,y (n, m) 2 (∂ x + i∂ y ) ψ g n -,m + (x) + b x,y (n, -m) 2 (-∂ x + i∂ y ) ψ g n -,m -(x) +b z (n, m)∂ z ψ g n -,m (x) + b z (n + , m)∂ z ψ g n + ,m (x) + Σ g t (x)ψ g n,m (x) = G g =1 Σ g →g s,n (x)ψ g n,m (x) + λ -1 G g =1 χ g,g ν g Σ g f (x)δ n,0 δ m,0 ψ g 0,0 (x).

The Multigroup SP N Transport Equations

The P N transport equations can be simplified with the SP N method. The transport is decomposed with three main directions as shown in figure 1.3. We suppose that the transport is mainly in one direction, e z , thus, the flux ψ can be described as a function of the axial coordinate z and of the angle between Ω and e z . Under this hypothesis, we denote = Ω • e z and ψ g,z (z, ) := ψ g (x, Ω). Moreover, this hypothesis implies that the flux gradient is collinear with e z . Then the neutron transport equation (1.26) becomes:

∂ z ψ g,z (z, ) + Σ g t (x)ψ g,z (z, ) = G g =1 1 -1 Σ g →g s (x, c)ψ g ,z (z, )d +λ -1 1 4π G g =1 χ g,g ν g Σ g f (x) 1 -1 ψ g ,z (z, )d .
(1.31)

As in the P N transport derivation, we project the scattering macroscopic cross sections on Legendre polynomials (1.28). The angular flux is also projected on this polynomial family but we introduce some weight, (α n ) ∞ n=0 , that we specify later in order to simplify the computation ultimately:

ψ g,z (z, ) = ∞ n=0 α n ψ g,z n (z)P n ( ),
where

ψ g,z n (z) = 2n + 1 2α n 1 -1 ψ g,z (z, )P n ( )d .
Moreover, one can rewrite (1.29) with the help of the addition theorem (B.5) as:

Σ g→g s (x, c) = ∞ n=0 2n + 1 4π Σ g→g s,n (x) P n ( )P n ( ) +2 n m=1 (n -m)! (n + m)! P m n ( )P m n ( ) cos(m(ϑ -ϑ )) ,
where ϑ (resp. ϑ ) is given by the spherical coordinates of Ω (resp. Ω ).

As we suppose that the angular flux does not depend on the angle ϑ, we remark that the term 2π 0 cos(m(ϑϑ ))ψ g ,z (z, )dϑ is null. Thus, it stands

2π 0 1 -1 Σ g →g s (x, c)ψ g ,z (z, )d dϑ = ∞ n=0 α n Σ g→g s,n (x)ψ g ,z n (z)P n ( ).
By injecting those two projections in (1.31), multiplying it by αn 4π P n and integrating it over S 2 , for n ∈ N, one obtains the following equation:

∞ n =0 α n α n 2 ∂ z ψ g,z n (z) 1 -1 P n ( )P n ( )d +t n Σ g t (x)ψ g,z n (z) = G g =1 t n Σ g→g s,n (x)ψ g ,z n (z) + λ -1 G g =1 χ g,g ν g t 0 Σ g f (x)δ n,0 ψ g ,z 0 (z),
where we denote

t n = α 2 n 2n + 1 .
The integral term in the previous equation can be evaluated thanks to the recursive relations (B.1). Then, one can obtain:

h n+1 ∂ z ψ g,z n+1 (z) + h n ∂ z ψ g,z n-1 (z) + t n Σ g t (x)ψ g,z n (z) = G g =1 t n Σ g→g s,n (x)ψ g ,z n (z) + λ -1 G g =1 χ g,g ν g t 0 Σ g f (x)δ n,0 ψ g ,z 0 (z),
where we denote

h n = nα n α n-1 (2n + 1)(2n -1)
for n ≥ 1.

It is customary to choose the family (α n ) ∞ n=0 such that, for n ≥ 1, h n = 1. We choose here:

   α 0 = 1;

α n = 4n 2 -1 nα n-1 for n ≥ 1.
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • t n • α n Figure 1.4: Values of (α n ) ∞ n=0 and (t n ) ∞ n=0 .
With the same methodology one can use e x and e y as transport main axis and obtain the same equations for the moments of ψ g,x (x, Ω • e x ) := ψ g (x, Ω) and ψ g,y (y, Ω • e y ) := ψ g (x, Ω). For any d ∈ {x, y, z} it stands:

h n+1 ∂ d ψ g,d n+1 (x d ) + h n ∂ d ψ g,d n-1 (x d ) + t n Σ g t (x)ψ g,d n (x d ) = G g =1 t n Σ g→g s,n (x)ψ g ,d n (x d ) + λ -1 G g =1 χ g,g ν g t 0 Σ g f (x)δ n,0 ψ g ,d 0 (x d ).
(1.32)

As we decompose the transport directions in three directions and we neglect the couplings between these three directions, the SP N approximation solution does not converge towards the transport equation solution when N tends to infinity, whereas the P N approximation solution does. We denote φ g 0 (x) = d∈{x,y,z}

ψ g,d 0 (x d ).
Then, for n = 0, we sum equations (1.32) over d ∈ {x, y, z}:

h 1 d∈{x,y,z} ∂ d ψ g,d 1 (x d ) + t 0 Σ g t (x)φ g 0 (x) = G g =1 t 0 Σ g→g s,0 (x)φ g 0 (x) + λ -1 G g =1 χ g,g ν g t 0 Σ g f (x)φ g 0 (x).
The derivative sum is like a divergence, thus we denote p g 1 (x) = d∈{x,y,z} ψ g,d 1 (x d )e d and we obtain:

h 1 div p g 1 (x) + t 0 Σ g t (x)φ g 0 (x) = G g =1 t 0 Σ g→g s,0 (x)φ g 0 (x) + λ -1 G g =1 χ g,g ν g t 0 Σ g f (x)φ g 0 (x).
For n = 1, one can obtain by summing equation (1.32)e d over d ∈ {x, y, z}:

h 2 grad x φ g 2 (x) + h 1 grad x φ g 0 (x) + t 1 Σ g t (x)p g 1 (x) = G g =1 t 1 Σ g→g s,1 (x)p g 1 (x),
where φ g 2 (x) = d∈{x,y,z}

ψ g,d 2 (x d ).
One can continue to higher n value and define odd moments as vectors:

p g 2n+1 (x) = d∈{x,y,z} ψ g,d 2n+1 (x d )e d ,
and even moments as scalars:

φ g 2n (x) = d∈{x,y,z} ψ g,d 2n (x d ).
Let N ∈ N odd denotes the order of the SP N method. From now on, we denote the vector of the odd moment of the flux by p = (p

1 , • • • , p G ) T with p g = (p g 1 , p g 3 , • • • , p g N )
T and the vector of the even moment of the flux φ

= (φ 1 , • • • , φ G ) T with φ g = (φ g 0 , φ g 2 , • • • , φ g N -1
) T . We denote also N = N -1 2 . We truncate the expansion of the angular flux to the N th Legendre polynomial. Therefore, the multigroup SP N transport equations read as coupled diffusion-like equations set in a mixed formulation:

T o p + grad x Hφ = 0; H T div p + T e φ = λ -1 M f φ.
(1.33)

In the following, we describe the operators T o , T e , M f and H used in equation (1.33). The operators T o and T e are defined by energy blocks:

T o =      T 1 o -S 2→1 o • • • -S G→1 o -S 1→2 o . . . . . . . . . . . . . . . . . . -S G→G-1 o -S 1→G o • • • -S G-1→G o T G o      , (1.34)
and In order to define the blocks of these two operators, we introduce the moments of the removal macroscopic cross sections:

T e =      T 1 e -S 2→1
Σ g r,n (x) = Σ g t (x) -Σ g→g s,n (x) (1.36)
The blocks of T o and T e are given by:

T g o = diag t 2n+1 Σ g r,2n+1 N n=0 ; S g →g o = diag t 2n+1 Σ g →g s,2n+1 N n=0 ; T g e = diag t 2n Σ g r,2n N n=0 ; S g →g e = diag t 2n Σ g →g s,2n N n=0
.

The fission operator reads

M f =    M 1→1 f • • • M G→1 f . . . . . . . . . M 1→G f • • • M G→G f    ,
where the blocks are

M g →g f = χ g,g t 0 ν g Σ g f δ 2m,0 δ 2n,0 N n,m=0
.

Finally the moment coupling matrix is defined as:

H = diag ( Ĥ) G g=1 ,
where

Ĥ =         1 1 0 • • • 0 0 . . . .
• • • • • • 0 1        
.

The SP N approximation has less angular moments than the P N approximation. For a problem with d spatial dimensions and a given N ∈ N odd, the SP N approximation has (d + 1) N +1 2 moments whereas the P N approximation has (N + 1) 2 moments. Moreover, as we suppose that there is no coupling between the transport directions, the SP N linear system is sparser than the P N one. Indeed, for a problem with d spatial dimensions and G groups and a given N ∈ N odd, the SP N system has G 2 N +1 2 + GN non null coefficients and the P N has G 2 (N + 1) 2 + 6GN 2 -1 non null coefficients. Thus, the advantage of the SP N method over the P N approximation is that the first one has less angular moment coupling than the second, see table 1.1. On the other hand, due to the hypothesis of no coupling between the transport directions, the SP N solution does not converge towards the transport solution but the P N solution does.

Spherical Harmonic Moments of the Cross Sections

In this work we are interested in the simulation of pressurized water reactor (PWR). As a macroscopic section can not be infinite, they are supposed to be bounded over R. The space of the bounded functions over R is represented by L ∞ (R). Moreover, the macroscopic cross sections are homogenized on a smaller scale than the reactor geometry, so we supposed that the removal macroscopic cross sections are piecewise regular. We introduce the function space:

W 1,∞ (R) = {f ∈ L ∞ (R), grad x f ∈ L ∞ (R)} .
Let N h be an integer and (R i ) N h i=1 be the partition of R on which the macroscopic cross sections are homogenized, i.e the homogenized macroscopic cross sections are regular on each R i , for 1 ≤ i ≤ N h . Thus, for all 1 ≤ g, g ≤ G and for all 1 ≤ n ≤ N , we suppose that:

(Σ g r,n , Σ g →g s,n , νΣ g f ) ∈ PW 1,∞ (R) × L ∞ (R) × L ∞ (R), (1.37)
where

PW 1,∞ (R) = D ∈ L ∞ (R), D |R i ∈ W 1,∞ (R i ), 1 ≤ i ≤ N h .
In PWR reactor, the scattering macroscopic cross sections are nearly isotropic, that means that the scattering macroscopic cross sections admit only small variations over the c variables. Therefore the Legendre moments of a scattering macroscopic cross section of order n, with 1 ≤ n ≤ N , are smaller than the Legendre moment of order 0. For all 1 ≤ g, g ≤ G, for all 0 ≤ n ≤ N , it stands:

|Σ g→g s,n | ≤ Σ g→g s,0
a.e. in R.

(1.38) Furthermore, from equation (1.27), for all 1 ≤ g ≤ G, it holds:

Σ g t = G g =1 Σ g→g s,0 + Σ g f + Σ g a.e. in R.
(1.39)

We recall that, from inequality (1.4), for all 1 ≤ g ≤ G, Σ g > 0. According to (1.36) and (1.39), for all 1 ≤ g, g ≤ G and for all 0 ≤ n ≤ N , it exists (Σ r,n ) * > 0 such that 0 < (Σ r,n ) * ≤ t n Σ g r,n . Moreover, the macroscopic cross sections are bounded, thus, for all 1 ≤ g, g ≤ G and for all 0 ≤ n ≤ N , it exists (Σ r,n ) * > 0 such that t n Σ g r,n ≤ (Σ r,n ) * . To summarize, we have the following hypothesis on the removal macroscopic cross sections, for all 1 ≤ g ≤ G and for all 0 ≤ n ≤ N :

∃(Σ r,n ) * , (Σ r,n ) * > 0, 0 < (Σ r,n ) * ≤ t n Σ g r,n ≤ (Σ r,n ) * a.e. in R. (1.40)
By combining the definition of the removal macroscopic cross section in (1.36), equation (1.39) and inequality (1.38), one obtains for all 1 ≤ g ≤ G and for all 0 ≤ n ≤ N :

Σ g r,n > g =g Σ g→g s,n
a.e. in R.

More precisely, when modelling a PWR, the scattering macroscopic cross sections are such that for all 1 ≤ g, g ≤ G, for all 0 ≤ n ≤ N :

∃0 < ε < 1 G -1 , |Σ g→g s,n | ≤ εΣ g r,n
a.e. in R.

(1.41)

Assumption (1.40) is less restrictive than assumption (1.41), but both are satisfied for PWR. Finally, for each energy group 1 ≤ g ≤ G, we denote by νΣ g f the product ν g Σ g f which is positive:

0 ≤ νΣ g f a.e. in R, (1.42)
and the matrix M f is non null, there exists g and g such that χ g,g νΣ g f = 0. We regroup the assumptions in (1.37),(1.40),(1.41) and (1.42) to obtain the following physical hypotheses on the coefficients: Hypothesis 1.6. For all energy groups 1 ≤ g, g ≤ G, g = g and for all 0 ≤ n ≤ N ,it stands:

               (Σ g r,n , Σ g →g s,n , νΣ g f ) ∈ PW 1,∞ (R) × L ∞ (R) × L ∞ (R), ∃(Σ r,n ) * , (Σ r,n ) * > 0, 0 < (Σ r,n ) * ≤ t n Σ g r,n ≤ (Σ r,n ) * a.e. in R, ∃0 < ε < 1 G -1 , |Σ g→g s,n | ≤ εΣ g r,n a.e. in R 0 ≤ νΣ g f a.
e. in R, ∃g, g s.t. χ g,g νΣ g f = 0.

(1.43a)

(1.43b) (1.43c) (1.43d)
As in this work, we suppose that the coefficients are valid for PWR simulation, hypothesis 1.6 is always satisfied. Thanks to assumption (1.43c) in condition 1.6, the transposed matrices of T o and T e are diagonally dominant and thus T o and T e are invertible. Equations (1.33) can be read as a set of coupled diffusion-like equations written in a primal formulation:

-T Hdiv (T -1 o grad x (Hφ)) + T e φ = λ -1 M f φ. (1.44)
In this study, we need a positive property on the matrices T o and T -1 e for the mixed formulation (1.33) or on the matrices T -1 o and T e for the primal formulation (1.44). The methodology used to have these properties is first to exhibit a condition on the macroscopic cross sections which ensures that T h , h ∈ {e, o}, has a positive property, then we show another condition on the macroscopic cross sections such that T -1 h , h ∈ {e, o}, has a positive property. To do that, we need first to introduce the following operators:

• :      R d × R d → R (x, y) → x • y = d i=1 x i y i • :        (R N ) G × (R N ) G → R (x, y) → x • y = G g=1 N n=1 x g n y g n :        ((R d ) N ) G × ((R d ) N ) G → R (x, y) → x y = G g=1 N n=1 x g n • y g n We define J o = {2h + 1 | h = 0, N } and J e = {2h |, h = 0, N }. Let h be in {e, o}, we also define T h,r = diag (T g h ) G g=1 and U h = (U g,g h ) G g,g =1
where U g,g h = 0 when g = g and

U g,g h = 1 ε (T g h,r ) -1 S g →g h = diag (Σ g →g s,n /εΣ g r,n
) n∈J h when g = g . Thus the matrix T h can be decomposed into:

T h = T h,r (I -εU h ).
(1.45) Moreover, we define the extrema of the removal macroscopic cross section.

Definition 1.7. For h ∈ {e, o}, we denote the supremum of the removal macroscopic cross section by:

(Σ h r ) * = max n∈J h (Σ r,n ) * ,
and the infimum of the removal macroscopic cross section:

(Σ h r ) * = min n∈J h (Σ r,n ) * .
Definition 1.8. For h ∈ {e, o}, we define:

α s,h = max n∈J h max g,g =g sup x∈R |Σ g →g s,n | Σ g r,n
, and:

α r,h = (Σ h r ) * (Σ h r ) * .
The following proposition gives us a bound on the norm of U h x, for all x in (R N ) G . This bound is used to obtain the conditions on the macroscopic cross sections such that T h and T -1 h have a positive property.

Proposition 1.9. For h ∈ {e, o}, it stands, for all x in (R N ) G :

U h x ≤ α s,h ε (G -1) x a.e. in R.
(1.46)

Proof. Let x be in (R N ) G , then it stands, for all 1 ≤ g ≤ G and for all n in J h :

(U h x) g n = g =g Σ g →g s,n εΣ g r,n
x g n .

As it stands U h x 2 = U h x • U h x, from the definition of •, one obtains:

U h x 2 = G g=1 n∈J h g =g Σ g →g s,n εΣ g r,n
x g n 2 .

We recall that for any real sequence (a n ) n∈N and A ∈ N * , it holds

A n=1 a n 2 ≤ A A n=1
a 2 n . Then, using this result and the definition of α s,h , it stands:

U h x 2 ≤ α 2 s,h ε 2 (G -1) G g=1 n∈J h g =g |x g n | 2 a.e. in R.
In the previous inequality, we rewrite

g =g |x g n | 2 as G g =1 |x g n | 2 -|x g n | 2
, and we obtain:

U h x 2 ≤ α 2 s,h ε 2 (G -1) G g=1 n∈J h G g =1 |x g n | 2 -|x g n | 2 a.e. in R.
Using the definition of the norm on (R N ) G , we find:

U h x 2 ≤ α 2 s,h ε 2 (G -1) G g=1
x 2x 2 a.e. in R.

Finally, it stands:

U h x 2 ≤ α 2 s,h
ε 2 (G -1) 2 x 2 a.e. in R.

Under the condition 1.6, the matrix T h,r , h ∈ {e, o}, is diagonal with strictly positive coefficient, thus it is invertible and we have the following bounds, for all x in (R N ) G :

T h,r x • x ≥ (Σ h r ) * x 2 a.e. in
α s,h α r,h (G -1) < 1, (1.49)
then for all x ∈ (R N ) G it stands:

T h x • x ≥ (T h ) * x 2 a.e. in R,
where

(T h ) * = (Σ h r ) * (1 -α s,h α r,h (G -1)) .
(1.50)

Proof. Let x be in (R N ) G . Then, from equation (1.45), it holds:

T h x • x = T h,r x • x -ε(T h,r U h x • x).
(1.51)

As T h,r is diagonal, one obtains:

-ε(T h,r U h x • x) = -ε(U h x • T h,r x).
Thanks to Cauchy-Schwarz inequality, (1.48) and (1.46), one obtains:

-ε(U h x • T h,r x) ≥ -α s,h (G -1)(Σ h r ) *
x 2 a.e. in R. Using this inequality and inequality (1.47) in equation (1.51), we obtain:

T h x • x ≥ (Σ h r ) * (1 -α s,h α r,h (G -1))
x 2 a.e. in R. Thanks to (1.49) we have (T h ) * > 0.

Under the hypothesis 1.6, the matrix T h,r , h ∈ {e, o}, is invertible and we have the following bound, for all x in (R N ) G :

T -1 r,h x ≤ 1 (Σ h r ) *
x a.e. in R.

(1.52)

Proposition 1.11. Let h ∈ {o, e}. Under hypothesis 1.6 and if α s,h and α r,h satisfy:

α s,h (G -1) < 1 1 + α r,h , (1.53)
then for all x ∈ (R N ) G it stands:

T -1 h x • x ≥ (T -1 h ) * x 2 a.e. in R, where (T -1 h ) * = 1 (Σ h r ) * 1 -α r,h α s,h (G -1) 1 -α s,h (G -1) (1.54)
Proof. As the matrices T h and T h,r are invertible, the matrix (I -εU h ) is invertible too, and it stands:

T -1 h = (I -εU h ) -1 T -1 h,r , (1.55)
where ε is defined in hypothesis 1.6. From proposition 1.9 and assumption (1.43c), (I -εU h ) -1 can be rewritten as:

(I -εU h ) -1 = (I + l>0 ε l U l h ).
Substituting this expansion in equation (1.55), one obtains:

T -1 h = (I + l>0 ε l U l h )T -1 h,r .
(1.56) Let x be in (R N ) G , we look for a lower bound to:

T -1 h x • x = T -1 h,r x • x + l>0 ε l U l h T -1 h,r x • x.
First we bound each term in the sum over l. From Cauchy-Schwartz inequality, we obtain for all l > 0:

U l h T -1 h,r x • x ≥ -U h T -1 h,r x x .
Using equations (1.46) and (1.52), one obtains:

U l h T -1 h,r x • x ≥ - α l s,h ε l (G -1) l 1 (Σ h r ) * x 2 a.e. in R.
Furthermore, it stands:

T -1 h,r x • x = G g=1 n∈J h 1 t n Σ g r,n |x g h | 2 ; ≥ 1 (Σ h r ) * x 2 a.e. in R.

Finally we obtain:

T -1 h x • x ≥ 1 (Σ h r ) * - 1 (Σ h r ) * l>0 α l s,h (G -1) l x 2 a.e. in R; ≥ 1 (Σ h r ) * - 1 (Σ h r ) * α s,h (G -1) 1 -α s,h (G -1) x 2 a.e. in R; ≥ 1 (Σ h r ) * 1 -α r,h α s,h (G -1) 1 -α s,h (G -1)
x 2 a.e. in R.

Thanks to inequality (1.53) we have (T -1 h ) * > 0. At last, we bound the matrices T h and T -1 h uniformly.

Proposition 1.12. Let h be in {e, o}. For all x in (R N ) G , it stands:

T h x ≤ (T h ) * x a.e. in R, (1.57)
where

(T h ) * = (Σ h r ) * (1 + α s,h (G -1)) ; (1.58)
and

T -1 h x ≤ (T -1 h ) * x a.e. in R, (1.59) where (T -1 h ) * = 1 (Σ h r ) * (1 -α s,h (G -1))
.

(1.60)

Proof. Let x be in (R N ) G . Using the decomposition of T h in (1.45) and the triangle inequality, one obtains:

T h x ≤ T r,h x + ε T h,r U h x .
Thanks to inequalities (1.48) and (1.46), we obtain:

T h x ≤ (Σ h r ) * (1 + α s,h (G -1)) x a.e. in R.
Now, we prove inequality (1.59). Thanks to equation (1.56) and inequality (1.52), one obtains:

T -1 h x ≤ 1 (Σ h r ) * I + l>0 ε l U l h x .
Thanks to proposition 1.9, we know that the series ε l U h x l converges, we can use the triangle inequality on x + l>0 ε l U l h x . Thus, we obtain:

T -1 h x ≤ 1 (Σ h r ) * x + l>0 ε l U l h x .
From inequality (1.46) and the fact that

1 + l>0 α l s,h (G -1) l = (1 -α s,h (G -1)) -1 , one obtains: T -1 h x ≤ 1 (Σ h r ) * (1 -α s,h (G -1))
x .

For the mixed formulation of the multigroup SP N equations we need a positive property on T o and T -1 e , thus for this formulation we make the following assumption:

Condition 1.13. For the mixed formulation of the multigroup SP N transport equations, in addition to hypothesis 1.6, we suppose that:

   α s,o α r,o (G -1) < 1; α s,e (G -1) < 1 1 + α r,e .
Therefore, for the mixed formulation the matrices T o and T -1 e are bounded and have a positive property.

Proposition 1.14. Under condition 1.13, it exists (T o ) * > 0 and (T -1 e ) * > 0 such that for all x ∈ (R N ) G it holds:

T o x • x ≥ (T o ) * x 2 ; (1.61) and T -1 e x • x ≥ (T -1 e ) * x 2 . (1.62) Moreover, it exists (T o ) * > 0 and (T -1 e ) * > 0 such that for all x ∈ (R N ) G it holds: T o x ≥ (T o ) * x ; (1.63) and T -1 e x ≥ (T -1 e ) * x .
(1.64)

For the primal formulation (1.44), we need a positive property on T -1 o and T e , thus for this formulation we make the following assumption: Condition 1.15. For the primal formulation of the multigroup SP N transport equations, in addition to hypothesis 1.6, we suppose that:

   α s,o (G -1) < 1 1 + α r,o ; α s,e α r,e (G -1) < 1.
Therefore, for the primal formulation the matrices T -1 o and T e are bounded and have a positive property.

Proposition 1.16. Under condition 1.15, it exists (T -1 o ) * > 0 and (T e ) * > 0 such that for all x ∈ (R N ) G it holds:

T -1 o x • x ≥ (T -1 o ) * x 2 .
(1.65)

and

T e x • x ≥ (T e ) * x 2 ;
(1.66)

Moreover, it exists (T e ) * > 0 and (T -1 o ) * > 0 such that for all x ∈ (R N ) G it holds: T -1 o x ≥ (T -1 o ) * x .
(1.67)

and

T e x ≥ (T e ) * x ; (1.68)
In general, in PWR simulation 2 ≤ G 36 and N ∈ {1, 3, 5}.

Special Instances of these Equations

In this subsection, we explicit the multigroup SP 1 and SP 3 neutron transport equations, and the one-group neutron diffusion equations. We will solve these models mathematically and numerically in the next chapters.

The Multigroup SP 1 Neutron Transport Equations

If the angular flux is considered to be independent from the angular direction, we only need to truncate the Legendre polynomial expansion of the angular flux and the macroscopic cross sections to the first term. We first need to calculate the coefficients t 0 and t 1 , for that the values of α 0 and α 1 are required:

α 0 = 1; α 1 = 3; t 0 = 1; t 1 = 3.
Moreover, we define the diffusion coefficient, for all energy group g ∈ {1, • • • , G}:

D g = 3Σ g r,1 .
With this truncation, one obtains the multigroup SP 1 neutron transport equations: For all g ∈ {1, • • • , G}, solve in (p g , φ g , λ) such that:

                 D g p g 1 (x) + grad x φ g 0 (x) = g =g 3Σ g →g s,1 (x)p g 1 (x); div p g 1 (x) + Σ g r,0 φ g 0 (x) = g =g Σ g →g s,0 (x)φ g 0 (x) +λ -1 G g =1 χ g,g νΣ g f (x)φ g 0 (x).
(1.69)

The Multigroup SP 3 Neutron Transport Equations

When the Legendre polynomial expansion of the angular flux and the cross sections are truncated to three terms, one obtain the multigroup SP 3 neutron transport equations: For all g ∈ {1, • • • , G}, solve in (p g , φ g , λ) such that:

                                     div p g 1 (x) + Σ g r,0 φ g 0 (x) = g =g Σ g →g s,0 (x)φ g 0 (x) +λ -1 G g =1 χ g,g νΣ g f (x)φ g 0 (x); 3Σ g r,1 p g 1 (x) + grad x φ g 0 (x) + grad x φ g 2 (x) = g =g 3Σ g →g s,1 (x)p g 1 (x); div p g 1 (x) + div p g 3 (x) + 5 2 Σ g r,2 φ g 2 (x) = g =g 5 2 Σ g →g s,2 (x)φ g 2 (x) 14 3 Σ g r,3 p g 3 (x) + grad x φ g 2 (x) = g =g 14 3 Σ g →g s,3 (x)p g 3 (x).
(1.70)

The Neutron Diffusion Equations

Here we only take one energy group (G = 1), so with an abuse of notation we omit the energy superscript, we denote the scalar flux by

φ(x) = S 2 ψ(x, Ω)dΩ,
and we denote the scalar current by

p(x) = S 2
Ωψ(x, Ω)dΩ.

With this notation, one can integrate equation (1.26) over S 2 and obtain:

div p(x) + Σ r,0 (x)φ(x) = λ -1 νΣ f (x)φ(x) (1.71)
This is the first equation of the neutron diffusion equations. The second equation is a Fick's law which reads:

1 D(x) p(x) + grad x φ(x) = 0, (1.72)
where D(x) is the diffusion coefficient. Using the physical assumptions that the medium is homogeneous, the variations of the flux are weak, and the spatial dependence of the flux is linear, computing the neutrons flux across a small surface gives

D(x) = 1 3Σ t (x)
.

Thus the neutron diffusion equations reads:

   1 D(x) p(x) + grad x φ(x) = 0; div p(x) + Σ t (x)φ(x) = λ -1 νΣ f (x)φ(x).
(1.73)

Boundary Conditions

In order to close the multigroup SP N transport equations (1.33), one need to add some boundary conditions. There is three main boundary conditions:

• The flux is null on the boundary (Dirichlet boundary conditions):

φ = 0 on ∂R;
(1.74)

• The neutrons are reflected by the boundary (Neuman boundary conditions):

p • n = 0 on ∂R;
(1.75)

• Albedo boundary conditions:

φ + c a p • n = 0 on ∂R;
(1.76) with c a strictly positive.

In the case of the domain R can be written as

R = d d=1 ]0; L d [, where L d > 0, 1 ≤ d ≤ d,
we can define periodic boundary conditions:

∀1 ≤ d ≤ d, ∀x ∈ d-1 d =1 ]0; L d [{0} × d d =d+1 ]0; L d [, φ(x) = φ(x + L d e d ) p(x) = p(x + L d e d ) .
(1.77)

Part I Continuous Problem

In this part we study the well-posedness of multigroup SP N equations under their primal formulation and under their mixed formulation. First, we consider the one-group diffusion model.

The diffusion equations are studied in [START_REF] Crank | The Mathematics of Diffusion[END_REF][START_REF] Ern | Theory and Pratice of Finite Elements[END_REF] for the primal setting and in [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF] for the mixed setting. Moreover, the study of the eigenvalue diffusion problem is done thanks to the spectral linear compact operator theory [START_REF] Dunford | Linear Operators[END_REF][START_REF] Brezis | Analyse Fonctionelle -Théorie et Applications[END_REF]. This part is organized as follow:

• in chapter 2, we recall the results for the diffusion problem in both primal and mixed settings;

• in chapter 3, we extend them to the multigroup SP N transport problem.

Chapter 2

Neutron Diffusion Equations

In order to do the study of the neutron diffusion equations, we need to introduce some function spaces. Let R be an bounded, connected and open subset of R d (d = 1, 2, 3), having a Lipschitz boundary which is piecewise smooth. First, we denote L 2 (R) the Hilbert space of all the square-integrable functions over R. The inner product (resp. the norm

) of L 2 (R) is denoted by (•, •) 0,R (resp. • 0,R
) and is defined by, for all f and g in L 2 (R):

(f, g) 0,R = R f g.
For any s ∈ R, the Sobolev space of order s over R is denoted by H s (R) and it is a Hilbert space with the inner product (resp. norm) (•, •) s,R (resp. • s,R ). The broken Sobolev spaces are defined by:

PH s (R) = ψ ∈ L 2 (R), ψ |R i ∈ H s (R i ), 1 ≤ i ≤ N h , s > 0.
We recall that (R i ) N h i=1 is the homogenization partition defined in § 1.5.3. If s is a non-negative integer, the inner product and the norm of H s (R) are defined by, for all f and g in H s (R):

(f, g) s,R = |α|≤s R D α f D α g
If s is a positive real and is not an integer, then the inner product of H s (R) are defined by, for all f and g in H s (R):

(f, g) s,R = (f, g) [s],R + |α|≤s R R (D α f (x) -D α f (y))(D α g(x) -D α g(y)) |x -y| d+2(s-[s]) dxdy,
where [s] is the integer part of s.

The set of all functions in C ∞ (R) with compact support is denoted by D(R). For any s > 0, the closure of

D(R) in H s (R) is denoted by H s 0 (R).
We define H -s (R) as the dual space of H s 0 (R), for any s > 0 and the dual product is denoted by

•, • H -s (R),H s 0 (R) .
The set of all the continuous linear applications from E into E, E a Hilbert space, is denoted by L(E) and • L(E) represents its norm. Finally, we denote

V = H 1 0 (R), L = L 2 (R) and V = H -1 (R).
From the several models presented in § 1, we are focusing here on the stationary one-group (G = 1) SP 1 /diffusion equations (1.73). Equations (1.73) are written under a mixed form. We first write it under its primal form. From the first equation of (1.73), one can deduce that:

p(x) = D(x)grad x φ(x).
By substituting the value of p in the second equation of (1.73), one obtains the primal form of equation (1.73):

-div (Dgrad x φ) + Σ r,0 φ = 1 λ νΣ f φ. (2.1)
We recall that the difference between the diffusion problem and the SP 1 is the value of the coefficient D:

• Diffusion: D(x) = 3Σ t (x); • SP 1 : D(x) = 3Σ r,1 (x).
In order to close equation (2.1), we add some boundary conditions. We choose here to impose the flux to be null on the boundary of the reactor R.

Primal Approach for the Diffusion Equations

In this section we study the diffusion equation under its primal form, the problem is written with only the scalar flux and the current does not appear. The SP 1 /diffusion problem written under its primal form is:

Problem 2.1. Find (λ, φ) ∈ R × V \{0} such that: -div (Dgrad x φ) + Σ r,0 φ = λ -1 νΣ f φ in R. (2.2)
Remark 2.2. In the case of other boundary conditions, the following work can be applied with small modifications. Before studying the eigenvalue problem 2.1, we consider the associated source problem given by: Problem 2.3. For a given source S f in V , find φ ∈ V such that:

-div (Dgrad x φ) + Σ r,0 φ = S f in R. (2.
3)

The physical hypothesis on the coefficients given in hypothesis 1.6 for the multigroup SP N transport equations reads in our case (G = 1, N = 1):

Hypothesis 2.4.                    (D, Σ r,0 , νΣ f ) ∈ PW 1,∞ (R) × PW 1,∞ sym (R) × L ∞ (R), ∃(D) * , (D) * > 0, 0 < (D) * ≤ D ≤ (D) * a.e. in R, ∃(Σ r,0 ) * , (Σ r,0 ) * > 0, 0 < (Σ r,0 ) * ≤ Σ r,0 ≤ (Σ r,0 ) * a.e. in R, 0 ≤ νΣ f a.e. in R, νΣ f = 0.

: Primal Approach for the Diffusion Equations

Under hypothesis 2.4, the solution φ has some extra regularity (see [START_REF] Costabel | Singularities of Electromagnetic Fields in Polyhedral Domains[END_REF][START_REF] Bonito | Regularity of the Maxwell Equations in Heterogeneous Media and Lipschitz Domains[END_REF] and [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF], proposition 1) as stated below:

Proposition 2.5. Let hypothesis 2.4 hold. There exists r max ∈]0, 1], called the regularity exponent, such that for all source terms S f ∈ L, the solution of problem 2.3 φ ∈ V belongs to 0≤r<rmax PH 1+r (R) (r max < 1) or PH 2 (R) (r max = 1) with continuous dependence:

∀r ∈ [0, r max [, φ PH 1+r (R) S f 0,R (r max < 1) or φ PH 2 (R) S f 0,R (r max = 1).
In all applications, the source term S f will be always in L, but in the theoretical study of the primal setting the source S f is taken in the larger function space V .

Source Problem

One way to deal with problem 2.3 is to solve its variational formulation. To establish the variational form, one can take ψ in V and use the dual product of (2.3) against ψ.

-div

(Dgrad x φ) + Σ r,0 φ, ψ V ,V = S f , ψ V ,V .
One can use Green's formula and remark that every term is in L.

R

Dgrad x φ.grad x ψ + R Σ r,0 φψ = S f , ψ V ,V
From now on, we define the following bilinear forms:

c p :    V × V → R (φ, ψ) → R Dgrad x φ.grad x ψ + R Σ r,0 φψ f p : V × V → R (S f , ψ) → S f , ψ V ,V
Therefore, the variational source problem reads:

Problem 2.6. For a given S f in V , find φ in V such that for all ψ in V c p (φ, ψ) = f p (S f , ψ) (2.4)
Remark 2.7. The solution here is looked for in H 1 0 (R) and div (Dgrad

x φ) is only in H -1 (R) because the source is in H -1 (R)
The uniqueness of the solution of problem 2.6 is ensured by the following theorem.

Theorem 2.8. For any source term S f in V , there exists a unique solution φ in V satisfying (2.4). Moreover, it stands

φ 1,R S f -1,R .
Proof. To prove this theorem, we use Lax-Milgram theorem. We verify all the hypothesis. Let ψ and ψ be in V , we can bound with the triangle inequality and hypothesis 2.4:

|c p (ψ, ψ )| ≤ D * R grad x ψ.grad x ψ + (Σ r,0 ) * R ψψ ≤ max(D * , (Σ r,0 ) * ) ψ 1,R ψ 1,R .
Thus, c p is continuous. In the same way, we prove that

f p (S f , •) is continuous |f p (S f , ψ)| = | S f , ψ V ,V | ≤ S f -1,R ψ 1,R .
The last hypothesis is the coercivity of c p . Let ψ be in V , thanks to hypothesis 2.4, it stands:

c p (ψ, ψ) = R D|grad x ψ| 2 + R Σ r,0 |ψ| 2 ≥ min(D * , (Σ r,0 ) * ) ψ 2 1,R .
Thus, one can apply Lax-Milgram theorem. Moreover we have

φ 1,R ≤ S f -1,R min(D * , (Σ r,0 ) * ) .
We already know that every solutions of the continuous problem 2.3 is solution of the variational problem 2.6. In fact, these two problems are equivalent.

Theorem 2.9. The solution of the variational problem 2.6 satisfies the continuous problem 2.3.

Proof. Let take ψ in D(R). As D(R) ⊂ V it stands R Dgrad x φ.grad x ψ + R Σ r,0 φψ = S f , ψ V ,V
Integrating by part, one can obtain

div (Dgrad x φ), ψ D (R),D(R) + Σ r,0 φ, ψ D (R),D(R) = S f , ψ D (R),D(R) .
Thus, φ verifies in D (R):

-div (Dgrad x φ) + Σ r,0 φ = S f .

We recall that φ is in V thus the boundary conditions are satisfied.

Eigenvalue Problem

In the same way as for the source problem, one can derive the eigenvalue variational formulation from problem 2.1, which reads:

Problem 2.10. Find (λ, φ) in R × V \{0} such that for all ψ in V c p (φ, ψ) = λ -1 f p (νΣ f φ, ψ).
(2.5)

Theorem 2.11. The eigenvalue problem 2.10 is equivalent to problem 2.1.

Proof. From problem 2.10, we know that (λ, φ) satisfies for all ψ in V :

c p (φ, ψ) = f p (λ -1 νΣ f φ, ψ).
We denote S f = λ -1 νΣ f φ, from hypothesis 2.4, S f is in L thus in V . Therefore, φ satisfies the source problem associated to S f . Then, one can conclude that (λ, φ) satisfies problem 2.1.

The following theorem allows us to introduce an operator that admits the same eigenpairs as the initial problem 2.1.

Theorem 2.12. There exists a unique compact operator T from V to V such that for all φ and ψ in V c p (T φ, ψ) = f p (νΣ f φ, ψ).

(2.6)

Proof. To prove this theorem, we apply the work of Babuška and Osborn in [START_REF] Babuška | Eigenvalue Problem[END_REF]. We already know that c p is a bilinear continuous coercive form onto V × V . The bilinear form

f p (i L 2 →H -1 (νΣ f •), •) is a continuous bilinear form on L × V . Moreover, one can remark that V is included in L with a compact embedding.
Thus there is an equivalence between the eigenpairs of the variational formulation and the ones of the operator T . The couple (λ -1 , φ) is solution of problem 2.10 if and only if the couple (λ, φ) is an eigenpair of T . We denote σ(T ) the spectrum of the operator T . The properties of σ(T ) are detailed in [START_REF] Brezis | Analyse Fonctionelle -Théorie et Applications[END_REF] (Section VI.3), we summarize them here. The set σ(T

) is a compact set included in [-T L(V ) , T L(V ) ]
. Moreover, 0 is in σ(T ) but is not an eigenvalue of T and every λ in σ(T )\{0} is an eigenvalue of T . Finally, T has no eigenvalue or it has a countable set of eigenvalues or its eigenvalues form a sequence which converges towards 0. One can remark that c p is equivalent to the natural inner product on V . From the symmetry of c p , the operator T is selfadjoint. Then, T is diagonalisable. We recall that from Krein-Rutman theorem 1.3, we know that the only physical solution of problem 2.1 is the fundamental mode.

Mixed Approach for the Diffusion Equation

The natural form of the neutron diffusion model is its mixed form. In this section, we study this setting. We need to introduce the Sobolev space H(div , R) which is defined as: H(div , R) = q ∈ L d , div q ∈ L . From now on Q stands for H(div , R). The function space Q is an Hilbert space with the inner product (•, •) div ,R which is defined for all q and p in Q:

(p, q) div ,R = R p • q + R div pdiv q.
The norm on Q induced by the inner product is denoted by • div ,R . Finally, we denote X the function space Q × L, this space is normed with the product norm • X , which for (q, ψ) in X is defined by: (q, ψ) 2 X = q 2 div ,R + ψ 2 0,R . From now on, we use the notations: ζ := (p, φ) and ξ := (q, ψ). As we did for the primal formulation, we study first the diffusion problem with null flux boundary conditions. We recall from §1.5 that this problem reads:

Problem 2.13. Find (λ, p, φ) ∈ R × Q × V \{0} such that D -1 p + grad x φ = 0 in R; div (p) + Σ r,0 φ = λ -1 νΣ f φ in R.
(2.7)

The source problem associated to the problem 2.13, for any source term S f ∈ L, is given by: Problem 2.14.

Find (p, φ) ∈ Q × V such that D -1 p + grad x φ = 0 in R; div (p) + Σ r,0 φ = S f in R. (2.8)
We suppose in this section, that the coefficients satisfy hypothesis 2.4.

Source Problem

We already know that problem 2.14 is equivalent to problem 2.3. Thus, from proposition 2.5, one can conclude that p has an extra regularity.

Proposition 2.15.

Let (D, Σ r,0 , νΣ f ) satisfy hypothesis 2.4. The current p in Q belongs also in 0≤r<rmax PH r (R) (r max < 1) or PH 1 (R) (r max = 1) with continuous dependence on the source term: ∀r ∈ [0, r max [, p PH r S f 0,R (r max < 1) or φ PH 1 (R) S f 0,R (r max = 1).
In order to study the well-posedness of problem 2.14, we consider its variational formulation. Let (q, ψ) be in X, by summing the first equation of (2.8) times -q and the second equation of (2.8) times ψ, and integrating the results over R, one obtains:

- R D -1 p • q + grad x φ • q + R (ψ div p + Σ r,0 φψ) = R S f ψ.
(2.9)

The second term, with the gradient of the flux, can be transformed with Green's formula. Indeed, for any q in Q, it stands:

R grad x φ • q + R φ div q = q • n, φ H -1/2 (∂R),H 1/2 (∂R)
.

(2.10) Moreover, as the flux φ is looked in V , φ is null on ∂R. Thus, it stands:

- R D -1 p • q + R φ div q + R ψ div p + R Σ r,0 φψ = R S f ψ. (2.11)
We define the bilinear forms:

a :    Q × Q → R (p, q) → R -D -1 p • q . (2.12) b :    Q × L → R (q, ψ) → R ψ div q . (2.13) t :    L × L → R (φ, ψ) → R Σ r,0 φ ψ .
(2.14) and:

c : X × X → R (ζ, ξ) → a(p, q) + b(q, φ) + b(p, ψ) + t(φ, ψ) .
(2.15)

We consider the linear form:

f :    X → R ξ → R S f ψ . (2.16)
Using equation (2.10) in equation (2.9) removes grad x φ, hence the regularity requirement on the flux φ can be lowered to φ in L. Therefore the variational problem reads:

Problem 2.16. For a given S f in L, find ζ ∈ X such that ∀ξ ∈ X:

c(ζ, ξ) = f (ξ).
(2.17)

Remark 2.17. In the primal setting of the problem, the source is taken in H -1 (R) and the flux φ is looked for in H 1 0 (R) (see remark 2.7). Here we need a source more regular, in L 2 (R) in order to have some regularity on the current p = -Dgrad x φ, p is in H(div , R) which is not the case for the primal setting.

Theorem 2.18. For a given S f in L, there exists a unique solution ζ in X to problem 2.16. moreover, it stands:

||ζ|| X ||S f || L .
Proof. To prove the well-posedness of the variational problem 2.16, we prove that c verifies an inf-sup condition: There exists β > 0 such that

inf ζ∈X ζ =0 sup ξ∈X ξ =0 |c(ζ, ξ)| ζ X ξ X ≥ β.
Given ζ = (p, φ) in X, one can remark that ξ = (q, ψ) = (-p, 1 2 φ + 1 2Σ r,0 div p) in X satisfies:

ξ 2 X = p 2 Q + 1 2 φ + 1 2Σ r,0 div p 2 0,R ≤ p 2 Q + 1 4 φ 2 0,R + 1 4 Σ -1 r,0 div p 2 0,R ≤ 1 + 1 4 ((Σ r,0 ) * ) -2 ζ 2 X .
The last line is obtained thanks to hypothesis 2.4.

One can obtain the bound:

c(ζ, ξ) = R D -1 |p| 2 + 1 2Σ r,0 |div p| 2 + Σ r,0 2 |φ| 2 ≥ (D * ) -1 p 2 0,R + 1 2(Σ r,0 ) * div p 2 0,R + (Σ r,0 ) * 2 φ 2 0,R ≥ min 1 D * ; 1 2(Σ r,0 ) * ; (Σ r,0 ) * 2 ζ 2 X .
To obtain the second line, we use hypothesis 2.4. With the bound of ξ X by ζ X , one can conclude the proof.

We already know that every solutions of the continuous problem 2.14 is solution of the variational problem 2.16. In fact, these two problems are equivalent.

Theorem 2.19. The solution of the variational problem 2.16 satisfies the continuous problem 2.14.

Proof. Let ζ = (p, φ) be the solution of problem 2.16. For any q in D(R) d , from equation (2.17) with (q, 0) as test function, one can obtain:

- R D -1 p • q + R φ div q = 0. (2.18)
Thanks to Green's formula, the second integral can be transformed into:

R φ div q = -grad x φ, q D ,D .
Thus one can conclude that grad x φ = -D -1 p in the sense of distribution. The set

D(R) d is dense in L d [GiRa86], then grad x φ = -D -1 p in L d . Thus φ is in H 1 (R).
In order to find the boundary condition on the flux φ, we apply (2.17) to (q, 0), with q in Q, and we use Green's formula

- R D -1 p • q + R φ div q = - R D -1 p • q - R grad x φ • q + q • n, φ H -1 2 ,H 1 2 = 0.
As grad x φ = -D -1 p, for all q ∈ Q, it stands:

q • n, φ H -1 2 (∂R),H 1 2 (∂R) = 0. The application q ∈ Q → q |∂R • n ∈ H -1 2 (∂R)
is surjective, thus for any τ ∈ H -1/2 (∂R) it stands:

τ, φ H -1 2 ,H 1 2 = 0.
Therefore, φ |∂R is null and the flux is in H 1 0 (R). By applying (2.17) to (0, ψ), with ψ in L 2 (R), one can obtain easily the second equation of system (2.8).

Eigenvalue Problem

Let us focus on the generalized eigenvalue problem 2.13. To approximate the solution (λ, p, φ) of this problem, one can study the direct operator, or the inverse operator, which associates the solution ζ ∈ X to a source S f ∈ L. As the fission macroscopic cross section can be null somewhere, we need to study the last one. As the flux solution of problem 2.13 is more regular than L, we define the inverse operator onto a more regular space. Let 0 < µ < 1/2 be given, we introduce the inverse operator B µ associated to the source problem 2.14:

given f ∈ H µ (R), we call B µ f = φ in H 1 (R) the flux that solves problem 2.14 with source S f = νΣ f f . Lemma 2.20. B µ is a compact operator from H µ (R) to H µ (R). Proof. Since νΣ f satisfy hypothesis 2.4, it holds S f µ,R f µ,R because µ < 1/2. Hence, B µ is a bounded operator from H µ (R) to itself. Indeed, first by continuous embedding of H 1 (R) in H µ (R) it stands: B µ f µ,R B µ f 1,R .
Moreover, the solution φ depends continuously on the source term, the it comes:

B µ f 1,R = ||φ|| 1,R S f 0,R .
To finish, H µ (R) is compactly embedded in L:

S f 0,R S f µ,R .
Then, one can conclude that:

B µ f µ,R f µ,R .
In addition, since the eigenfunction actually belongs to H 1 (R) (proposition 2.5) with continuous dependence ( φ 1,R S f µ,R ), it follows that B µ is a compact operator.

Looking for the eigenpair (λ -1 , (p, φ)) of problem 2.13 is the same as looking for the eigenpair of B µ (λ, φ), setting p = -D -1 grad x φ. Moreover as B µ is a compact operator, its eigenvalues are bounded, non-null and countable. We recall that from Krein-Rutman theorem 1.3, we know that the only physical solution of problem 2.13 is the fundamental mode.

On the Regularity of the Solution

In the primal formulation of the neutron diffusion source problem 2.3, the source term is taken in V . Therefore the flux solution of problem 2.3 is in V , and the current, p = -Dgrad x , is only in L d and its divergence is in (V ) d . If the source term is taken in L, the current is more regular, p is in Q.

Under the mixed setting of the neutron diffusion problem 2.14, the source is taken is L, and the solution (φ, p) is in V × Q. But one can lower the regularity requirement in the variational formulation on φ to φ ∈ L.

The choice of the solving the primal or the mixed setting depends on which component, φ or p, one wants to have the best precision of its discrete counterpart.

Chapter 3

Extension to the Multigroup SP N Transport Equations

Introduction

This model takes into account the energy dependence and the angular dependence. This is done by using the multigroup approximation and by projecting the problem on the spherical harmonics § 1.5. We recall that G is the number of energy groups and N the number of even and odd moments. The unknowns are the moments of the current and the flux. To describe them, we introduce the following function spaces:

Q = (Q N ) G ; V = (V N ) G ; L = (L N ) G ; X = Q × L and L = (L N ) G .
For W one of these previous function spaces, we denote the natural product norm • W . We recall from (1.33), the multigroup SP N transport problem, under its mixed formulation, reads:

Problem 3.1. Find (λ, p, φ) ∈ R × Q × V \{0}, such that:    T o p + grad x (Hφ) = 0 in R; T Hdiv (p) + T e φ = λ -1 M f φ in R; (3.1)
The source problem associated to problem 3.1 reads:

Problem 3.2. For a given source S f ∈ L, find (p, φ) ∈ Q × V such that:    T o p + grad x Hφ = 0 in R; T Hdiv p + T e φ = S f in R. (3.2)
We introduce three lemmas which are used in this chapter.

Lemma 3.3. For any function ψ in V and q in Q, it stands:

grad x (Hψ) = Hgrad x ψ;
T Hdiv (q) = div ( T Hq).

The matrix

H = diag (H d ) G g=1 where H d is in R d N ×d N
and is defined by:

H d =         I d I d 0 • • • 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 . . . . . . . . . I d 0 • • • • • • 0 I d         .
Proof. This can be proved easily by evaluating the two sides of the equality.

Lemma 3.4. There exists α H , β H > 0 such that for any ψ in L:

α H ψ L ≤ Hψ L ≤ β H ψ L . (3.3)
Lemma 3.5. There exists α H , β H > 0 such that for any q in L:

α H ψ L ≤ Hψ L ≤ β H q L . (3.4)
The proof of lemma 3.4 is given in appendix C. We denote To = H -1 T o T H -1 . As the matrices H and T o are invertible, the matrix To is invertible too and T-1

o = T HT -1
o H. Thanks to lemma 3.4, the primal form of the multigroup SP N transport equations (1.44) reads:

Problem 3.6. Find (λ, φ) ∈ R × V \{0}, such that: -div ( T-1 o grad x φ) + T e φ = λ -1 M f φ. (3.5)
The source problem associated to problem 3.6 reads:

Problem 3.7. For a given S f ∈ V , find φ ∈ V , such that:

-div ( T-1 o grad x φ) + T e φ = S f . (3.6)
We recall that the flux φ is equal to

(φ 1 , • • • , φ G ) T with φ g = (φ g 0 , φ g 2 , • • • , φ g N -1 ) T , for 1 ≤ g ≤ G.
Hence, we cannot generalize the regularity property given for the one-group diffusion problem 2.3. As a matter of fact, in proposition 2.5 the symmetry of the tensor is required (hypothesis 2.4), whereas the diagonal energy blocks of the matrix T-1 o are not symmetric. Instead, in order to obtain error estimates, we suppose that the solution of problem 3.7 has an extra regularity.

Hypothesis 3.8. We suppose that there exists r max ∈]0, 1] such that for all source terms S f ∈ L, the flux solution of problem 3.7 φ ∈ V belongs to 0≤r<rmax P((

H 1+r (R)) N ) G (r max < 1) or P((H 2 (R)) N ) G (r max = 1) with continuous dependence: ∀r ∈ [0, r max [, φ P((H 1+r (R)) N ) G S f L (r max < 1) or φ P((H 2 (R)) N ) G S f L (r max = 1).

Primal Approach

In this section, we study the primal formulation of the multigroup SP N equations. For this study, we suppose that the macroscopic cross sections satisfy condition 1.15 on page 27, which allows to ensure the continuity and the coercivity of the bilinear form associated to problem 3.7 (see proposition 1.16.

3.2 : Primal Approach

Source Problem

As we did for the diffusion one-group model in primal form, we study problem 3.7 through its variational formulation. Let ψ be in V , we multiply equation (3.6) by ψ and we integrate it over R:

R -div ( T-1 o grad x φ) • ψ + R T e φ • ψ = S f , ψ V ,V .
One can use Green's formula and obtain:

R T-1 o grad x φ grad x ψ + R T e φ • ψ = S f , ψ V ,V .
From now on, we define the following bilinear forms:

c s,p :    V × V → R (φ, ψ) → R T-1 o grad x φ grad x ψ + R T e φ • ψ f s,p : V × V → R (S f , ψ) → S f , ψ V ,V
Therefore, the variational source problem reads:

Problem 3.9. For a given S f in V , find φ in V such that for all ψ in V it stands:

c s,p (φ, ψ) = f s,p (S f , ψ). (3.7)
The uniqueness of the solution of problem 3.9 is ensured by the following theorem.

Theorem 3.10. Under condition 1.15, for any source term S f in V , there exists a unique solution φ in V satisfying problem 3.9. Moreover, it stands

φ V S f V .
Proof. To prove this theorem, we use Lax-Milgram theorem. We verify all the hypothesis. One can easily prove that f s,p (S f , •) is continuous. Let φ be in V . From the definition of c s,p and the definition of To , it stands:

c s,p (φ, φ) = R T HT -1 o Hgrad x φ grad x φ + R T e φ • φ. (3.8)
The first term in the left hand side of (3.8) is bounded by:

R T HT -1 o Hgrad x φ grad x φ = R T -1 o Hgrad x φ Hgrad x φ; ≥ (T -1 o ) * Hgrad x φ 2 L ; ≥ (T -1 o ) * α 2 H grad x φ 2 L .
In the second line we use proposition 1.16 and in the third line we use lemma 3.5. Using proposition 1.16, one can bound the second term in the left hand side of (3.8) as:

R T e φ • φ ≥ (T e ) * φ 2 L .
Thus, the bilinear form c s,p is coercive. We recall from proposition 1.16 that T -1 o is bounded and from lemma 3.5 that H is bounded too. Therefore, T-1 o is bounded. Furthermore, we know from proposition 1.16, that T e is bounded. Then, c s,p is continuous. Finally we can apply Lax-Milgram theorem.

Moreover, we have the following equivalence.

Theorem 3.11. The solution of problem 3.9 satisfies problem 3.7.

Proof. The proof use the same methodology of the proof of theorem 2.9.

Eigenvalue Problem

Now we consider the eigenvalue problem 3.6. We use the same methodology as for the diffusion model (cf. 2). The variational formulation of problem 3.6 reads: Problem 3.12. Find (λ, φ) ∈ R × V \{0} such that for all ψ in V it stands:

c s,p (φ, ψ) = λ -1 f s,p (M f φ, ψ).
(3.9)

The following theorem allows us to introduce an operator that admits the same eigenpairs as the initial problem 2.1.

Theorem 3.13. There exists a unique compact operator T from V to V such that for all φ and ψ in V it stands: c s,p (T φ, ψ) = f s,p (M f φ, ψ).

(3.10)

Proof. We use the same argument as in the proof of theorem 2.12 (for the primal form of the neutron diffusion equations).

Thus there is an equivalence between the eigenpairs of the variational formulation and the ones of the operator T . The couple (λ -1 , φ) is solution of problem 3.12 if and only if the couple (λ, φ) is an eigenpair of T . We recall from §2.1.2 that for the neutron diffusion model under its primal setting the operator T is diagonalisable because the bilinear form c p is symmetric. In the case of the multigroup SP N model under its primal setting, c s,p is not symmetric thus T is not diagonalisable. Thus there is a countable number of non-null eigenvalues for problem 3.6, moreover the eigenvalues are bounded. We recall that from Krein-Rutman theorem 1.3, we know that the only physical solution of problem 3.6 is the fundamental mode.

Mixed Approach

In this section, we study the mixed formulation of the multigroup SP N equations. For this study, we suppose that the macroscopic cross sections satisfy condition 1.13, which allows to ensure the inf-sup condition on the bilinear form associated to problem 3.2.

Source Problem

To derive the variational formulation of problem 3.2, we take (q, ψ) in X then we multiply by -q and integrate over R the first equation in (3.2), we multiply by ψ and integrate over R the second equation in (3.2) and we sum the volume integrals.

-

R T o p q + grad x (Hφ) q + R ψ • T Hdiv p + T e φ • ψ = R S f • ψ. (3.11)
Thanks to Green's formula, we can change the integral with grad x φ.

Proposition 3.14. For any ψ in V and q in Q, it stands: R grad (Hψ) q + R ψ • T Hdiv q = 0.

(3.12) Therefore, equation (3.11) can be transformed in:

- R T o p q + R φ • T Hdiv q + R ψ • T Hdiv p + R T e φ • ψ = R S f • ψ. (3.13)
Thus one can look for the flux φ in L instead of V . We define the bilinear forms:

a s :    Q × Q → R (p, q) → R -T o p q . (3.14) b s :    Q × L → R (q, ψ) → R ψ • T Hdiv q .
(3.15)

t s :    L × L → R (φ, ψ) → R T e φ • ψ .
(3.16) and:

c s : X × X → R (ζ, ξ) → a s (p, q) + b s (q, φ) + b s (p, ψ) + t s (φ, ψ) .
(3.17)

We consider the linear form:

f s :    X → R ξ → R S f • ψ . (3.18)
Therefore the variational problem reads:

Problem 3.15. Find ζ ∈ X such that for all ξ ∈ X c s (ζ, ξ) = f s (ξ). (3.19)
Theorem 3.16. For any source term S f in L, there exists a unique solution ζ in X satisfying problem 3.15. Moreover, it stands

ζ X S f L .
Proof. To prove the claim, one looks for an inf-sup condition and a solvability condition [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Ern | Theory and Pratice of Finite Elements[END_REF] to ensure well-posedness of problem 3.15. The solvability condition writes The set ξ ∈ X | ∀ζ ∈ X, c(ζ, ξ) = 0 is equal to {0}.

. First, we prove the following inf-sup condition: There exists β > 0 such that

inf ζ∈X ζ =0 sup ξ∈X ξ =0 |c(ζ, ξ)| ζ X ξ X ≥ β. (3.20) Given ζ in X, one can remark that ξ = (-p, 1 2 φ + 1 2 T (T -1 e ) T Hdiv p) in X satisfies: ξ 2 X = p 2 Q + 1 2 φ + 1 2 T (T -1 e ) T Hdiv p 2 L .
Using the triangular inequality, one obtains:

ξ 2 X ≤ p 2 Q + 1 4 φ 2 L + 1 4 T (T -1 e ) T Hdiv p 2 L .
According to proposition 1.14, T -1 e satisfies the following boundedness condition, for all ψ ∈ L:

T -1 e ψ L ≤ (T -1 e ) * ψ L . Furthermore, the matrices T -1
e and T (T -1 e ) satisfy the same boundedness condition, thus for all ψ ∈ L, it stands

T (T -1 e )ψ L ≤ (T -1 e ) * ψ L .
Moreover, thanks to lemma 3.4, one obtains:

ξ 2 X ≤ 1 + 1 4 ((T -1 e ) * ) 2 β 2 H ζ 2 X . Now we look for a lower bound to c(ζ, ξ). c(ζ, ξ) = R T o p p + 1 2 T (T -1 e ) T Hdiv p • T Hdiv p + 1 2 T e φ • φ .
From proposition 1.14, one obtains the following bound:

c(ζ, ξ) ≥ (T o ) * p 2 L + 1 2 (T -1 e ) * T Hdiv p 2 L + 1 2 (T e ) * φ 2 L . Using lemma 3.3, one obtains: c(ζ, ξ) ≥ min (T o ) * ; 1 2 (T -1 e ) * α 2 H ; 1 2 (T e ) * ζ 2 X .
With the bound of ξ X by ζ X , one obtains the inf-sup condition in 3.20. Let ξ be in ξ ∈ X | ∀ζ ∈ X, c(ζ, ξ) = 0 and ξ = 0, we denote ζ = (p, φ) with p = -q and φ = 1 2 ψ + 1 2 (T e ) -1T Hdiv q. In the same way as for the inf-sup condition, we find the following bound:

c(ζ, ξ) ≥ min (T o ) * ; 1 2 (T -1 e ) * α 2 H ; 1 2 (T e ) * ξ 2 X .
As ξ is non null, we have c(ζ, ξ) > 0. By definition of ξ, c(ζ, ξ) = 0. This contradiction implies that ξ = 0.

We already know that every solution of the continuous problem 3.2 is solution of the variational problem 3.15. In fact, these two problems are equivalent.

Theorem 3.17. The solution of the variational problem 3.15 satisfies the continuous problem 3.2 (for the mixed form of the neutron diffusion equations).

Proof. The proof is the same as the one for theorem 2.19.

Eigenvalue Problem

Let us focus on the approximation of the generalized eigenvalue problem 3.1. In our lowregularity setting, we supplement the assumptions of condition 1.13 with the following condition:

Condition 3.18. For all 1 ≤ g ≤ G, it holds: νΣ g f ∈ PW 1,∞ (R). We use the same methodology as for the one-group diffusion model (cf. 2.2): we study the inverse operator since the right hand side of problem 3.1 may vanish locally. For 0 ≤ µ < 1/2, we denote H µ = ((H µ (R)) N ) G . As the flux solution of problem 3.1 belongs to a more regular space than L, we define the inverse operator onto a more regular space. Let 0 ≤ µ < 1/2 be given, we introduce the inverse operator B µ associated to the source problem 3.2: given f ∈ H µ , we call

B µ f = φ in V the flux that solves (3.2) with source S f = M f f . Therefore, under condition 3.18, for all f ∈ H µ , M f f is in H µ . Lemma 3.19. B µ is a compact operator from H µ to H µ . Proof. We recall that PH µ (R) = H µ (R) because µ < 1/2. Since M f satisfies condi- tion 3.18, it holds S f H µ f H µ .
Hence, B µ is a bounded operator from H µ to itself. Indeed, first by continuous embedding of V in H µ it stands:

B µ f H µ B µ f V .
Moreover, the solution φ depends continuously on the source term, the it comes:

B µ f V S f L .
To finish, H µ is continuous embedded in L:

S f L S f H µ .
Then, one can conclude that:

B µ f H µ f H µ .
In addition, since the eigenfunction actually belongs to V with continuous dependence ( φ V S f H µ ), it follows that B µ is a compact operator.

Then, looking for the eigenvalue λ of problem 3.1 is the same as looking for the inverse of the eigenvalue of B µ which has non-null, bounded and countable eigenvalues. We recall that from Krein-Rutman theorem 1.3, we know that the only physical solution of problem 3.1 is the fundamental mode.

Part II

Discrete Problem

The study of the convergence of the multigroup SP N approximation let us know the behaviour of the implementation and it helps us to calibrate the solver. Moreover, the behaviour of the method is important to know in order to study the uncertainty propagation [START_REF] Diop | Cross-Section Uncertainty Propagation in the Multigroup Slowing-Down Equations Using Probability Table Formalism[END_REF]. Indeed, in [START_REF] Charrier | Strong and weak error estimates for the solutions of elliptic partial differential equations with random coefficients[END_REF], the author proposes an a priori error estimate for diffusion problem with random coefficients. This estimate is composed of two parts, one is the stochastic error which measures the uncertainty on the solution due to the uncertainty on the coefficient, the second part is a deterministic error due to the discretization of the diffusion. The numerical analysis proposed here can be used in order to determine the deterministic error of the method. The approximation of the diffusion problem with a source term, written under its primal setting (see problem 2.3) has been widely studied. The reader can refer to [START_REF] Crank | The Mathematics of Diffusion[END_REF][START_REF] Ern | Theory and Pratice of Finite Elements[END_REF] for instance. Concerning the eigenvalue problem, the primal approximation was studied by Osborn et al in [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF][START_REF] Babuška | Eigenvalue Problem[END_REF].

Concerning the mixed setting (see problem 3.2), one can refer to the work of Brezzi and Fortin in [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] for the source problem. Recently, the authors of [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF] consider the eigenvalue diffusion problem without perturbation term, which corresponds to vanishing even removal cross sections i.e. T e = 0 in problem 3.1 (obviously this never occurs in our case), and in [START_REF] Boffi | Optimal Convergence of Adaptive FEM for Eigenvalue Clusters in Mixed Form[END_REF] the authors proposed a method to derive an optimal rate of convergence for mixed eigenvalue problem without pertubation term. The numerical analysis of the mixed approximation of the diffusion with low-regularity solution has been carried out recently in [START_REF] Ciarlet | Numerical Analysis of the Mixed Finite Element Method for the Neutron Diffusion Eigenproblem with Heterogeneous Coefficients[END_REF].

We will treat the use of a domain decomposition method with the L 2 -jump domain decomposition method, which is described in [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF]. In case of using RTN finite elements, its numerical analysis is carried out in [START_REF] Ciarlet | Numerical Analysis of the Mixed Finite Element Method for the Neutron Diffusion Eigenproblem with Heterogeneous Coefficients[END_REF]. This part is organized as follow:

• in chapter 4, we extend the classical results of the approximation of the diffusion equations to the multigroup SP N transport equations under their primal settings;

• in chapter 5, we present the study published in the first part of [START_REF] Ciarlet | Numerical Analysis of the Mixed Finite Element Method for the Neutron Diffusion Eigenproblem with Heterogeneous Coefficients[END_REF] and we extend it to the multigroup SP N transport equations under their mixed setting;

• in chapter 6, we present the last part of [START_REF] Ciarlet | Numerical Analysis of the Mixed Finite Element Method for the Neutron Diffusion Eigenproblem with Heterogeneous Coefficients[END_REF] which concerns the numerical analysis of the L 2 -jump domain decomposition method for the diffusion problem.

Chapter 4

Primal Resolution

In this chapter we discretize the continuous source problem 3.7 and the continuous eigenvalue problem 3.6, set in the function space V = ((H 1 0 (R)) N ) G . We recall the condition 1.15 on page 27: Condition 4.1. For all energy groups 1 ≤ g, g ≤ G, g = g and for all 0 ≤ n ≤ N ,it stands:

           (Σ g r,n , Σ g →g s,n , νΣ g f ) ∈ PW 1,∞ (R) × L ∞ (R) × L ∞ (R), ∃(Σ r,n ) * , (Σ r,n ) * > 0, 0 < (Σ r,n ) * ≤ t n Σ g r,n ≤ (Σ r,n ) * a.e. in R, ∃0 < ε < 1 G -1 , |Σ g →g s,n | ≤ εΣ g r,n
a.e. in R 0 ≤ νΣ g f a.e. in R, ∃g, g s.t. χ g,g νΣ g f = 0.

Moreover, we suppose that:

   α s,o (G -1) < 1 1 + α r,o , α s,e α r,e (G -1) < 1.
We suppose that condition 4.1 is satisfied. We recall from proposition 1.16 that under condition 4.1 the matrices T -1 o and T e are bounded and have positive properties which ensure the well-posedness of problem 3.12. The discretization proposed here is a H 1 -conforming finite element method, i.e. the space of discretization is included in the space where the solution is looked for. As the diffusion model and the multigroup diffusion model are specific cases of the SP N model, we only prove the well-posedness of the discrete problem and find a priori error estimates in the more general model. First in § 4.1, we introduce some notation and definition, then in § 4.2 we study the approximation of the source problem. Finally in § 4.3 we consider the eigenvalue problem.

Discrete Spaces

We consider here a H 1 -conforming finite element method. We denote by (T h ) a family of triangulations of the reactor R. The subscript h of a triangulation T h represents the greatest diameter of its elements:

h = diam K∈T h (K). (4.1)
We assume that the triangulation family (T h ) h is shape-regular (cf definition 1.107, [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF]).

Definition 4.2. A triangulation family (T h ) h is said to be shape-regular if there exists σ 0 > 0 such that:

∀h > 0, ∀K ∈ T h , h K ≤ σ 0 ρ K ,
where ρ K is the diameter of the largest ball that can be inscribed in K.

We also suppose that the macroscopic cross sections satisfy the following condition:

Condition 4.3. A function Σ ∈ PW 1,∞ (R)
is said to be smooth on each element of a triangulation T h if it satisfies:

∀K ∈ T h , Σ ∈ W 1,∞ (K).
For k ∈ N * , the set of all the polynomials with degree less than or equal to k is denoted P k . Moreover, we define the space V k h ⊂ V by:

V k h = ψ ∈ V, ∀K ∈ T h , ψ |K ∈ P k . (4.2) We define V k h = ((V k h ) N ) G
, which is the approximation space for the space V . Moreover, we denote π k p the orthogonal projection from V to V k h , π k p is defined by:

∀ψ ∈ V , ∀ψ h ∈ V k h , (ψ -π k p ψ, ψ h ) V = 0,
From ([ErGu04], proposition 1.134) it stands, for all r > 0:

∀ψ ∈ ((H 1+r (R)) N ) G , ψ -π k p ψ V h min(k,r) ψ V (4.3)
This estimation give us the following approximability property:

Proposition 4.4. The finite dimension space V k h has the following approximability property:

∀ψ ∈ V , lim h→0 inf ψ h ∈V k h ψ -ψ h V = 0.
The solution of the variational problem 3.9 and problem 3.12 are approximated in V k h .

Source Problem

From the variational formulation problem 3.9, we derive the discrete variational formulation problem, which reads:

Problem 4.5. For a given source

S f in V , find φ h in V k h such that for all ψ h in V k h , such that: c s,p (φ h , ψ h ) = f s,p (S f , ψ h ). (4.4)
As V k h is included in V , all the properties of the bilinear forms c s,p and f s,p are still true on V k h . Therefore, problem 4.5 is well-posed. The proof is the same as theorem 3.10. Now we look for an error estimate on the discrete solution.

Theorem 4.6. Suppose that there exists r max in [0, 1] such that ∀r ∈ [0, r max [, φ ∈ ((H 1+r (R)) N ) G . Then the solution of problem 4.5, φ h converges towards the solution φ in V and it holds:

φ -φ h V h ω S f V , (4.5)
where ω = min(r max , k).

Proof. From Céa's lemma as detailed in ([ErGu04], §2.3).

We also find an a priori error estimate in the norm of L:

Theorem 4.7. Suppose that there exists

r max in [0, 1] such that ∀r ∈ [0, r max [, φ ∈ ((H 1+r (R) N ) G .
Then the solution of problem 4.5, φ h converges in L towards the solution φ and it holds:

φ -φ h L h 2ω S f V , (4.6)
where ω = min(r max , k).

Proof. From Aubin-Nitsche lemma as detailed in ([ErGu04], §2.3).

Eigenvalue Problem

The discrete counterpart of problem 3.12 reads:

Problem 4.8. Find (λ h , φ h ) in R × V k h \{0} such that for all ψ h in V k h , it stands: c s,p (φ h , ψ h ) = λ -1 h f s,p (M f φ h , ψ h ). (4.7)
As we consider here a H 1 -conforming finite element method, the discrete counterpart of theorem 3.13 reads:

Theorem 4.9. There exists a unique compact operator T h from V k h to V k h such that for all φ h and ψ h in V k h it stands:

c s,p (T h φ h , ψ h ) = f s,p (M f φ h , ψ h ). (4.8)
Proof. We use the same methodology as the proof of theorem 2.12 in the continuous case.

We denote by σ(T h ) the spectrum of the operator T h . As for the continuous problem, λ h is in σ(T h ) if and only if λ -1 h is an eigenvalue of (4.7). Moreover, the operator T can be written P h T where P h is the projection from V onto V k h such that it stands:

∀φ ∈ V , ∀ψ h ∈ V k h , c s,p (P h φ, ψ h ) = c s,p (φ, ψ h ).
The sequence of the operator (P h ) h is pointwise converging in L(V ) towards the identity operator, so the sequence (T h ) h is pointwise converging towards T . Moreover T h is a compact operator thus, the sequence (T h ) h is converging in L(V ) towards T . The norm convergence guarantees that there is no spectral pollution (see [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF]):

• Given any closed, non-empty disk D ⊂ C such that D ∩ σ(T ) = ∅, there exists h 0 such that, for all h < h 0 , D ∩ σ(T h ) = ∅.

• Given any closed, non-empty disk D ⊂ C such that D ∩ σ(T ) = {λ}, with λ of multiplicity m λ , there exists h 0 such that, for all h < h 0 , D ∩ σ(T h ) contains exactly m λ discrete eigenvalues. Now, thanks to the work of Babuška and Osborn in [START_REF] Babuška | Eigenvalue Problem[END_REF], we find an a priori error estimate on the eigenvalues.

Theorem 4.10. We denote by (λ, φ) (resp. (λ h , φ h )) the solution of problem 3.12 (resp 4.8).

Moreover, we denote by ω λ the regularity of the eigenfunction φ (φ ∈ ((

H 1+ω λ (R)) N ) G ),
and ω = min(ω λ , k). The following a priori error estimate holds:

|λ -λ h | h 2ω . (4.9)
Proof. We apply theorem 8.3 in [START_REF] Babuška | Eigenvalue Problem[END_REF].

Chapter 5

Mixed Resolution

In this chapter, we study in section 5.1 the discretization of the diffusion equations (1.73). This section comes from [START_REF] Ciarlet | Numerical Analysis of the Mixed Finite Element Method for the Neutron Diffusion Eigenproblem with Heterogeneous Coefficients[END_REF]. In section 5.2, we extend the results on the diffusion equations to the multigroup SP N neutron transport equations.

The Neutron Diffusion Equations

As we did for the continuous study, we first consider the neutron diffusion model. We recall that the functional spaces used for the mixed setting are defined as follows:

Q = H(div , R) , q Q := q div ,R ; X = { ξ := (q, ψ) ∈ Q × L} , ξ X := q 2 div ,R + ψ 2 0,R 1/2 .
The conditions on the coefficients defining the model are those in hypothesis 2.4.

Discretization

We study conforming discretizations of the variational formulation (2.17). To fix ideas, we use a family of triangulations (T h ) h , indexed by a parameter h, which is classically chosen as the largest diameter of elements of the triangulation. We introduce discrete, finite-dimensional, spaces indexed by h as follows:

Q h ⊂ H(div , R), and L h ⊂ L.
For approximation purposes, and following Definition 2.14 in [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF], we assume that (Q h ) h , resp. (L h ) h have the approximability property in the sense that

∀q ∈ H(div , R), lim h→0 inf q h ∈Q h q -q h div ,R = 0, ∀ψ ∈ L, lim h→0 inf ψ h ∈L h ψ -ψ h 0,R = 0, (5.1)
and also that L h includes the subspace L 0 h of piecewise constant fields on the triangulation. We impose: div

Q h ⊂ L h .
We endow Q h with the norm • div ,R , while L h is endowed with • 0,R . We finally define:

X h = { ξ h := (q h , ψ h ) ∈ Q h × L h } , endowed with • X .
The conforming discretization of the variational formulation (2.17) reads:

Find (p h , φ h ) ∈ X h , such that ∀(q h , ψ h ) ∈ X h : a(p h , q h ) + b(q h , φ h ) + b(p h , ψ h ) + t(φ h , ψ h ) = (S f , ψ h ) 0,R .
(5.2)

Or equivalently:

Find ζ h ∈ X h such that ∀ξ h ∈ X h , c(ζ h , ξ h ) = f (ξ h ).
(5.3)

For later use, we denote π 0 the L orthogonal projector on its subspace L 0 h . By construction, it holds range(π 0 ) = L 0 h where π 0 is defined by:

∀ψ ∈ L, ∀ψ h ∈ L 0 h , (π 0 ψ -ψ, ψ h ) 0,R = 0. According to [ErGu04, Proposition 1.135]: ∀z ∈ L, z -π 0 z 0,R z L , ∀z ∈ PH 1 (R), z -π 0 z 0,R h z PH 1 (R) , ∀z ∈ PW 1,∞ (R), z -π 0 z ∞,R h z PW 1,∞ (R) .
(5.4)

Similar results hold on subsets of R, provided the discretizations are conforming.

Discrete inf-sup Condition

The discrete inf-sup condition to be found writes:

∃ η h > 0, inf ζ h ∈X h sup ξ h ∈X h c(ζ h , ξ h ) ζ h X ξ h X ≥ η h .
(5.5) Once (5.5) is achieved, one obtains existence and uniqueness of the discrete solution ζ h , hence the corresponding linear system is well-posed. Generally, η h depends on h, but our aim is to obtain that (η h ) h is uniformly bounded away from 0. In this sense, one has at hand a uniform discrete inf-sup condition (udisc), from which the error analysis can classically be derived. In the case where (η h ) h is uniformly bounded, we keep the index h in order to keep in mind that we consider the discrete udisc.

Theorem 5.1. Let D, resp. Σ a ∈ PW 1,∞ (R), satisfy hypothesis 2.4. The discrete infsup condition (5.5) is fulfilled. Moreover, it is a uniform discrete inf-sup condition with respect to h and k.

Proof. In order to prove the discrete inf-sup condition, we use the same method as for the continuous inf-sup condition (cf. proof of Theorem 2.18). One can remark that if Σ a is piecewise-constant, 1 2 Σ -1 a div p h is automatically in L h . Otherwise, we project Σ -1 a on the piecewise-constant functions. One can choose:

q h = -p h ∈ Q h , ψ h = 1 2 φ h + 1 2 π 0 ((Σ a ) -1 )div p h ∈ L h .
Using (5.4) with z = (Σ a ) -1 yields (Σ a ) -1π 0 ((Σ a ) -1 ) ∞,R h, which allows us to derive again a udisc in this more general case.

Numerical Analysis of the Source Problem

We consider the neutron diffusion equation assuming that D, resp. Σ a ∈ PW 1,∞ (R), satisfy hypothesis 2.4. Under the assumptions of §5.1.1, it follows from the previous study that lim h→0 ζζ h X = 0. We find below a sharper bound of the error ζζ h X by using Proposition 2.5. In order to obtain optimal a priori error estimates, we must know the regularity of the solution to problem initial equation. Since we have assumed that the source term S f belongs to L, we already know that φ 1,R S f 0,R . Moreover, under the assumptions of Proposition 2.5, the solution φ has some extra regularity, and the low-regularity case corresponds to r max < 1/2 there. This is the case that we are focusing on now. In this setting, the field p := -D grad φ automatically belongs to PH r (R), for 0 ≤ r < r max . We suppose in addition that

∃µ ∈]0, r max [, S f ∈ PH µ (R).
Then we have div p ∈ PH µ (R) (recall PH µ (R) = H µ (R) for µ < 1/2). We will use this hypothesis on S f to carry on the calculations of the error estimates. We recall below the definition of the Raviart-Thomas-Nédélec (or RTN) finite element [START_REF] Raviart | Primal Hybrid Finite Element Methods for 2nd Order Elliptic Equations[END_REF][START_REF] Nédélec | Mixed Finite Elements in R 3[END_REF]. Let (K ) 1≤ ≤L be a conforming mesh, or triangulation, of R made of parallelepipeds (a mesh, or triangulation, is said to be conforming if in every K , D and Σ a are smooth). Let P (K ) be the set of polynomials defined over K . For integer values l, m, p ≥ 0, we consider the following subspace of P (K ):

Q l,m,p (K ) = q(x, y, z) ∈ P (K ) | q(x, y, z) = l,m,p e,j,k=0
a e,j,k x e y j z k , a e,j,k ∈ R .

For integer k ≥ 0, let us set k = k + 1 and introduce the vector polynomial space:

D k (K ) = [Q k ,k,k (K ) × 0 × 0] ⊕ [0 × Q k,k ,k (K ) × 0] ⊕ [0 × 0 × Q k,k,k (K )].
We can now define the RTN [k] finite element subspace of H(div , R) × L:

Q k h = q ∈ H(div , R) | ∀ ∈ {1, ..., L}, q |K ∈ D k (K ) , L k h = ψ ∈ L | ∀ ∈ {1, ..., L}, ψ |K ∈ Q k,k,k (K ) .
(5.6)

As required, it holds div

Q k h ⊂ L k h and L 0 h ⊂ L k h . We recall that for any q in H(div , R), its RTN [k] -interpolant q k R ∈ Q k h satisfies: ∀ψ h ∈ L k h , b(q -q k R , ψ h ) = 0.
(5.7)

In addition thanks to the commuting diagram property, cf. [BoBF13, §2.5.2], it holds ∀q ∈ H(div , R), div q 0 R = π 0 (div q).

(5.8)

Let q ∈ H r (R), such that div q ∈ H s (R), 0 < r, s < r max . According to [BGNR06, Lemma 3.3]: q -q 0 R 0,R (h r |q| r,R + h div q 0,R ) , div (q -q 0 R ) 0,R h s |div q| s,R .
(5.9)

Similar results hold on subsets of R, provided the discretizations are conforming.

Remark 5.2. If one chooses another discretization, all results presented hereafter hold provided the estimates (5.9) remain true. For instance, for the RTN [k] finite element defined on tetrahedral triangulations of R, cf. [BoBF13, §2.3.1]. To prove (5.9) in this case, one has simply to apply the results of [BGNR06, §3.2]. On the other hand, provided that the field q and its divergence are "smooth" in the sense that they belong to PH m+1 (R) for some integer m ≥ 0, using the RTN [m] finite element one can recover interpolation estimates in O(h m+1 ), cf. [BoBF13, §2.5.5]. For meshes made of affine elements such as tetrahedra or parallelepipeds, the approximation estimate (5.9-top) does not require the term with the divergence (see, e.g. [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF], §2.5.1).

A Priori Error Estimates

Since we focus on the low-regularity case, we choose the RTN [0] finite element, i.e.

X h = Q 0 h × L 0 h .
If the solution is "smooth", one can increase the order of the RTN finite element. This will be used in particular in §5.1.4 for the study of the error on the eigenvalues. According to first Strang's Lemma [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF] and because (1 + c (η h ) -1 ) 1, the error reads:

ζ -ζ h X inf ξ h ∈X h ζ -ξ h X .
(5.10)

Theorem 5.3. Under the assumptions of Proposition 2.5, it holds, with r max < 1/2:

∀µ ∈]0, r max [, ∀S f ∈ H µ (R), p -p h div ,R + φ -φ h 0,R h µ S f µ,R .
(5.11)

Remark 5.4. In particular, for "smooth data" S f , i.e. S f ∈ H rmax (R), one expects a convergence rate at least in h rmax-η for η > 0 arbitrary small: by a slight abuse of notation there and in the sequel, we shall write h rmax . Also, the previous analysis can be extended to the case where r max is in [1/2, 1] and µ < r max (or µ ≤ 1 if r max = 1). Furthermore, for a "smooth" solution, one may recover a convergence rate like O(h m+1 ) for an RTN [m] discretization of order m ≥ 0.

Proof. Choosing ξ h = (p 0 R , π 0 φ) ∈ X h , then thanks to the a priori estimates (5.4) and (5.9), it follows that:

ζ -ξ h 2 X = p -p 0 R 2 div ,R + φ -π 0 φ 2 0,R h 2µ (|p| 2 H µ (R) + div p 2 H µ (R) ) + h 2 φ 2 H 1 (R) h 2µ S f 2 µ,R .

Aubin-Nitsche-type Estimates

To derive improved estimates on the error φ-φ h 0,R in X h = Q 0 h ×L 0 h , we shall rely on the illuminating work of Falk-Osborn [START_REF] Falk | Error Estimates for Mixed Methods[END_REF]. Interestingly, one can obtain an improvement of the convergence rate, contrary to the case where the solution is 'smooth". From the previous analysis, for all µ < r max , we already have the estimate (5.11).

Lemma 5.5. Let (p, φ) (resp. (p h , φ h )) the solution of continuous (resp. discrete) variational problem (2.11) (resp. (5.2)). For all (q h , ψ h ) in X h , it holds:

a(p -p h , q h ) + b(q h , φ -φ h ) = 0, (5.12) b(p -p h , ψ h ) + t(φ -φ h , ψ h ) = 0.
(5.13)

Proof. Let (q h , ψ h ) be in X h . The subtraction of (2.11) from (5.2), with (q, ψ) = (q h , ψ h ) in the former, gives

a(p -p h , q h ) + b(q h , φ -φ h ) + b(p -p h , ψ h ) + t(φ -φ h , ψ h ) = 0
We obtain the first equality (5.12) (resp. the second equality (5.13)) with ψ h = 0 (resp. q h = 0).

Before improving the estimate, we need to introduce the adjoint problem: For d ∈ L, find (y d , η d ) ∈ X such that ∀(q, ψ) ∈ X:

a(q, y d ) + b(q, η d ) + b(y d , ψ) + t(ψ, η d ) = (ψ, d) 0,R .
(5.14)

Theorem 5.6. Under the assumptions of Proposition 2.5, it holds, with r max < 1/2:

∀µ ∈]0, r max [, ∀S f ∈ H µ (R), φ -φ h 0,R h 2µ S f µ,R .
(5.15)

Proof. Adapting the methodology of [START_REF] Falk | Error Estimates for Mixed Methods[END_REF] and by using (0, φφ h ) as a test function in the adjoint problem (5.14), we remark:

φ -φ h 0,R = sup d∈L\{0} b(y d , φ -φ h ) + t(φ -φ h , η d ) d 0,R .
(5.16)

We now look for an upper bound of the supremum in (5.16). We find that the numerator is successively equal to:

b(y d -(y d ) 0 R , φ -φ h ) + b((y d ) 0 R , φ -φ h ) + t(φ -φ h , η d ) ; using (5.7), for any ψ * h , ψ h in L h : b(y d -(y d ) 0 R , φ -ψ * h ) + b((y d ) 0 R , φ -φ h ) + t(φ -φ h , η d -ψ h ) + t(φ -φ h , ψ h ) ; using (5.12) with q h = (y d ) 0 R : b(y d -(y d ) 0 R , φ -ψ * h ) -a(p -p h , (y d ) 0 R ) + t(φ -φ h , η d -ψ h ) + t(φ -φ h , ψ h ) ; now we use (5.13) with ψ h = ψ h : b(y d -(y d ) 0 R , φ -ψ * h ) -a(p -p h , (y d ) 0 R ) + t(φ -φ h , η d -ψ h ) -b(p -p h , ψ h ) ; we add (5.14) with (p -p h , 0) as a test function: b(y d -(y d ) 0 R , φ -ψ * h ) + a(p -p h , y d -(y d ) 0 R ) +t(φ -φ h , η d -ψ h ) + b(p -p h , η d -ψ h ).
(5.17)

All terms( * ) in the previous relation can be bounded with an h-dependent term:

inf ψ * h ∈L h |b(y d -(y d ) 0 R , φ -ψ * h )| div (y d -(y d ) 0 R ) 0,R inf ψ * h ∈L h φ -ψ * h 0,R div y d 0,R h φ 1,R h S f µ,R d 0,R ; |a(p -p h , y d -(y d ) 0 R )| p -p h 0,R y d -(y d ) 0 R 0,R h µ S f µ,R (h µ |y d | µ,R + h div y d 0,R ) h 2µ S f µ,R d 0,R .
The last two terms in (5.17) are considered together.

inf

ψ h ∈L h |b(p -p h , η d -ψ h ) + t(φ -φ h , η d -ψ h )| ( div (p -p h ) 0,R + φ -φ h 0,R ) inf ψ h ∈L h η d -ψ h 0,R h µ S f µ,R inf ψ h ∈L h η d -ψ h 0,R h µ S f µ,R h η d 1,R h µ+1 S f µ,R d 0,R .
Thus, for low-regularity solutions (µ < 1/2), we conclude that it holds:

φ -φ h 0,R max(h, h 2µ , h µ+1 ) S f µ,R h 2µ S f µ,R .
Corollary 5.7. In the case of "smooth data" S f , i.e. S f ∈ H rmax (R), the error estimate gives: φφ h 0,R h 2rmax S f rmax,R .

Numerical Analysis of the Generalized Eigenvalue Problem

Let us focus on the approximation of the generalized eigenvalue problem (2.7) in our low-regularity setting, under the assumptions of Proposition 2.5, supplemented with

νΣ f ∈ PW 1,∞ (R).
The approximation of mixed eigenvalue problem have been studied in [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF] and references therein. However, the framework developed there can only be applied in our case when Σ r,0 vanishes in problem 2.14, that is when one solves:

Find (λ, p, φ) ∈ R × Q × V \{0} such that D -1 p + grad x φ = 0 in R; div (p) = λ -1 νΣ f φ in R.
More precisely, if one compares the solution of the neutron diffusion source problem (p, φ) to the discrete neutron diffusion source problem (p h , φ h ), this framework crucially relies on the fact that, for all ψ h in L h , it holds that b(p Rp h , ψ h ) = 0. Whereas in our case, one has by definition b

(p R -p h , ψ h ) = -t(φ -φ h , ψ h ) = -R Σ r,0 (φ -φ h )ψ h ,
and this term does not vanish in general because Σ r,0 = 0. Thus, we propose another approach to address this difficulty.

Let 0 ≤ µ < r max be given, we introduce an operator B µ associated to the source problem (2.11): given f ∈ H µ (R), we call B µ f = φ ∈ H 1 (R) the second component of the couple (p, φ) that solves (2.11) with source S f = νΣ f f . Since νΣ f belongs to PW 1,∞ (R), it holds S f µ,R f µ,R because µ < 1/2. Hence, B µ is a bounded operator from H µ (R) to itself:

B µ f µ,R B µ f 1,R = φ 1,R S f 0,R S f µ,R f µ,R ;
we write B µ ∈ L(H µ (R)) for short. In addition, since the second component of the solution actually belongs to H 1 (R) with continuous dependence (

φ 1,R f µ,R ), it
Let F h be the restriction of E h to W . One has the following simple results as a consequence of the gap property.

Lemma 5.15. For h small enough, F h is a bijection from W to W h . Moreover

∀ϕ ∈ W, (νΣ f ) 1 2 (ϕ -F h ϕ) 0,R h ω ϕ W .
(5.28)

Let S h = F -1 h E h -I ∈ L(L) for h small enough.
Lemma 5.16. For h small enough, W ⊂ ker(S h ) ; (S h ) h is uniformly bounded.

One can then prove an "orthogonality" result involving S h .

Proposition 5.17. For all f in L and ϕ h in W h , one has for h small enough

(νΣ f S h f, ϕ h ) 0,R = 0.
Proof. Let f be in L and ϕ h be in W h . We find:

(νΣ f S h f, ϕ h ) 0,R = (νΣ f (F -1 h E h f -f ), ϕ h ) 0,R = (νΣ f (F -1 h E h f -E h f ), ϕ h ) 0,R = (νΣ f (F -1 h E h f -F h F -1 h E h f ), ϕ h ) 0,R .
The second equality uses (5.27) with ϕ = f . One concludes by remarking that ψ =

F -1 h E h f ∈ W so (νΣ f (ψ -F h ψ), ϕ h ) 0,R = 0 using again (5.27), because F h ψ = E h ψ.
To obtain an optimal rate of convergence we restrict the operators B µ and B h µ to the eigenspace W . We denote finally by Bµ and Bh µ the operators, from W to itself,

Bµ = B µ | W and Bh µ = F -1 h B h µ F h . Let us estimate Bµ -Bh µ L(W ) = sup ϕ,ϕ ∈W \{0} |(ϕ, ( Bµ -Bh µ )ϕ ) W | ϕ W ϕ W .
Theorem 5.18. Let ω be as in definition 5.11. Then for h small enough, the following estimate holds true Bµ -Bh

µ L(W ) h 2ω .
(5.29)

Proof. Using the definition of F h , Lemma 5.14 and finally Lemma 5.16, one checks that for all ϕ ∈ W :

( Bµ -Bh µ )ϕ = B µ ϕ -F -1 h B h µ F h ϕ = B µ ϕ -F -1 h B h µ E h ϕ = B µ ϕ -F -1 h E h B h µ ϕ = (B µ -B h µ )ϕ + B h µ ϕ -F -1 h E h B h µ ϕ + S h B µ ϕ = (B µ -B h µ )ϕ + S h (B µ -B h µ )ϕ .
(5.30)

Hence, given ϕ, ϕ ∈ W , we can bound |(ϕ, ( Bµ -Bh µ )ϕ ) W | = |(νΣ f ϕ, ( Bµ -Bh µ )ϕ ) 0,R | by |(νΣ f ϕ, (B µ -B h µ )ϕ ) 0,R | + |(νΣ f ϕ, S h (B µ -B h µ )ϕ ) 0,R |.
Condition 5.21. For all energy groups 1 ≤ g, g ≤ G, g = g and for all 0 ≤ n ≤ N ,it stands:

               (Σ g r,n , Σ g →g s,n , νΣ g f ) ∈ PW 1,∞ (R) × L ∞ (R) × L ∞ (R), ∃(Σ r,n ) * , (Σ r,n ) * > 0, 0 < (Σ r,n ) * ≤ t n Σ g r,n ≤ (Σ r,n ) * a.e. in R, ∃0 < ε < 1 G -1 , |Σ g →g s,n | ≤ εΣ g r,n a.e. in R 0 ≤ νΣ g f a.
e. in R, ∃g, g s.t. χ g,g νΣ g f = 0.

(5.31a)

(5.31b)

(5.31c)

(5.31d)

Moreover, it stands:

   α s,o α r,o (G -1) < 1, α s,e (G -1) < 1 1 + α r,e
.

The coefficients are supposed to satisfy condition 5.21. We recall from proposition 1.14 that the matrices T o and T -1 e are bounded and have a positive property. Moreover, we suppose that hypothesis 3.8 holds. We recall here this hypothesis.

Hypothesis 5.22. There exists r max ∈]0, 1] such that for all source terms S f ∈ L, the flux solution of problem 3.7 φ ∈ V belongs to 0≤r<rmax PH 1+r (r max < 1) or PH 2 (r max = 1) with continuous dependence:

∀r ∈ [0, r max [, φ PH 1+r S f L (r max < 1) or φ PH 2 S f L (r max = 1).
We recall from problem 3.15 that the multigroup SP N source problem reads Problem 5.23. Find ζ = (p, φ) in X such that for all ξ = (q, ψ) ∈ X:

c s (ζ, ξ) = a s (p, q) + b s (q, φ) + b s (p, φ) + t s (φ, ψ) = (S f , ψ) L = f s (ξ), (5.32)
where a s , b s , t s and f s are define on page 49.

As we still are focused on the low-regularity case, when r max < 1/2, the discrete spaces for the approximation are:

Q h = (Q N h ) G ; L h = (L N h ) G ; abnd X h = Q h × L h . The discrete SP N variational source problem reads: Problem 5.24. Find ζ h ∈ X h such that for all ξ h ∈ X h c s (ζ h , ξ h ) = a s (p h , q h ) + b s (q h , φ h ) + b s (p h , φ h ) + t s (φ h , ψ h ) = (S f , ψ h ) L = f s (ξ h ). (5.33)
Thanks to proposition 1.14, one can show that the SP N variational problem 5.24 is wellposed.

Theorem 5.25. Under condition 5.21 and hypothesis 5.22, for any source term S f in L, there exists a unique solution ζ h in X h satisfying problem 5.24.

Proof. The proof is the same as the one for theorem 3.16 with (q h , ψ h ) = (-p h , 1 2 φ h + 1 2 T -1 e,h div p), where T -1 e,h is the matrix of the projection of all the components of T -1 e on L 0 h . Using estimation (5.4) on each component, one can prove that for all ψ ∈ L, it yields T -1 e,h ψ -T -1 e ψ L h ψ L .

Numerical Analysis of the Source Problem

For q in Q, its RTN [k] -interpolant q k R is the vector of all the RTN [k] -interpolant of its components, for all 1 ≤ g ≤ G, for n ∈ J o , (q k R ) g n = q g,k n,R . Lemma 3.3 in [BGNR06], estimates (5.9), can be extended to Q. Proposition 5.26. Let q be in H r , such that div q is in H s , 0 < r, s < r max , it stands:

q -q 0 R L h r |q| H r + h div q L , div (q -q 0 R ) 0,R h s |div q| H s .
(5.34)

Proof. Let q be in H r . q -q 0 R 2 L = G g=1 n∈Jo q g n -q g,0 n,R 2 L .
Using the first inequality in (5.9), one obtains:

q -q 0 R 2 L G g=1 n∈Jo (h r |q g n | H r + h div q g n L ) 2 .
For all a, b in R, (a + b) 2 ≤ 2(a 2 + b 2 ), therefore, it stands:

q -q 0 R 2 L h 2r G g=1 n∈Jo |q g n | 2 H r + h 2 G g=1 n∈Jo div q g n 2 L .
For all a, b in R positive, √ a + b ≤ √ a + √ b, therefore, we obtain:

q -q 0 R 2 L h r |q| H r + h div q 2 L .
We use the same methodology to obtain the second inequality in (5.34).

We define the operator π 0 from L to L 0 h = ((L 0 h ) N ) G by:

∀ψ ∈ L, π 0 ψ = ((π 0 ψ g n ) N n=0 ) G g=1 . According to [ErGu04, Proposition 1.135]: ∀ψ ∈ L, ψ -π 0 ψ L ψ L , ∀ψ ∈ PH 1 , ψ -π 0 ψ L h ψ PH 1 , ∀ψ ∈ P((W 1,∞ (R)) N ) G , ψ -π 0 ψ ((L ∞ (R)) N ) G h ψ P((W 1,∞ (R)) N ) G .
(5.35) Moreover, equations (5.7)-(5.8) can be extend for the multigroup SP N case:

Proposition 5.27. Let q be in Q, its RTN [k] -interpolant q k R in Q k h satisfies: ∀ψ h ∈ L k h , b s (q -q k R , ψ h ) = 0.
(5.36)

Moreover, it stands: div q 0 R = π 0 (div q).

(5.37)

Proof. Let q be in Q and ψ h be in L h . From the definition of b s (3.15) on page 49, one can remark that:

b s (q -q k R , ψ h ) = G g=1 n∈Jo b(q g n -q g,k n,R , ψ g h,n ).
From equation (5.7) on page 63, one obtains equation (5.36). For equation (5.37), we remark that: div q 0 R = ((div q g n ) n∈Jo ) G g=1 . Finally, using equation (5.8) on each component and using the definition of π 0 , one obtains equation (5.37).

A Priori Error Estimates

Theorem 5.28. Under condition 5.21 and hypothesis 5.22, it holds, with r max < 1/2:

∀µ ∈]0, r max [, ∀S f ∈ H µ , p -p h Q + φ -φ h L h µ S f H µ .
(5.38)

Proof. We use the same methodology as in the proof of theorem 5.3 for the mixed neutron diffusion equations, using inequalities (5.34) and inequalities (5.35). From

As, for the diffusion equations, we can obtain a better a priori error estimates thanks to an Aubin-Nitsche-type estimates. The matrices T o and T e are not symmetric, thus problem 3.2 is not symmetric, therefore its adjoint problem is different from the direct one.

The Adjoint Problem

The adjoint problem associated to the problem 3.15 reads:

Problem 5.29. For any d in L, find (y d , η d ) in X such that for all (q, ψ) ∈ X:

c s ((q, ψ), (y d , η d )) = a s (q, y d ) + b s (y d , ψ) + b s (q, η d ) + t s (ψ, η d ) = (ψ, d) L , (5.39)
where a s , b s and t s are define on page 49.

In the idea of the proof of theorem 2.19, one can prove that problem 5.29 is equivalent to the following problem.

(y d , η d ) X d L .
(5.41)

Proof. We recall from proposition 1.14 on page 27, that the matrices T o and T -1 e have a positive property and are bounded. Therefore, T T o and T T -1 e have a positive property and are bounded, which ensure the inf-sup condition for the adjoint problem 5.29. Moreover, in order to derive a uniform rate of convergence, we suppose that a property similar to hypothesis 5.22 holds (with the same regularity exponent r max for simplicity).

Hypothesis 5.32. We suppose that there exists r max ∈]0, 1] such that for all source terms d ∈ L, the solution of problem 5.30 η d ∈ V belongs to 0≤r<rmax P((H 1+r (R)) N ) G (r max < 1) or P((H 2 (R)) N ) G (r max = 1) with continuous dependence.

The discrete counterpart of the adjoint variational problem 5.29 reads: Problem 5.33. For any d in L, find (y d,h , η d,h ) in X h such that for all (q h , ψ h ) ∈ X h :

c s ((q h , ψ h ), (y d,h , η d,h )) = (ψ h , d) L .
(5.42)

As, under condition 5.21, the matrices T T o and T T -1 e have a positive property and are bounded, the discrete adjoint problem 5.33 is well-posed.

Theorem 5.34. Under condition 5.21 and hypothesis 5.32, for any source term d in L, there exists a unique solution (y d,h , η d,h ) in X h satisfying problem 5.24.

And we have the following a priori error estimate.

Theorem 5.35. Under condition 5.21 and hypothesis 5.32, it holds, with r max < 1/2:

∀µ ∈]0, r max [, ∀d ∈ H µ , y d -y d,h Q + η d -η d,h L h µ d H µ .
(5.43)

Aubin-Nitsche-type Error Estimates

To derive an improved error estimates, we work in the same way as for the mixed neutron diffusion equations in § 5.1.3.

Lemma 5.36. Let (p, φ) (resp. (p h , φ h )) the solution of continuous (resp. discrete) variational problem (3.19) (resp. (5.33)). For all (q h , ψ h ) in X h , it holds:

a s (p -p h , q h ) + b s (q h , φ -φ h ) = 0, (5.44) b s (p -p h , ψ h ) + t s (φ -φ h , ψ h ) = 0.
(5.45)

Theorem 5.37. Under condition 5.21 and hypothesis 5.22, it holds, with r max < 1/2:

∀µ ∈]0, r max [, ∀S f ∈ H µ , φ -φ h L h 2µ S f H µ .
(5.46)

Proof. Using (0, φφ h ) as a test function in the adjoint problem 5.29, we remark:

φ -φ h L = sup d∈L\{0} b s (y d , φ -φ h ) + t s (φ -φ h , η d ) d L .
(5.47)

We adapt the methodology of the proof of theorem 5.6. The numerator of the supremum in (5.47) is successively equal to:

b s (y d -(y d ) 0 R , φ -φ h ) + b s ((y d ) 0 R , φ -φ h ) + t s (φ -φ h , η d ); using equation (5.36), for any ψ * h , ψ h in L h : b s (y d -(y d ) 0 R , φ -ψ * h ) + b s ((y d ) 0 R , φ -φ h ) + t s (φ -φ h , η d -ψ h ) + t s (φ -φ h , ψ h ); using (5.44) with q h = (y d ) 0 R : b s (y d -(y d ) 0 R , φ -ψ * h ) -a s (p -p h , (y d ) 0 R ) + t s (φ -φ h , η d -ψ h ) + t s (φ -φ h , ψ h ); using (5.45) with ψ h = ψ h : b s (y d -(y d ) 0 R , φ -ψ * h ) -a s (p -p h , (y d ) 0 R ) + t s (φ -φ h , η d -ψ h ) -b s (p -p h , ψ h )
; we add equation (5.39) with (pp h , 0) as a test function:

b s (y d -(y d ) 0 R , φ -ψ * h ) + a s (p -p h , y d -(y d ) 0 R ) +t s (φ -φ h , η d -ψ h ) + b s (p -p h , y d ) -ψ h );
(5.48)

Now, we bound each term in the previous equality. The first term is bounded as:

inf ψ * h ∈L h |b s (y d -(y d ) 0 R , φ -ψ * h )| div (y d -(y d ) 0 R ) L inf ψ * h ∈L h φ -ψ * h L div y d L h φ V h S f H µ d L ;
where we use (5.37) and (5.35) for div (y d -(y d ) 0 R ) L and (5.35) for inf ψ * h ∈L h φψ * h L . The second term in (5.48) is bounded as:

|a s (p -p h , y d -(y d ) 0 R )| p -p h L y d -(y d ) 0 R L h µ S f H µ h r |y d | H r + h div y d L h 2µ S f H µ d L ;
where we use (5.38) for pp h L and (5.34) for y d -(y d ) 0 R L . The second term in (5.48) are bounded together as:

inf ψ h ∈L h |t s (φ -φ h ,η d -ψ h ) + b s (p -p h , y d ) -ψ h )| div (p -p h ) L + φ -φ h L inf ψ h ∈L h y d ) -ψ h L h µ S f H µ h η d V h µ+1 S f H µ d L

Therefore, it holds:

φφ h L min(h, h 2µ , h µ+1 ) S f H µ Thus for low regularity solution (µ < 1/2), we have min(h, h 2µ , h µ+1 ) h 2µ .

Corollary 5.38. In the case of "smooth data" S f , i.e. S f ∈ H rmax , the error estimate gives:

φ -φ h L h 2rmax S f H rmax .
We remark that ϕ 2 W = (M f ϕ, ϕ) L is a norm over W (see the definition of • W in § 5.1). However, this norm does not derive from a inner product since M f is not symmetric. For this reason, we cannot use the same methodology as in § 5.1.4. Proposition 5.42. Let ω be as in definition 5.41. For every ϕ in W , the following inequalities hold:

(B µ -B h µ )ϕ L h ω ϕ W .
Proof. This inequality comes from the first Strang's lemma. The method is the same as for theorem 5.3. Here, we use the equivalence of all norms on W to state the result. Chapter 6 Domain Decomposition Method

The Neutron Diffusion Equations

We continue by considering the neutron diffusion problem under its mixed form using a Domain Decomposition Method (DDM). The diffusion problem with low-regularity solution in a mixed, multi-domain form has been analyzed in [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF]: to solve the problem, the authors introduce a method called the DD+L 2 -jumps. This method can handle globally non-conforming triangulation with locally conforming triangulations. As we need a Lagrange multiplier, it can be considered as a generalization of some hybrid finite element method [START_REF] Roberts | Mixed and hybrid finite element methods[END_REF], or as a mortar finite element method [START_REF] Bernardi | Domain decomposition by the mortar element method[END_REF][START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF]. Mortar finite elements were introduced in the late eighties by Y. Maday et al to couple spectral and finite element methods [START_REF] Bernardi | Couplage de methodes spectrale et d'elements finis: Premiers resultats d'approximation[END_REF], and to handle complex geometries with spectral discretization [START_REF] Maday | Nonconforming mortar element methods: Application to spectral discretizations[END_REF]. These method consists in adding an interface condition between the subdomains directly in the discrete functional space. The interface conditions can be imposed in a constrained, piecewise regular space [START_REF] Bernardi | Domain decomposition by the mortar element method[END_REF][START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF]. In that case, the support of the constraint is called the non-mortar side, and is chosen to be the finest triangulation on the interface. This condition brings also a Lagrange multiplier in the discrete trace space of the discrete solution, which ensures the weak continuity of the solution between the subdomains [ Ben 99]. We refer to for an overview of these methods [Wohl01]. In the first subsection 6.1.1, we first define some notations and spaces, and then we recall some results that ensure the well-posedness of the discrete DDM. The numerical analysis of the Domain Decomposition Method is carried out in the next subsection §6.1.2. This section comes from the article [START_REF] Ciarlet | Numerical Analysis of the Mixed Finite Element Method for the Neutron Diffusion Eigenproblem with Heterogeneous Coefficients[END_REF].

Setting of the Domain Decomposition Method

Let us consider a partition { R i } 1≤i≤ N of R which can be independent from the physical partition of the materials in R (see e.g. [START_REF] Jamelot | Fast Non-Overlapping Schwarz Domain Decomposition Methods for Solving the Neutron Diffusion Equation[END_REF][START_REF] Brenner | A multigrid algorithm for the lowest-order Raviart-Thomas mixed trangular finite element method[END_REF][START_REF] Braess | A posteriori error estimators for the Raviart-Thomas element[END_REF]). In other words, it can happen that { R i } 1≤i≤ N = {R i } 1≤i≤N . We denote by Γ ij the interface between two subdomains R i and R j , for i = j:

if the Hausdorff dimension of R i ∩ R j is d -1, then Γ ij = int( R i ∩ R j ); otherwise, Γ ij = ∅. By construction, Γ ij = Γ ji .
We define the interface Γ S , respectively the wirebasket ∂Γ W by

Γ S = N i=1 N j=i+1 Γ ij , ∂Γ W = N i=1 N j=i+1 ∂Γ ij .
It is stressed that the resulting interface Γ S needs not necessarily coincide with the physical interface between cells. When d = 2, the wirebasket consists of isolated crosspoints. When d = 3, the wirebasket consists of open edges and crosspoints. For a field v defined over R, we shall use the notation v i = v | R i , for 1 ≤ i ≤ N . Let us define the function space with zero Dirichlet boundary condition:

PH 1 0 (R) = ψ ∈ L 2 (R) | ψ i ∈ H 1 ( R i ), ψ |∂ R i \Γ S = 0, 1 ≤ i ≤ N . When Γ ij = ∅, let H 1/2 Γ ij be the set of H 1/2 (Γ ij ) functions whose continuation by 0 to ∂R i belongs to H 1/2 (∂R i ). On can prove that H 1/2 Γ ij = H 1/2
Γ ji . We also introduce the space of piecewise H(div ) vector-valued functions:

PH(div , R) = q ∈ L 2 (R) | q i ∈ H(div , R i ), 1 ≤ i ≤ N , ||q|| PH(div ,R) = i q i 2 H(div , R i ) 1/2 . For p ∈ PH(div , R), let us set [p • n] ij := k=i,j p k • n k|Γ ij the jump of the normal component of p on Γ ij when Γ ij = ∅. [p • n] ij is well defined in (H 1/2 Γ ij ) the dual space of H 1/2
Γ ij (see e.g. [START_REF] Fernandes | Magnetostatic and electrostatic problems in inhomogeneous anisotropic media with irregular boundary and mixed boundary conditions[END_REF]). The global jump [p • n] of the normal component on the interface is defined by:

[p • n] |Γ ij := [p • n] ij , for 1 ≤ i, j ≤ N .
By definition, it holds

[p • n] ∈ i<j (H 1/2 Γ ij )
. We recall that for p ∈ H(div , R), the global jump vanishes: [p • n] = 0 (see e.g. [CiJK17, Lemma 1]). We introduce finally the following Hilbert spaces:

M = ψ S ∈ i<j L 2 (Γ ij ) , ψ S M = i<j ψ S 2 0,Γ ij 1/2 ; H 1/2 -(Γ S ) = ψ S ∈ M | ψ S|Γ ij ∈ H 1/2 (Γ ij ), ∀i < j , with graph norm ; Q = q ∈ PH(div , R) | [q • n] ∈ M , q Q = ||q|| 2 PH(div ,R) + ||[q • n]|| 2 M 1/2 ; X = ξ := (q, ψ) ∈ Q × L 2 (R) , ξ X := q 2 Q + ψ 2 0,R 1/2 ; W = w := (ξ, ψ S ) ∈ X × M , w W := ξ 2 X + ψ S 2 M 1/2
.

By construction, one has M ⊂ i<j (H 1/2 Γ ij ) . We will next define a variational formulation which is conforming in Q × L 2 (R). The mixed form of the neutron diffusion problem (2.8) is now given by (see [CiJK17, §3.2]):

Find (p, φ, φ S ) ∈ Q × PH 1 0 (R) × M such that:          -D -1 i p i -grad φ i = 0 in R i , for 1 ≤ i ≤ N , div p i + Σ a,i φ i = S f,i in R i , for 1 ≤ i ≤ N , φ i = φ S on ∂ R i ∩ Γ S , for 1 ≤ i ≤ N , [p • n] = 0 on Γ S .
(6.1)

To solve this problem, we are looking for a solution

((p, φ), φ S ) in W. Find ((p, φ), φ S ) ∈ W, such that ∀((q, ψ), ψ S ) ∈ W: R -D -1 p • q + φ div q + ψ div p + Σ a φ ψ + Γ S [p • n] ψ S - Γ S [q • n] φ S = R S f ψ. (6.2)
In (6.1)-(6.2), φ S , ψ S play the role of Lagrange multipliers, with M the space of those Lagrange multipliers. To be mathematically precise, we should be integrating on ∪ i<j Γ ij instead of Γ S . We make this slight abuse of notations from now on. This approach is called the DD+L 2 -jumps method.

From now on, we use the notations:

• u = (ζ, φ S ), ζ = (p, φ), p = (p i ) 1≤i≤ N and φ = (φ i ) 1≤i≤ N ;
• w = (ξ, ψ S ), ξ = (q, ψ), q = (q i ) 1≤i≤ N and ψ = (ψ i ) 1≤i≤ N ; and we define the bilinear forms:

S :    W × W → R (u, w) → Γ S [p • n] ψ S , (6.3)
and:

c S : W × W → R (u, w) → c(ζ, ξ) + S (u, w) -S (w, u) . (6.4)
We consider the linear form:

f S : W → R w → f (ξ) . (6.5)
Above, we extended the definition (2.15) (resp. (2.16)) of the form c (resp. f ), to elements of X × X (resp. X). We may rewrite the variational formulation (6.2) as:

Find u ∈ W such that ∀w ∈ W: c S (u, w) = f S (w). (6.6)
We recall that c S satisfies an inf-sup condition, so the variational problem is well-posed (see [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF]§4]), and that, under the assumptions of proposition 2.5, the global jump of p vanishes:

[p • n] = 0 in M (see [CiJK17, Lemma 1]).
We study abstract, conforming, discretization of the variational formulation (6.6) as it is done in [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF]§5]. To that aim, we introduce discrete, finite-dimensional, spaces indexed by a (small) parameter h as follows:

Q i,h ⊂ H(div , R i ) and L i,h ⊂ L 2 ( R i ), for 1 ≤ i ≤ N .
We impose the following requirements, for all 1 ≤ i ≤ N :

• q i,h • n |∂ R i ∈ L 2 (∂ R i ) for all h > 0, for all q i,h ∈ Q i,h ; • div Q i,h ⊂ L i,h for all h > 0 ;
• (Q i,h ) h and (L i,h ) h satisfy the approximability property (5.1) in R i .

Then, let

Q h = 1≤i≤ N Q i,h and L h = 1≤i≤ N L i,h .
In particular, the discretization

Q h × L h is globally conforming in Q × L 2 (R). We endow Q h with the norm • Q , while L h is endowed with • 0,R .
We then define T i,h as the space of the normal traces of vectors of

Q i,h on ∂ R i ∩ Γ S : T i,h := q i,h ∈ L 2 (∂ R i ∩ Γ S ) | ∃ q i,h ∈ Q i,h , q i,h = q i,h • n i|∂ R i ∩Γ S . (6.7)
Classically, several situations can occur on a given interface Γ ij , 1 ≤ i, j ≤ N :

1. non-nested meshes:

T i,h|Γ ij ⊂ T j,h|Γ ij and T j,h|Γ ij ⊂ T i,h|Γ ij ; 2. nested meshes: T i,h|Γ ij ⊂ T j,h|Γ ij or T j,h|Γ ij ⊂ T i,h|Γ ij ;
3. matching meshes: nested meshes with T i,h|Γ ij = T j,h|Γ ij .

Usually, the term nested meshes is used to describe a family of successively refined meshes.

In this paper, we will use this expression to express that on all interfaces Γ ij , case (2) described above holds. As an illustration, see the interfaces between the subdomains, in figure 6.1, for nested non-matching mesh (left) and non-nested mesh (right).

Nested Mesh Non-nested Mesh Let us denote by M h ⊂ M the discrete space of the Lagrange multipliers. We assume that M h includes the subspace M 0 h of piecewise constant fields. We introduce the discrete projection operators [CiJK17, §5] from the spaces of normal traces T i,h to M h , and vice versa, which are defined by:

∀q i,h ∈ T i,h , ∀ψ S,h ∈ M h        ∂ R i ∩Γ S (Π i (q i,h ) -q i,h ) ψ S,h = 0 ∂ R i ∩Γ S (π i (ψ S,h ) -ψ S,h ) q i,h = 0 . (6.8)
As the operators Π i and π i are orthogonal projections, they are continuous, with a continuity modulus equal to 1. We also introduce the orthogonal projection operator Π 0 S : M → M 0 h . According to [ErGu04, Proposition 1.135], if we denote by h S the meshsize on Γ S :

∀ψ S ∈ H 1/2 -(Γ S ), ψ S -Π 0 S (ψ S ) M h 1/2 S ψ S H 1/2 -(Γ S ) .
(6.9)

Next, let p h ∈ Q h . We define the discrete jump of the normal component of p h on the interface

Γ ij as [p h • n] h,ij := l=i,j Π l (p l,h • n l|Γ ij ). The discrete global jump of the normal component, [p h • n] h ∈ M h , is defined by: [p h • n] h|Γ ij := [p h • n] h,ij , for 1 ≤ i, j ≤ N .
We finally define:

X h = ξ h := (q h , ψ h ) ∈ Q h × L h , endowed with • X , W h = w h := (ξ h , ψ S,h ) ∈ X h × M h , endowed with • W .
In the DD+L 2 -jumps setting, the conforming discretization of the variational formulation (6.6) reads:

Find u h ∈ W h such that ∀w h ∈ W h , c S (u h , w h ) = f S (w h ). (6.10)
It is shown in [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF]§5] that c S verifies a discrete inf-sup condition if the following conditions hold:

∃β h > 0, ∀q h ∈ Q h , Γ S [q h • n] h [q h • n] ≥ β h Γ S [q h • n] 2 (6.11) and ∃γ h > 0, ∀ψ S,h ∈ M h , N i=1 N j=i+1 Γ ij π i (ψ S,h ) 2 + π j (ψ S,h ) 2 ≥ γ h ψ S,h 2 M , (6.12)
Moreover, if β h and γ h can be chosen independently of h, the form c S satisfies a udisc. For instance, conditions (6.11)-(6.12) are uniformly fulfilled when M h is chosen as

M h = N i=1
T i,h . (6.13) Last, under (6.11), one easily checks that [p h • n] = 0. In other words:

p h ∈ H(div , R) ∩ Q h . (6.14)
For the DDM, we define

Q h = H(div , R) ∩ Q h .

Numerical Analysis of the Domain Decomposition Method

To carry out the numerical analysis in the low-regularity case, we first introduce a suitable discretization of the DD problem, and then we carry out the numerical analysis on this discretization. Again, if one chooses another discretization that fulfills those properties detailed in the previous section, one may recover similar convergence results.

Discretization

We consider (6.10) where the RTN finite element is used on each subdomain with a conforming mesh, or triangulation. For 1 ≤ i ≤ N , let h i denote the local meshsize in R i , and h = max i h i the global meshsize. Let us denote by k i ≥ 0 the order of the discretization in R i , and k = min i k i , the minimal order of the RTN finite element. The local RTN finite element subspace of

H(div , R i ) × L 2 ( R i ) is defined as Q k i i,h i × L k i i,h i . With this choice, we have div Q k i i,h i ⊂ L k i i,h i as required: local consistency is ensured. Now, if we set Q k h = 1≤i≤ N Q k i i,h i and L k h = 1≤i≤ N L k i i,h i , we have q i,h • n |∂ R i ∈ L 2 (∂ R i ) for all q i,h ∈ Q k i i,h i , hence it follows that Q k h ⊂ Q: the discretization Q k h × L k h is globally conforming in Q × L 2 (R).
For the reader's convenience, we omit the superscript k i in the analysis below. Finally, we choose M h so that on the one hand (6.11)-(6.12) hold uniformly, and on the other hand it holds h S h: we refer to [CiJK17, §5.2] for an extended discussion on suitable choices. According to the first Strang's Lemma [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF] and because c S verifies a udisc, the error reads:

u -u h W inf w h ∈W h u -w h W . (6.15)
As a consequence lim h→0 uu h W = 0. This result holds for nested and non-nested meshes. We study below how to improve the bound on the error, how to derive an Aubin-Nitsche estimate, and finally how to prove convergence for the generalized eigenvalue problem, for nested meshes( * ). As previously, those results hold under hypothesis 2.4 (plus νΣ f ∈ PW 1,∞ (R) for the eigenproblem). We focus again on the low-regularity case.

A Priori Error Estimates

Let q ∈ H(div , R) ∩ PH µ (R), with 0 < µ. A global RTN interpolant of q is defined on every subdomain R i via its restriction q i , and denoted by qi,R for 1 ≤ i ≤ N .

Definition 6.1. Let q ∈ H(div , R) ∩ PH µ (R), with 0 < µ. The global RTN interpolant qR of q is defined by, for 1 ≤ i ≤ N :

qR | R i = qi,R .
Below, we also use the orthogonal projection operators π 0 : L 2 (R) → L 0 h (see §5.1.3) and Π 0 S : M → M 0 h (see § 6.1.1). One has the following result, whose proof is given in the Appendix. Lemma 6.2. Assume that the meshes are nested, non-matching, on the interface Γ f c , and that they are quasi-uniform on Γ f c . To fix ideas, we assume T c,h|Γ f c ⊂ T f,h|Γ f c with

T c,h|Γ f c = T f,h|Γ f c ( † ). Let q ∈ H(div , R) ∩ H µ (R) with 0 < µ < 1/2, it holds: [q R • n] 0,Γ f c h 1/2 f q f H(div , R f ) .
Theorem 6.3. Let the assumptions of proposition 2.5 hold, with r max < 1/2. One has for matching meshes:

∀µ ∈]0, r max [, ∀S f ∈ H µ (R), p -p h H(div ,R) + φ -φ h 0,R + φ S -φ S,h M h µ S f µ,R .
(6.16)

For nested, non-matching meshes, the result holds under the assumption that on an interface Γ ij where the meshes T i,h|Γ ij and T j,h|Γ ij are non-matching (T i,h|Γ ij = T j,h|Γ ij ), the families of triangulations of T i,h|Γ ij and T j,h|Γ ij are quasi-uniform.

Proof. We bound the different contributions in the right-hand side of (6.15) for some appropriately chosen discrete field w h . Recall that u = ((p, φ), φ S ).

Matching meshes. We know that [p•n] = 0. For matching meshes, one has also

[p R •n] = 0, so [(p -pR ) • n] = 0.
Starting from (6.15), the conclusion follows. Indeed, according to the a priori estimates (5.4), (5.9) and (6.9),

w h = (p R , π 0 φ, Π 0 S (φ S )) ∈ W h is such that u -w h 2 W = N i=1 p i -p i,R 2
H(div , R i ) + φ -π 0 φ 2 0,R + φ S -Π 0 S (φ S ) 2 M h 2µ (|p| 2 µ,R + div p 2 µ,R ) + h 2 φ 2 PH 1 (R) + h S φ S 2 H 1/2 -(Γ S ) h 2µ S f 2 µ,R .
Hence we conclude that for matching meshes it holds:

u -u h W h µ S f µ,R .
(6.17)

Nested meshes. In this case, [p R • n] = 0 in general. Nonetheless, one can use the result of Lemma 6.2, to find that

[(p -pR ) • n] M h 1/2 p H(div ,R) ,
provided that the meshes are quasi-uniform on the part of the interface where they are non-matching. One concludes that the estimate (6.17) still holds for nested meshes under this condition.

Conclusion. Noting that it always holds [p • n] = [p h

• n] = 0 (cf. (6.14)), developing the norm uu h W , one concludes:

p -p h H(div ,R) + φ -φ h 0,R + φ S -φ S,h M h µ S f µ,R .
In other words, we have the a priori error estimate (6.16).

As in section 5.1, for "smooth data" S f , i.e. S f ∈ H rmax (R), one expects a convergence rate at least in h rmax .

Aubin-Nitsche-type Estimates

To derive improved estimates on the error φ-φ h 0,R , we adapt the calculations of § 5.1.3 to the DDM. Recall that Q h = Q h ∩ H(div , R). We already know that when conditions (6.11)-(6.12) hold, the solution

((p h , φ h ), φ S,h ) ∈ X h × M h of (6.10) (discrete DDM) is such that (p h , φ h ) ∈ X h , since p h ∈ Q h .
Then restricting the test-fields in (6.10) to elements of X h × M h we observe that (p h , φ h ) satisfies (5.3) too (discrete mixed problem in § 5), because all interface terms vanish. Hence, to estimate φ-φ h 0,R in the DDM, we explicitly consider that the discrete fields (p h , φ h ) are also the solution to the variational formulation (5.3). Let us begin by a technical result, whose proof is given in the Appendix.

Lemma 6.4. Let the assumptions of Lemma 6.2 hold. Let q ∈ H(div , R) ∩ H µ (R) with 0 < µ < 1/2, and define

δq f c ∈ Q f,h by δq f c • n |Γ f c = (q c,R • n -qf,R • n) |Γ f c and zero extension in R f \ Γ f c . It holds δq f c H(div , R f ) h µ q f µ, R f + div q f 0, R f .
Theorem 6.5. Under the assumptions of Theorem 6.3 with r max < 1/2, one has for nested meshes:

∀µ ∈]0, r max [, ∀S f ∈ H µ (R), φ -φ h 0,R h 2µ S f µ,R . (6.18)
Proof. Matching meshes. In this case, one can use the theory already developed in § 5.1.3, to conclude that (6.18) holds. Nested meshes. The difficulty for non-matching meshes is that one can not define the global RTN-interpolant of p directly. Instead it is defined via its subdomain interpolants ( p i,R ) 1≤i≤ N . Introduce, for 1 ≤ i ≤ N , I i as the set of indices j such that T j,h|Γ ij ⊂ T i,h|Γ ij (since we are dealing with nested meshes, it holds T j,h|Γ ij ⊂ T i,h|Γ ij or T i,h|Γ ij ⊂ T j,h|Γ ij ). We proceed as follows to obtain an H(div , R)-conforming approximant, i.e. an element of Q h . On all interfaces Γ ij , introduce

δp ij • n = pc,R • n |Γ ij -pf,R • n |Γ ij
where pf,R is the interpolant from the finer discretization on Γ ij , resp. pc,R is the interpolant from the coarser discretization on Γ ij . By construction, δp ij • n = 0 when T i,h|Γ ij = T j,h|Γ ij . Then δp ij • n is extended by zero in R i to define an element of Q i,h ; with a slight abuse of notation, we still denote the extension by δp ij . The H(div , R)-conforming approximant p R ∈ Q h is then defined subdomain by subdomain as

p i,R = p i,R + j∈I i δp ij for 1 ≤ i ≤ N . Indeed, [p R • n] Γ ij = 0 for 1 ≤ i, j ≤ N by direct inspection. It remains to evaluate p -p R 2 H(div ,R) = 1≤i≤ N p i -p i,R 2
H(div , R i ) , with p i -p i,R H(div , R i ) ≤ p i -p i,R H(div , R i ) + j∈I i δp ij H(div , R i ) for 1 ≤ i ≤ N .
Above, the fact that the index j belongs to I i implies that if δp ij = 0, then the finer discretization on Γ ij automatically originates from R i . To evaluate δp ij H(div , R i ) , one uses the results of Lemma 6.4 to find

δp ij H(div , R i ) h µ p i µ, R i + div p i 0, R i .

About Non-nested Meshes

We recall that, for general non-nested meshes, one has convergence without explicit convergence rate, as soon as (6.11)-(6.12) hold uniformly. In the most general case however, it seems difficult to obtain a convergence error that depends explicitly on h. On the other hand, let us consider the case where the meshes are non-nested, with some structure. By structure, it is understood that the non-nestedness can be described by a finite number of configurations (e.g. 3-face mesh vs. 5-face mesh, etc.) that are reproduced at smaller and smaller scales when the meshsize diminishes. We note first that a result similar to Lemma 6.2 can be recovered. Going back to the reference configurations (by assumption there are a finite number of them) and taking the supremum in the upper bounds among all these configurations, we infer from

(C.6) that [q R • n] 0,Γ f c h c|Γ f c q f,h 0,Γ f c , i.e.
one can conclude the proof as before. As a consequence, an explicit convergence rate may be derived for the source problem as in Theorem 6.3. Then, one may proceed in a similar fashion to prove Lemma 6.4, so as to derive an Aubin-Nitsche estimate as in Theorem 6.5. Finally, because interface terms are absent in the analysis of the convergence rate of the eigenvalues (see in particular (6.20)-(6.21)), such estimates can also be proved for non-nested meshes, with some structure.

Here we summarize the results on the different types of meshes. In the case where the meshes are matching or nested non-matching one has convergence and furthermore, we provide a convergence rate (theorem 6.5). If the meshes are non-nested but there exists a structure (see above), one has convergence and one can derive a convergence rate along the same lines as in the proof of theorem 6.5. For the more general case, we only prove the abstract convergence of the discrete solution to the continuous one cf. (6.15).

Extension to the Multigroup SP N Transport Equations

In this section, we extend the results of the domain decomposition method proposed in § 6.1 to the multigroup SP N transport equations. We recall that G is the number of groups and N is the number of odd and even moments. We still suppose the same conditions as in § 5.2, we recall here this condition: Condition 6.11. For all energy groups 1 ≤ g, g ≤ G, g = g and for all 0 ≤ n ≤ N ,it stands:

               (Σ g r,n , Σ g →g s,n , νΣ g f ) ∈ PW 1,∞ (R) × L ∞ (R) × L ∞ (R), ∃(Σ r,n ) * , (Σ r,n ) * > 0, 0 < (Σ r,n ) * ≤ t n Σ g r,n ≤ (Σ r,n ) * a.e. in R, ∃0 < ε < 1 G -1 , |Σ g →g s,n | ≤ εΣ g r,n a.e. in R 0 ≤ νΣ g f a.e. in R, ∃g, g s.t. χ g,g νΣ g f = 0. (6.25a) (6.25b) (6.25c) (6.25d)
Moreover, it stands:

   α s,o α r,o (G -1) < 1, α s,e (G -1) < 1 1 + α r,e
.

We still consider the same partition { R i } 1≤i≤ Ñ of R as in § 6.1. We denote

M = (M N ) G ; Q = ( Q N ) G ; X = ( X N ) G ; and W = (W N ) G .
Moreover, we denote, for

1 ≤ i < j ≤ Ñ , L i = ((L 2 ( R i )) N ) G , L ij = ((L 2 (Γ ij )) N ) G , Q i = ((H(div , R i ) N ) G and H µ i = ((H µ ( R i ) N ) G
, for all µ > 0. We also introduce the operator:

• : ((R d ) N ) G × R d → (R N ) G (x, y) → x • y = ((x g n • y) N n=0 ) G g=1
The multigroup SP N transport problem with a source, problem 3.2, can be rewritten as: Problem 6.12. For a given source

S f ∈ L, find ((p, φ), φ S ) ∈ Q × V × M such that:                      T o p i + grad x Hφ i = 0 in R i ; T Hdiv p i + T e φ i = S f,i in R i ; Hφ i = φ S on ∂ R i ∩ Γ S ; [p • n] = 0 on Γ S .
(6.26)

The variational formulation of problem 6.12 reads: Problem 6.13. For a given source S f ∈ L, find ((p, φ), φ S ) ∈ W such that for all ((q, ψ), ψ S ) ∈ W:

- R T o p q + R φ • T Hdiv q + R ψ • T Hdiv p + R T e φ • ψ + Γ S [p • n] • ψ S - Γ S [q • n] • φ S = R S f • ψ (6.27)
From now on, we use the notations:

• u = (ζ, φ S ), ζ = (p, φ), p = (p i ) 1≤i≤ N and φ = (φ i ) 1≤i≤ N ;
• w = (ξ, ψ S ), ξ = (q, ψ), q = (q i ) 1≤i≤ N and ψ = (ψ i ) 1≤i≤ N ;

and we define the bilinear forms:

s,S :

   W × W → R (u, w) → Γ S [p • n] ψ S , (6.28) and:
c s,S :

W × W → R (u, w) → c s (ζ, ξ) + s,S (u, w) -s,S (w, u) . (6.29)
We consider the linear form:

f s,S : W → R w → f s (ξ) . (6.30)
Above, we extended the definition (3.17) (resp. (3.18)) of the form c s (resp. f s ), to elements of X × X (resp. X). We may rewrite the variational problem 6.13 as:

A Priori Error Estimates

Let q ∈ Q ∩ PH µ , with 0 < µ. Each component q i , 1 ≤ i ≤ Ñ , has a global RTN interpolant q i,R , see definition 6.1 on page 84. Definition 6.18. Let q ∈ Q ∩ PH µ , with 0 < µ. The global RTN interpolant qR of q is defined by, for 1 ≤ g ≤ G and for 0 ≤ n ≤ N :

(q R ) g n = qg n,R .
Below, we also use the orthogonal projection operators π 0 : L → L 0 h (see § 5.2.1) and Π 0 S : M → M 0 h . The projection Π 0 S correspond to the projection of all the components on M h by the projection Π 0 S (see § 6.1.1). One has the following result.

Lemma 6.19. Assume that the meshes are nested, non-matching, on the interface Γ f c , and that they are quasi-uniform on Γ f c . To fix ideas, we assume

T c,h|Γ f c ⊂ T f,h|Γ f c with T c,h|Γ f c = T f,h|Γ f c . Let q ∈ Q ∩ H µ with 0 < µ < 1/2, it holds: [q R • n] L f c h 1/2 f q f Q f . Proof. Let q ∈ Q ∩ H µ with 0 < µ < 1/2
, by definition of the norm on L f c , one has:

[q R • n] 2 L f c = G g=1 N n=0 [q g n,R • n] 2 0,Γ f c .
Therefore, using lemma 6.2 on each components, one can finish the proof in the same idea as in the proof of property 5.26.

Theorem 6.20. Under condition 6.11 and hypothesis 5.22 with r max < 1/2, one has for matching meshes:

∀µ ∈]0, r max [, ∀S f ∈ H µ , p -p h Q + φ -φ h L + φ S -φ S,h M h µ S f H µ . (6.34)
For nested, non-matching meshes, the result holds under the assumption that on an interface Γ ij where the meshes T i,h|Γ ij and T j,h|Γ ij are non-matching (T i,h|Γ ij = T j,h|Γ ij ), the families of triangulations of T i,h|Γ ij and T j,h|Γ ij are quasi-uniform.

Proof. The proof is the same as for theorem 6.3, using proposition 5.26, lemma 6.19 and inequalities (5.35).

As in section 5.2, for "smooth data" S f , i.e. S f ∈ H rmax , one expects a convergence rate at least in h rmax .

Aubin-Nitsche-type Estimates

To derive improved estimates on the error φφ h L , we adapt the calculations of § 5.2 to the DDM. Recall that Q h = Q h ∩ Q. We already know that when conditions (6.11)-(6.12) hold for each energy group and for each moment, the solution

((p h , φ h ), φ S,h ) ∈ X h × M h Part III

Implementation

This part is dedicated first to the implementation of the Raviart-Thomas finite element method on Cartesian meshes with the L 2 -jump domain decomposition method for the resolution of the multigroup SP N neutron transport equations. Secondly, we propose some adaptive methods to improve the quality of the solution. The Raviart-Thomas finite element method was already implemented in APOLLO3 R before this work [START_REF] Baudron | MINOS: A Simplified Pn Solver for Core Calculation[END_REF][START_REF] Baudron | Simplified p n Transport Core Calculations in the APOLLO3 System[END_REF]. An optimized Schwartz domain decomposition method is availlable fo parallel computation [JaBL12, JaCi13, JCBL14]. Some a posteriori error estimates are proposed in [Wang09, Lath11] for the resolution of the S N neutron transport equations with discontinuous Galerkin method in order to use a adaptive mesh refinement. In [START_REF] Ragusa | A simple Hessian-based 3D mesh adaptation technique with applications to the multigroup diffusion equations[END_REF][START_REF] Wang | Three-Dimensional h-Adaptivity for the Multigroup Neutron Diffusion Equations[END_REF], the authors propose an a posteriori error estimates for the multigroup diffusion equations under their primal setting. We propose here an a posteriori error estimate for the resolution diffusion equations under their mixed form. This part is composed of two chapters which is organized as follow:

• in chapter 7, in section 7.1, we describe the algorithm used in the neutronic platform APOLLO3 R to solve the multigroup SP N neutron transport;

• in section 7.2, we present a test case for the neutron diffusion equations where the L 2 -jump domain decomposition method is used with non-conforming triangulation;

• in section 7.3, we apply here L 2 -jump domain decomposition method is used with non-conforming triangulation for the resolution of the the multigroup SP N neutron transport on a PWR reactor.

• In chapter 8, in section 8.1, we derive an a posteriori error estimate for the neutron diffusion equations;

• in section 8.2, we present some possible improvements of the resolution of the multigroup SP N neutron transport.

Chapter 7

Numerical Applications

APOLLO3 R is a shared neutronic platform of CEA and EDF. It includes different deterministic solvers. Each solver is characterized by the angular discretization and the spatial discretization of the neutron transport equation. Our work is to integrate the L 2 -jump domain decomposition in the MINOS solver. The interest of this method is a way to parallelize the solver and make it able to use non-conforming triangulations. Let us mention that another numerical solver of the multigroup SP N neutron transport equations, named COCAGNE [START_REF] Calloo | COCAGNE: EDF new neutronic core code for ANDROMÈDE calculation chain[END_REF], has recently been developed at EDF. It is based on a domain decomposition method [START_REF] Lathuilière | Méthode de Décomposition de Domaine Pour Les Équations Du Transport Simplifié En Neutronique[END_REF] which is similar to the numerical implementation of the L 2 -jump domain decomposition. However, there is no mathematical theory provided in [START_REF] Lathuilière | Méthode de Décomposition de Domaine Pour Les Équations Du Transport Simplifié En Neutronique[END_REF]. Also, COCAGNE does not support globally non-conforming triangulations.

In order to simplify the implementation, the triangulation is created such that the macroscopic cross sections are constant on each cell, each color correspond to a material cell in Figure 2 middle. Indeed, the quadrature to evaluate the matrix coefficients do not consider the macroscopic cross sections. The geometry of the reactor is described by a grid of material (not necessarily conforming). Over this grid, we add another grid, called a partition. Each cell of the partition is called a subdomain. Then the triangulation is constructed from the sum of the material grid and the partition. We consider two types of triangulation, the globally conforming and the non-conforming. In the case of non-conforming, the subdomains have a conforming mesh or triangulation. One advantage of the non-conforming mesh is that each subdomain has its own mesh structure independent to the ones of its neighbours. Each subdomain can be meshed separately.

When one wants to have a better resolution of the problem in some interest points, the triangulation needs to be refined near these points. However, a globally Cartesian conforming mesh refinement costs a lot, because one has to create new elements along all the directions and it increases a lot the number of degrees of freedom. A less costed way is to use a non-conforming mesh and to refine it locally, in this case the new degrees of freedom are only create in the area of interest. First we present the MINOS solver, we show some results obtained with this solver for a checkerboard test case and for a pressurized water reactor (PWR). Those results are published in [START_REF] Giret | Criticality Computation with Finite Element Method on Non-Conforming Meshes[END_REF]. All the illustrations in this chapter are done with a Cartesian mesh, but all our study works fine with any type of triangulations.

MINOS Solver

The MINOS solver computes an approximated solution of the SP N multigroup equations with a mixed finite element method. More precisely, it uses Raviart-Thomas-Nédélec finite element method on Cartesian or hexagonal meshes [START_REF] Schneider | Éléments Finis Mixte Duaux Pour La Résolution Numérique de l'Équation de La Diffusion Neutronique En Géométrie Hexagonale[END_REF]. We present here our algorithm for the L 2 -jump domain decomposition method. An optimized Schwarz domain decomposition method is already available in MINOS [START_REF] Jamelot | Domain Decomposition for the Neutron SP n Equations[END_REF][START_REF] Jamelot | Fast Non-Overlapping Schwarz Domain Decomposition Methods for Solving the Neutron Diffusion Equation[END_REF].

Let us consider the multigroup SP N transport equation with G energy groups, N an odd integer. We denote R the domain of the reactor and D = 1, 2, 3 its dimension, R ⊂ R D .

Let Ñ be a non-null integer and ( R I ) Ñ I=1 be a partition of R into Ñ subdomains:

     R I ∩ R J = ∅, ∀I, J ∈ {1, • • • , Ñ }, R = Ñ I=1 R I .
We denote V I the set of all the indexes of the neighbour subdomains of the subdomain I:

V I = J | R I ∩ R J = ∅ and its Hausdorff dimension is D -1 .
For a subdomain R I , I = 1, Ñ , its interface with the subdomain R J , J ∈ V I , is denoted by Γ I,J . By construction, we have Γ I,J = Γ J,I . In Figure 7.1, we represent a Partition of a cuboid with 8 subcuboids with in 3 dimensions (D = 3). We want to determine the set of the neighbours of the subdomain R 8 . The intersection between R 8 and R 7 is a surface, second figure in Figure 7.1, thus its Hausdorff dimension is 2. But the intersection of R 8 and R 5 (resp. R 8 and R 1 ) is an edge (resp. a vertex), third (resp. fourth) figure in Figure 7.1, its Hausdorff dimension is 1 (resp. 0). Therefore, in this case V 8 is {4, 6, 7}.

R 1 R 2 R 3 R 4 R 5 R 6 R 7 R 8 Figure 7.1: Intersection of two different subdomains.
We recall that for a SP N method there are N = N +1 2 harmonics for the current p and for the flux φ. We denote, for I = 1, Ñ , by n p,I (resp. n φ,I ) the number of finite element degrees of freedom of the current p (resp. φ) in the subdomain R I . Moreover, we denote, for I = 1, ÑS , by n S,I the number of finite element degrees of freedom of the Lagrangre multiplier φ S on the interface I. We denote also P (resp. Φ and Λ) the vector of the degrees of freedom sorted by energy groups:

P = (P g ) g=1,G for p; Φ = (Φ g ) g=1,G for φ; Λ = (Λ g ) g=1,G for φ S .
For every energy groups g, the vector current P g is sorted by directions (d = 1, D), subdomains (I = 1, Ñ ), then every finite element degrees of freedom (i = 1, n p,I ) contains N harmonics:

P g = (P g d ) d=1,D ; P g d = (P g d,I ) I=1, Ñ ; P g d,I = (P g d,I,i ) i=1,n p,I ; where P g d,I,i ∈ R N .

Likewise, for every energy groups g, the scalar flux Φ g is sorted by subdomains (I = 1, Ñ ), then every finite element degrees of freedom (i = 1, n φ,I ) contains N harmonics:

Φ g = (Φ g I ) I=1, Ñ ; Φ g d,I = (Φ g d,I,i ) i=1,n φ,I ; where Φ g d,I,i ∈ R N .
We denote by ÑS the number of interfaces in the partition (R I ) I=1, Ñ . For every energy groups g, the scalar Lagrange multiplier Λ g is sorted by interfaces (I = 1, ÑS ), then every finite element degrees of freedom (i = 1, n S,I ) contains N harmonics:

Λ g = (Λ g I ) I=1, ÑS ; Λ g d,I = (Λ g d,I,i ) i=1,n S,I ; where Λ g d,I,i ∈ R N .
As the current p and the flux φ are defined over the domain R, P and Φ are indexed by subdomains, whereas the Lagrange multiplier φ S is defined on the interfaces of the partition ( R I ) I=1, Ñ , Λ is indexed by the interfaces. The harmonics and the finite element degrees of freedom are always considered together in the following.

From the variational formulation of the L 2 -jump domain decomposition method established in § 6, one obtains the following linear system:

   -A B -C B T T 0 -C T 0 0      P Φ Λ   = 1 k eff   0 0 0 0 M f 0 0 0 0     P Φ Λ   . (7.1)
The matrices obtained are also defined by blocks. Their blocks are defined as follow:

A = (A g,g ) g,g =1,G ; A g,g = (A g,g d,d ) d,d =1,D ; A g,g d,d = diag (A g,g d,d ,I ) I=1, Ñ ; B = diag [(B) g=1,G ] ; B = [(B d ) d=1,D ] T ; B d = diag [(B d,I ) I=1 Ñ ] ; T = (T g,g ) g,g =1,G ; T g,g = diag (T g,g I ) I=1, Ñ ; C = diag [(C) g=1,G ] C = (C d ) D d=1 T ; C d = (1 V I (J)C d,I,J ) I,J=1 Ñ .
Where 1 V I is the characteristic function of the set V I .

One can remark that there is no directional coupling in B and C because there are vectorial.

Remark 7.1. In MINOS, when d = 3, meshes are extrusion of 2D Cartesian meshes, thus it stands that for all d = 3, A g,g d,3 = A g,g 3,d = 0.

Remark 7.2. Moreover, in the case of a Cartesian mesh, one can remark that the matrices A g,g are block diagonals, for all d and d , d = d , A g,g d,d = 0, whereas in hexagonal meshes their exist a directional coupling in the matrices A g,g .

The Power Inverse Iteration

We look for the criticality which is the inverse of the greatest eigenvalue of the inverse transport operator that is to say the smallest eigenvalue of the transport operator. To evaluate this eigenvalue, we use the inverse power iteration which reads:

Algorithm 1 Power Inverse Iteration 1: initial state (P 0 , Φ 0 , Λ 0 , k eff,0 ) 2: S 0 ← M f Φ 0 3: n ← 1 4: until convergence do 5:
Solve :

6:    -A B -C B T T 0 -C T 0 0      P n Φ n Λ n   = 1 k eff,n-1   0 S n-1 0   7: S n ← M f Φ n 8: k eff,n ← k eff,n-1 S n |S n S n-1 |S n 9: n ← n + 1 10: end until
We denote the residual of the method at the n th iteration by:

n = k -1 eff,n S n -k -1 eff,n-1 S n-1 ∞ k -1 eff,n S n 1
The norm • ∞ (resp. • 1 ) is given, for any S in R N , N in N, by S ∞ = max m=1,N |S m | (resp. S 1 = N m=1 |S m |). Then, the convergence criteria is obtained when n is less than a given precision. The convergence rate of this algorithm is governed by the ratio between the first and second eigenvalue. In our resolution, it is the one that leads the convergence of all the algorithm, thus to accelerate it we use the Chebyshev acceleration [START_REF] Varga | Matrix Iterative Analysis[END_REF]. This acceleration consists in approximating the solution with a linear combination of the previous iteration fluxes Φ. The coefficients of this combination are computed from evaluating a Chebyshev polynomial on the approximated eigenvalue inverse. The inverse power iterations are called outer iterations.

Gauss-Seidel on the Energy Blocks

At each outer iteration, in order to solve the linear system in Algorithm 1 line 6 and thanks to the energy block structure of the system, we use the Gauss-Seidel iteration. This one is given in Algorithm 2.

Algorithm 2 Gaus-Seidel Iteration

1:   P n,0 Φ n,0 Λ n,0   ←   P n-1 Φ n-1 Λ n-1   2: n gs ← 0 3: until convergence do 4:
n gs ← n gs + 1 5:

for g = 1 to G do 6:

S g p,n,ngs ← g <g A g,g P g n,ngs + g >g A g,g P g n,ngs-1 7:

S g φ,n,ngs ← 1 k eff,n-1 S g n-1 - g <g T g,g Φ g n,ngs - g >g
T g,g Φ g n,ngs-1 8:

Solve :

9:   -A g,g B -C B T T g,g 0 -C T 0 0     P g n,ngs
  P n Φ n Λ n   ←   P n,ngs Φ n,ngs Λ n,ngs  
In practice, the inverse power iteration leads the convergence, so that a single iteration for the Gauss-Seidel iteration is actually required. In the case where there is no up-scattering,

A and T are lower triangular and one iteration of the Gauss-Seidel iteration corresponds to an exact resolution of the problem.

Remark 7.3. The hypothesis of weak up-scattering is generally valid for pressurized water reactor (PWR).

Flux Substitution and Alternative Direction Iteration (ADI)

We are interested now in the resolution of the linear system in Algorithm 2 line 9. We consider now that the matrices T g,g are easy to inverse. For instance, with well suited finite element degrees of freedom, T g,g can be diagonal.

Remark 7.4. In the case of RTN finite element on rectangles (d = 2) or cuboids (d = 3), T g,g are diagonal.

With this remark in mind, we rewrite our system in two problems as follow:

W g,g C C T 0 P g n,ngs Λ g n,ngs = S g n,ngs 0 ; (7.2)
and Φ g n,ngs = (T g,g ) -1 (S g φ,n,ngs -B T P g n,ngs );

(7.3)
where W g,g = A g,g + B(T g,g ) -1 B T and S g n,ngs = B(T g,g ) -1 S g φ,n,ngs -S g p,n,ngs . The matrices W g,g have the following block representation:

W g,g = (W g,g d,d ) d,d =1,D ; W g,g d,d = A g,g d,d + B T d (T g,g ) -1 B d .
Remark 7.5. In the case of Cartesian meshes, when d = d, it stands W g,g d,d = -B T d (T g,g ) -1 B d . In order to solve Equation (7.2), we use the direction with the alternative direction iteration (ADI) which is a Gauss-Seidel iteration over the directional blocks. This algorithm is given in Algorithm 3. The use of the ADI avoids us to invert W g,g globally but only its diagonal blocks W g,g d,d

(see Algotihme 3 line 10). The matrices W g,g d,d are symmetric positive definite. The iterations of the ADI are called inner iterations. As for the energy Gauss-Seidel iteration, one single iteration is enough for almost every applications. This is mainly the case in MINOS. One can remark that the ADI is exact when there is only one direction (D = 1).

Current Substitution

The linear system in Algorithm 3 line 10 can be simplified by substituting the current expression obtained with the first line of this system in the second one. This substitution leads to the following system: U npgc ← P g d R npgc+1

C T d (W g,g d,d ) -1 C d Λ g n,ngs,nadi = C T d (W g,g d,d ) -1 S g n,
13:

δ npgc+1 ← U npgc |R npgc+1
14:

β npgc ← δ npgc+1
δn pgc 15:

D npgc+1 ← D npgc -β npgc D npgc 16:
n pgc ← n pgc + 1 17: end while In Algorithm 4, the matrix P g d is the preconditioner. In our implementation, the preconditioner is taken as the inverse of the block diagonal matrix of C T d (W g d,d ) -1 C d , thus:

P g d = diag (C T d,I,J (W g,g d,d,I ) -1 C d,I,J + C T d,J,I (W g,g d,d,J ) -1 C d,J,I ) I=1, Ñ J∈V I -1 .

Special Case for the Cartesian Meshes

In the case where the partition is a Cartesian mesh, the matrices C T (W g,g ) -1 C are block diagonal due to that the normal of the interfaces between subdomains are collinear to the RTN basis functions of the current on the coincident borders of the subdomain. Indeed, the block corresponds to the line of subdomains in one direction. We show it on a 2dimension geometry with 3 × 2 domain decomposition, see Figure 7.2. Due to the form of the matrix, there is no coupling between some sets of Lagrange multiplier. In the example describe in Figure 7.2, those sets are: {Λ 1,2 , Λ 2,3 }, {Λ 4,5 , Λ 5,6 }, {Λ 1,4 }, {Λ 2,5 }, {Λ 3,6 }. We implement the gradient conjugate method in MINOS with these remarks in mind as it is done in [START_REF] Lathuilière | Méthode de Décomposition de Domaine Pour Les Équations Du Transport Simplifié En Neutronique[END_REF].

The final Algorithm

Here we regroup all the algorithm and we obtain the global algorithm used in MINOS. We denote by the precision required.

Algorithm 5 Global Algorithm

1: # Start Power Inverse Iteration 2: initial state (P 0 , Φ 0 , Λ 0 , k eff,0 ) 3: S 0 ← M f Φ 0 4: n ← 1 5: until convergence do 6:

# Start Gauss-Seidel Solve with PCG:

7:   P n,0 Φ n,0 Λ n,0   ←   P n-1 Φ n-1 Λ n-1   R 1 R 2 R 3 R 4 R 5 R 6 Γ 1,2 Γ 4,5 Γ 2,3 Γ 5,
24: # End Gauss-Seidel 37:

C T d (W g,g d,d ) -1 C d Λ g n,ngs,nadi = C T d (W
S n ← M f Φ n 38: k eff,n ← k eff,n-1 S n |S n S n-1 |S n 39:
n ← n + 1 40: end until 41: # End Power Inverse Iteration

Further comments on the Algorithm

The algorithm is based on nested iterative method. For some of these iterative methods, the number of iterations is blocked to one and thus the corresponding linear system is solved inexactly. In order to improve the convergence of the algorithm, the user can fix the maximun number of outer iterations, the maximum number of inner iterations, and the stopping criterium of the power inverse iteration. In the idea of [START_REF] Ern | Adaptive Inexact Newton Methods with A Posteriori Stopping Criteria for Nonlinear Diffusion PDEs[END_REF], one could use an a posteriori error estimate in order to optimize the number of iterations automatically. The parallelization of the L 2 -jump method has not been done in this work as we focus our work on the numerical analysis. According to [START_REF] Lathuilière | Méthode de Décomposition de Domaine Pour Les Équations Du Transport Simplifié En Neutronique[END_REF], the parallelization of this method brings some problem of load balancing. Indeed, in order to have a good load balancing, one have to consider each direction of the current independently. The non-overlaping Schwarz domain decomposition method (OSM) has been implemented in MINOS. This method consists in coupling each subdomain with some Robin interface conditions. Moreover, the coefficient appearing in these conditions can be optimized in order to improve the rate of convergence of the iterative DDM solver [START_REF] Jamelot | Fast Non-Overlapping Schwarz Domain Decomposition Methods for Solving the Neutron Diffusion Equation[END_REF]. The optimization is done by studying some asymptotic problems see [START_REF] Nataf | Convergence rate of some domain decomposition methods for overlapping and nonoverlapping subdomains[END_REF] for some order 1 approximation and [START_REF] Japhet | The optimized order 2 method: Application to convection-diffusion problems[END_REF] for the second order approximation. For a general overview on OSM, one can refer to [START_REF] Gander | Optimized Schwarz Methods[END_REF]. The spatial linear system to solve is not symmetric and it is solved by a Jacobi method. This resolution could be improved by using some more efficient Krylov methods like GMRES method [START_REF] Saad | A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems[END_REF]. Also, unlike the L 2 -jump method, OSM cannot treat non-conforming method. The linear system obtained by the OSM is known to have a better condition number than the one obtained by the L 2 -jump method. In [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions: The Finite Element Case[END_REF][START_REF] Japhet | A New Cement to Glue Non-Conforming Grids with Robin Interface Conditions: The p 1 Finite Element[END_REF], the authors proposed to combine these two Domain Decomposition Methods (DDM): the L 2 -jump method and the OSM. Their method consists in using a Lagrange multiplier on the interface which satisfies some Robin boundary conditions. The authors use a uniform coefficient in the Robin boundary condition which correspond to the optimized uniform coefficient for the OSM (see [START_REF] Gander | Optimized Schwarz Methods[END_REF]). This method could be an alternative to the L 2 -jump method when this one converges slowly. On the other hand, the numerical analysis done in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions: The Finite Element Case[END_REF][START_REF] Japhet | A New Cement to Glue Non-Conforming Grids with Robin Interface Conditions: The p 1 Finite Element[END_REF] only treats the case of regular solution. In order to treat non-conforming triangulation with Finite Element Method, one can also use the mortar element method. This method was first introduced for coupling Finite Element Method with Spectral Method [START_REF] Bernardi | Couplage de methodes spectrale et d'elements finis: Premiers resultats d'approximation[END_REF] and then it was adapted to couple two Finite Element Methods [START_REF] Bernardi | Domain decomposition by the mortar element method[END_REF]Wohl01]. This method consists in imposing the interface conditions directly in the discrete spaces and not with a Lagrange multiplier on the interfaces. This method cannot be parallelized, but it is more flexible than the L 2 -jump as it allows more general local refinements.

The solver MINOS also solves the kinetic neutron SP N equations, thus the L 2 -jump DDM could be used in this case. The parallelization in time of the solver was studied by Mula [START_REF] Baudron | The parareal in time algorithm applied to the kinetic neutron diffusion equation[END_REF]. This parallelization in time is based on the parareal method, which is an iterative technique where, at each iteration, a predictor corrector propagation is proposed based on two propagators :

-a coarse propagator (the predictor); -a fine propagator (the corrector).

Moreover, one can reduce the memory cost of the parareal method by using a reduced basis framework [START_REF] Maday | A generalized empirical interpolation method: Application of reduced basis techniques to data assimilation[END_REF]. A further work could be to couple this time parallelization with the L 2 -jump DDM. This time parallelization is also implemented in the solver MINARET in APOLLO3 R platform [START_REF] Baudron | MINARET: Towards a time-dependent neutron transport parallel solver[END_REF], which solves the multigroup transport equations with discrete ordinate (S N ) for the angular dependence (Ω) and discontinuous Galerkin method for the spatial dependence (x). This solver has a Schwarz DDM [START_REF] Odry | Méthode de Décomposition de Domaine Avec Parallélisme Hybride et Accélération Non Linéaire Pour La Résolution de l'équation Du Transport Sn En[END_REF]. According to [START_REF] Antonietti | A Class of Domain Decomposition Preconditioners for hp-Discontinuous Galerkin Finite Element Methods[END_REF][START_REF] Antonietti | Schwarz domain decomposition preconditioners for plane wave discontinuous Galerkin methods[END_REF], this could be improved thanks to a good preconditioner. Another way to parallelize in time, for the kinetic case, is the waveform relaxation proposed in [START_REF] Japhet | Optimized Schwarz waveform relaxation for porous media applications[END_REF]. The authors proposed a time DDM based on an optimized Schwarz method with relaxation. This method can handle different time steps between the subdomains.

In MINOS, the Power Inverse Iteration is accelerated with the Chebyshev acceleration which consists in taking the next approximation as a combination of the new approximation and the previous one. In MINARET, the Power Inverse Iteration is accelerated with the Diffusion Synthetic Acceleration (DSA) [START_REF] Moller | MINARET, a deterministic neutron transport solver for nuclear core calculations[END_REF]. This acceleration consists in solving an coarser angular problem, where for each discrete ordinate component of the flux we solve a diffusion problem. The DSA acceleration can reduce a lot the number of outer iterations. Nevertheless, this method is not parallel, thus is a bottleneck when using parallel capabilities of MINARET, according to [START_REF] Moller | Éléments Finis Courbes et Accélération Pour Le Transport de Neutrons[END_REF]. We could use the work in [START_REF] Antonietti | Domain decomposition preconditioning for discontinuous Galerkin approximations of convection-diffusion problems[END_REF] on DDM for Discontinuous Galerkin to parallelize the diffusion solver of the DSA.

Checkerboard testcase

The checkerboard testcase is an adaptation of the Maxwell eigenvalue benchmark proposed in [START_REF] Dauge | Benchmark Computations for Maxwell Equations for the Approximation of Highly Singular Solutions[END_REF] for the neutron diffusion eigenvalue problem with zero flux (Dirichlet) 7.2 : Checkerboard testcase boundary conditions, see [START_REF] Ciarlet | Numerical Analysis of the Mixed Finite Element Method for the Neutron Diffusion Eigenproblem with Heterogeneous Coefficients[END_REF] for Neumann boundary conditions. Set R =]0; 100[ 2 and divide it into four cells (R i ) i=1,4 as in Figure 7.3. We denote by D the diffusion coefficient. The coefficient D is piecewise constant such that:

D = D in R 1 ∪ R 4 ; D = 1 in R 2 ∪ R 3 ;
where D is a positive real. The nuclear data are set such that we solve the following problem:

R 1 R 3 R 2 R 4 D 1 = D D 2 = 1 D 3 = 1 D 4 = D
         1 D p + grad φ = 0 in R div p + φ = 1 k eff φ in R; φ = 0 on ∂R.
(7.5)

We compare the resolution to a fine resolution solution. The fine solution is computed on a 1000 × 1000 mesh with RT 4 elements and several inner iterations. We study the error with two different elements RT 0 and RT 1 on a conforming mesh. Those results are in table 7.1 for RT 0 element and in table 7.2 for RT 1 element, where the data are:

h: the meshsize, -N out : The number of outer iterations,

-ε k eff = |k h,eff -k eff | k eff
: the relative error on the criticity,

-ε φ = φ h -φ 0,R φ 0,R
: the relative error on the flux.

In the last line, we report the average rate of convergence of the computations. We use an uniform mesh for all this computations. Moreover, we set D = 5. Thanks to the numerical analysis done in Chapter 5, we know that the converging rate for the RT 0 elements is 2 min(ω ν , 1), where ω ν is the fundamental mode regularity. Thus one can conclude that the eigenfunction is regular. In the case of the RT 1 elements, the expected convergence rate is 4, the difference from the real one comes from that we compare to a fine mesh solution which is well described by RT 1 elements quickly. Now, we study the error using the L 2 -jump DDM with RT 0 elements on a non-conforming mesh,

Large Heavy Steel Reflector Reactor Core

We apply our method on a PWR-like reactor core with an heavy steel reflector. We use a benchmark similar as the one described in [START_REF] Sargeni | Coupling Effects in Large Reactor Cores: The Impact of Heavy and Conventional Reflectors on Power Distribution Perturbations[END_REF]. The reactor core is projected on a Cartesian mesh composed of 361 elements such that an assembly corresponds to one cell.

In the core there is 241 fuel assembly cells (in yellow in Figure 7.5) surrounded by 120 steel reflector cells (in blue in Figure 7.5). Each assembly is a set of 17 by 17 fuel rods, Figure 7.6. The coolant, here water, can pass through the core between the rods. As the physical nuclear data are the microscopique cross sections, first we have to evaluating the macroscopic cross sections. This evaluation is done as the same time as the homogenization of the core [START_REF] Sanchez | Assembly Homogenization Techniques for Core Calculations[END_REF][START_REF] Coste-Delclaux | Modélisation Du Phénomène d'autoprotection Dans Le Code de Transport Multigroupe APOLLO2[END_REF]. Each cell is homogenized on one of the three subgrids given in Figure 7.7.

After the numerical homogenization, we obtain piecewise constant macroscopic cross section. The macroscopic cross section are constant on each subcell. The finest homogenization meshes are placed where the assembly admits control rods and at the interface between the fuel and the reflector because it is the area of the core where we expect to have high variations of the flux. At the end of this process, the geometry is described by a mesh with 51 241 subcells and 229 different media. The mesh of the geometry is given in Figure 7.8. The modelization is done with 2 groups of energy and with SP 1 and SP 3 method. We compare the resolution with a conforming mesh against the resolution on a non-conforming mesh. The first one is created from the homogenized subgrid, instead the second one is the homogenized subgrid. The conforming mesh has 115 600 meshes and the non-conforming one has 51 241 meshes. In Table 7.5 and Table 7.4, we show the results for different uniform refinement on every meshes of the geometry for conforming meshes and non-conforming meshes. Chapter 8

Adaptive methods

A Posteriori Error Estimate

For low-regular solution, the increase of the finite element method order does not improve the approximation, the triangulation needs to be refined. In order to refine the triangulation without unnecessary degrees of freedom, one can use an a posteriori error estimate [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF] in order to refine only where the error is high. This method has been declined for different approximation: see [START_REF] Becker | A feed-back approach to error control in finite element methods: Basic analysis and examples[END_REF] for primal finite element, see [START_REF] Vohralik | A Posteriori Error Estimates for Lowest-Order Mixed Finite Element Discretization of Convection-Diffusion-Reaction Equations[END_REF][START_REF] Vohralik | Guaranteed and Fully Robust A Posteriori Error Estimates for Conforming Discretizations of Diffusion Problems with Discontinuous Coefficients[END_REF] for mixed finite element, see [START_REF] Antonietti | An hr-adaptive discontinuous Galerkin method for advection-diffusion problems[END_REF] for discontinuous Galerkin method, see [START_REF] Omnes | A posteriori error estimation for the discrete duality finite volume discretization of the Laplace equation[END_REF] for duality finite volume.

Another possible use of this technique occurs in the case of control rod cluster ejection. Indeed, in this situation the flux distribution evolves rapidly. In order to well approximate these flux modifications, one could use an adaptive mesh refinement method based on an a posteriori error estimate.

Derivation of an A Posteriori Error Estimate

We derive here an a posteriori error estimate for the one-group diffusion model for the mixed resolution with the Raviart-Thomas finite element method (see chapter 5). This work is an adaptation of [START_REF] Vohralik | A Posteriori Error Estimates for Efficiency and Error Control in Numerical Simulations[END_REF]. We recall that the variational source problem 2.16 associated to the one-group diffusion model reads:

Problem 8.1. For a given S f in L, find ζ in X, such that for all ξ in X, it stands:

c(ζ, ξ) = f (ξ).
As in § 2.2 and § 5.1, hypothesis 2.4 is supposed to be satisfied.

To approximate the solution of Problem 8.1, we use the Raviart-Thomas finite element method, which is described in § 5. We recall that this approximation reads:

Problem 8.2. For S f in L, find ζ h in X h such that for all ξ h in X h , it stands:

c(ζ h , ξ h ) = f (ξ h ).
As the discrete solution ζ h = (p h , φ h ) is in X the discrete flux φ h is not in V . Therefore, we denote ζh = (p h , φh ), where φh is a reconstruction of φ h in V . The obtention of φh is explained in § 8.1.2 below. We introduce the following bilinear forms:

B S : X × X → R (ζ, ξ) → -a(p, q) + t(φ, ψ) , and
B A : X × X → R (ζ, ξ) → b(p, ψ) -b(q, φ) ,
where a, b and t are defined on page 40. Thanks to hypothesis 2.4, the bilinear form B S is symmetric and coercive. Indeed, let ζ be in X, thus it stands:

B S (ζ, ζ) ≥ min((D * ) -2 , (Σ r,0 ) 2 * ) ζ 2 X .
Thus, B S helps to define a norm • S on X:

ζ 2 S = B S (ζ, ζ) + Σ -1 2 r,0 div p 2 0,R .
By construction the norm • S is equivalent to the natural norm • X . We recall that T h denote the triangulation, then it stands, for all ζ in X:

ζ 2 S = K∈T h ζ 2 S,K , where ζ S,K = ( D -1 2 p 2 0,K + Σ 1 2
r,0 φ 2 0,K + div p 2 0,K )

B(ζ, ξ) = c(ζ, (-q, ψ)) = f (ξ). (8.1)
We define also a norm on X, for all ζ in X:

|ζ| + = sup ξ∈X ξ S ≤1 B(ζ, ξ).
This norm is not equivalent to the full norm • X , due to the absence of the divergence term.

Lemma 8.3. Let ζ (resp. ζ h ) be the solution of problem 8.1 (resp. problem 8.2), and let ζh be a reconstruction of ζ h in Q × V . Thus, for all ξ ∈ X, it stands:

B(ζ -ζh , ξ) = R (S f -div p h -Σ r,0 φ h )ψ + R Σ r,0 (φ h -φh )ψ - R (D -1 p h + grad x φh ) • q. (8.2)
Proof. Let ξ be in X, thanks to equation 8.1, it stands:

B(ζ -ζh , ξ) = R S f ψ - R D -1 p h • q - R Σ r,0 φh ψ - R div p h ψ + R div q φh .
We recall that φh is in V , thus we can integrate by part the last integral:

B(ζ -ζh , ξ) = R (S f -div p h -Σ r,0 φh )ψ - R D -1 p h • q - R q • grad x φh .
To conclude the proof, one can add

R Σ r,0 φ h ψ - R Σ r,0 φ h ψ. Theorem
R,K = Σ -1 2 r,0 (S f -div p h -Σ r,0 φ h ) 0,K , the flux estimator by η F,K = D 1 2 (D -1 p h + grad x φh ) 0,K ,
and the non-conformity estimator by

η N C,K = ζ h -ζh S,K = Σ 1 2 r,0 (φ h -φh ) 0,K .
Then, it stands:

|ζ -ζ h | + ≤ 3 K∈T h η 2 N C,K 1 2 + K∈T h η 2 R,K + η 2 F,K 1 2 . (8.3)
Proof. Thanks to the triangle inequality, it stands:

|ζ -ζ h | + ≤ |ζ -ζh | + + | ζh -ζ h | + . (8.4)
First, we bound |ζ -ζh | + . From equation 8.2, it holds for all ξ in X:

B(ζ -ζh , ξ) = K∈T h K (S f -div p h -Σ r,0 φ h )ψ + K Σ r,0 (φ h -φh )ψ - K (D -1 p h + grad x φh ) • q .
Using Cauchy-Schwarz inequality on all the integrals, and the definition of the residual, the flux and the non-conformity estimators, it holds:

B(ζ -ζh , ξ) ≤ K∈T h η R,K Σ 1 2 r,0 ψ 0,K + η F,K D -1 2 q 0,K + K∈T h η N C,K Σ 1 2 r,0 ψ 0,K .
From the definition of the norm • S , one obtains:

B(ζ -ζh , ξ) ≤ K∈T h η 2 R,K + η 2 F,K 1 2 ξ S + K∈T h η 2 N C,K 1 2 ξ S .
From the definition of the norm | • | + and the previous inequality, we find that:

|ζ -ζh | + ≤ K∈T h η 2 R,K + η 2 F,K 1 2 + K∈T h η 2 N C,K 1 2
. Now, we bound the second term in the left hand side of inequality (8.4). Let ξ be in X, we look for an upper bound to:

B( ζh -ζ h , ξ) = B S ( ζh -ζ h , ξ) + B A ( ζh -ζ h , ξ). (8.5)
Using Cauchy-Schwarz in equation (8.5), one obtains:

B( ζh -ζ h , ξ) ≤ ζh -ζ h S ξ S + B A ( ζh -ζ h , ξ).
Moreover, for all ξ in X, it stands:

B A ( ζh -ζ h , ξ) = - R div q( φh -φ h )
Therefore, it holds:

B A ( ζh -ζ h , ξ) ≤ K∈T h η 2 N C,K 1 2 ξ S .
The a posteriori error estimate in theorem 8.4 is reliable (upper bound), however we have not proved its efficiency (lower bound).

Reconstruction of the Discrete Flux

To reconstruct the discrete flux φ h , we use the averaging operator as proposed in [START_REF] Burman | Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations[END_REF].

Other interpolation of non-smooth functions have been studied in [START_REF] Scott | Finite element interpolation of nonsmooth functions satisfying boundary conditions[END_REF]. Alternatively, one may use finite volume reconstruction [START_REF] Omnes | On the second-order convergence of a function reconstructed from finite volume approximations of the Laplace equation on Delaunay-Voronoi meshes[END_REF]. We recall that L k h is the approximation space of the flux for RTN [k] finite elements, k ≥ 0:

L k h = {ψ h ∈ L|∀K ∈ T h , ψ h|K ∈ Q k },
where Q k is the set of all polynomials of degree k in each direction. And on the other hand V k+1 h is:

V k+1 h = {ψ h ∈ V |∀K ∈ T h , ψ h|K ∈ Q k+1 },
which is conforming in V by definition. We consider a Lagrange basis of V k+1 h associated to the interpolation points (a i )

dimV k+1 h i=1 . We denote (ϕ i ) dimV k+1 h i=1 the basis of V k+1 h such that for all 1 ≤ i, j ≤ dimV k+1 h : ϕ i (a j ) = δ i,j .
For any x in R, we denote T x the set of all the elements of T h where x stands in:

T x = {K ∈ T h |x ∈ K}.
For any 1 ≤ i ≤ dimV k+1 h , we denote |T a i | the number of elements in T a i . For k ≥ 0, we define the average operator I k+1 from L k+1 h into V k+1 h defined by, for all ψ h in L k+1 h :

∀1 ≤ i ≤ dimV k+1 h , I k+1 (ψ h )(a i ) = 1 |T a i | K∈Ta i ψ h|K (a i ).
It is proven in [START_REF] Burman | Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations[END_REF] that, for all h, for all K ∈ T h and for all φ h ∈ L k h :

φ h -I k+1 (φ h ) 0,K √ h K k + 1 F ∈∂K [φ h ] 0,F . (8.6)
For fixed K ∈ T h and for all F = ∂K ∩ ∂K = ∅ (resp. F = ∂K ∩ ∂R), whose Hausdorff dimension is equal to d -1, where K ∈ T h , the jump of φ h across F is defined:

[φ h ] = φ h|K -φ h|K (resp. [φ h ] = -φ h|K ).
Remark 8.5. None of the operators I k+1 , for k ≥ 0, is defined on L 0 h . However one simply notices that L 0 h is a subset of k≥0 L k+1 h , so in practice one can choose which operator to use on L 0 h . Remark 8.6. From inequality (8.6), one can remark that the quality of the reconstruction of the discrete flux improves with k.

In order to evaluate the a posteriori error estimate described in theorem 8.4, for the RTN [k] finite element, k ≥ 0, we take

φ h = I k+1 (φ h ) = dimV k+1 h i=1   1 |T a i |   K∈Ta i ψ h|K (a i )   ϕ i   .
Therefore, φ h is a reconstruction of φ h on the Lagrange finite element Q k+1 in V .

In figure 8.1 we represent the reconstruction in Q 1 of a discrete flux given by a RTN [0] finite elements in one dimension (d = 1).

Application to the Resolution

We use the a posteriori error estimate described in (8.3) to improve the resolution of the neutron diffusion problem 2.13 with RTN finite elements. Indeed, this estimator is used to refine locally the triangulation in order to homogenized the error all over the triangulation. This method is called adaptive mesh refinement (AMR). The error estimate in (8.3), is a global error estimate, it can be localized if one defines the local error estimate on K, K ∈ T h , by

η K = 9η 2 N C,K + η 2 R,K + η 2 F,K . (8.7)
This local estimator is used to determine whether an element of the triangulation should be refined or not. We explain below where this estimator is used in our general algorithm. In section 7.1.1, we presented the power inverse iteration used to solve the neutron diffusion problem 2.13. We now modify this algorithm such that at each outer iteration,we determine whether or not the triangulation should be refined: for that we compute max K∈T h η K and compare it to a threshold value AM R . If this is the case, the triangulation refinement is done by ordering the elements according to the size of the error indicator η K and to refine the elements that make up a certain percentage of the total error estimate. The power inverse iteration with AMR is described in algorithm 6.

Obviously, if the triangulation has been refined in the previous iteration, the finite element matrices must be recalculated before line 6 in algorithm 6.

An alternative strategy for the triangulation refinement is to use error-balancing strategy proposed in [START_REF] Bangerth | Adaptive Finite Element Methods for Differential Equations[END_REF]. This strategy consists in refining some elements and derefining some groups of elements such that the local errors are comparable on all elements:

∀K ∈ T h , η K ∼ AM R N elt ,
where N elt is the number of elements of the triangulation T h . Solve :

6:    -A B -C B T T 0 -C T 0 0      P n Φ n Λ n   = 1 k eff,n-1   0 S n-1 0   7:
for K ∈ T h do 8:

Eval: η K 9:

end for 10:

if max K∈T h > AM R then 11:
Triangulation refinement 12:

Reconstruction of the finite element matrices 13:

end if 14:

S n ← M f Φ n 15: k eff,n ← k eff,n-1 S n |S n S n-1 |S n 16: n ← n + 1 17: end until

Application in One Dimension

We recall from § 1, that this problem reads: Find (λ, p, φ) in R × Q × V \{0}, such that:

D -1 p + ∂ x φ = 0 ∂ x p + Σ r,0 φ = λνΣ f φ .
The domain R is ]0, 1[ and the coefficient D is given by:

∀x ∈ R, D(x) =          5 x < 1 2 1 x ≥ 1 2
.

The reference solution, represented in figure 8.2 is computed on a fine triangulation with N elt = 1000 elements, and the eigenvalue is 0.040477. As the approximate flux φ h is computed in L 0 h , its gradient is not well approximated. Thus, the triangulation is expected to be refined in two regions: where the gradient of the flux is high, and where there is a jump of D, at x = 1 2 . In table 8.1, we compare the resolution with and without AMR. In the case with the triangulation refinement (I-AMR), we start with a 2-element triangulation. At each outer iteration, if the triangulation must be refined (if max K∈T h > AM R ), elements totalling 50% of the total error are refined. The final triangulation, in Figure 8.3, has N elt = 178 elements. The triangulation has been refined where we expected it to be. For the comparison, the case without AMR (II-No AMR), is done on a uniform triangulation with N elt = 178 elements. We remark that convergence without AMR is obtained after N it = 17 outer iterations, compared to N it = 21 outer iterations with AMR.

At last, in the column (III-AMR) of table 8.1, we present the error of the method with AMR stopped after N it = 17 outer iterations.

In figure 8.4, we plot the number of elements at each outer iteration. One can remark that the triangulation is refined in the first outer iterations only. In the next example, we change the value of coefficient D to:

∀x ∈ R, D(x) =                              1 x < 1 4 10 1 4 ≤ x < 1 2 2 1 2 ≤ x < 3 4 5 x ≥ 3 4
.

The reference solution is plotted in figure 8.5. The gradient of the flux is high where the coefficient D is small, thus we expect that the triangulation is refined there. In figure 8.6, we represent the final triangulation of the resolution with AMR: we see that the triangulation is refined where we expected. In case (I-AMR), the adaptive mesh refinement of the triangulation is refined from the 1 st iteration to the 14 th iteration. Thus, in case (III-AMR), the triangulations is not yet converged when the power inverse iteration is stopped. Therefore, the error in case (III-AMR) is much larger than in case (I-AMR). In these two examples, we note that first the triangulation is adapted to the problem, and then in a second step, the power inverse iteration converges on the final, well fitted, triangulation.

Application in Two Dimensions

In this subsection, we consider the source problem. We recall from problem 2.14, that this problem reads: For S f given, find (p, φ) in Q × V \{0}, such that:

D -1 p + grad x φ = 0 in R, div p + Σ r,0 φ = S f in R. (8.8)
The domain R is ]0, 1[×]0, 1[ and the coefficient D is given by: Above, the coefficient D is set to 100. In that case, the regularity exponent is roughly equal to 0.13 [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF]. We choose S f := 2x(x -1) + 2y(y -1) + x(x -1)y(y -1). The solution is plotted in Figure 8.9 on a uniformly refined grid. We use the a posteriori error estimate given in theorem 8.4 to perform an adaptive mesh refinement. We start with an initial mesh with 12 × 12 elements and we keep the logical structure during the AMR; namely we keep a grid structure like N x × N y elements grids. The final mesh obtained is given in figure 8.10 and has 43 × 35 elements. As predicted, the mesh is refined near the singularity in (0.5, 0.5). The loss of symmetry in the final mesh (figure 8.10) comes from the fact that only elements totalling 50% of the total error are refined at each iteration. The algorithm which selects the elements to be refined could be modified. Since the logical grid structure is chosen in this two dimension application, one could add all errors on a given line or column, and use these aggregated errors to perform the AMR on lines or columns directly. Finally the use of a DDM, with AMR by subdomain, should allow a better approximation at a lesser cost.

R 1 R 3 R 2 R 4 D 1 = D D 2 = 1 D 3 = 1 D 4 = D

Conclusion

We showed in § 8.1.3, that AMR improves the resolution of an eigenvalue problem with power inverse iteration (in one dimension). Increasing the resolution of the eigenfunction gives us have a better approximation of the eigenvalue. For the source problem (in two dimensions), the solution is refined where the largest variations occur. In the case of the multigroup equations, we recall that there are G neutron fluxes. In [START_REF] Ragusa | A simple Hessian-based 3D mesh adaptation technique with applications to the multigroup diffusion equations[END_REF][START_REF] Wang | Three-Dimensional h-Adaptivity for the Multigroup Neutron Diffusion Equations[END_REF], the authors propose an a posteriori error estimates for the multigroup diffusion equations under their primal setting, for which the G triangulations are refined independently from one another. Thus, one needs to project the neutron fluxes on the different triangulations. This can increase the memory cost or the computation time of the method. To avoid these projections, one can derive a fission source error estimate and compute the G neutron fluxes on the same (refined) triangulation.

Chapter 9 Conclusion

On the model. In this work, we study the approximation of the multigroup SP N equations, which are an approximation of the neutron transport equation, with mixed Finite Element Method. A continuation work could be to extend this study to the P N transport equation and to the transport equation. Moreover, we only study the spatial approximation error, without taking into account the other discretizations (multigroup and SP N ), a more complete study of the resolution of the neutron transport equation using the multigroup SP N equations would consider all the discretizations. In our model, the cross sections are homogenized, and thus they are piecewise regular. We say that the diffusion equations or the multigroup SP N equations are in a realistic configuration when the cross sections obtained after homogenization are piecewise regular. However, the cross sections oscillate rapidly, and the homogenization process removes these oscillations in the model. A way to keep these oscillations in the model would be to use a multiscale finite element method (MS-FEM, [START_REF] Engquist | The heterogeneous multiscale methods[END_REF]) or a heterogeneous multiscale method (HMM, [START_REF] Ciarlet | Finite Element Heterogeneous Multiscale Method for the Helmholtz Equation[END_REF]). These methods consist in solving the problem in two steps, a macro resolution and a micro resolution. The micro resolution is done on a fine triangulation of a small part of the geometry, then the solution is injected into the macro resolution in order to incorporate the fine scale effect of the cross sections into the resolution.

On the numerical analysis. We propose the numerical analysis of the mixed finite element method in order to approximate the solution of the multigroup SP N equations. In a realistic configuration, crosspoints (intersection points of more than 3 materials) are allowed and are common, the solution can have a low-regularity. The numerical analysis proposed in this work is done taking into account the low-regularity of the solution.

To our knowledge it is the first rigorous numerical analysis of the model with realistic configuration. Moreover, we extend this analysis to the associated eigenvalue problem. We prove the norm convergence of the discrete operator toward the continuous one which ensures that there is no spectral pollution. Therefore, in the limit, each discrete eigenvalue corresponds to only one continuous eigenvalue. To our knowledge it is the first rigorous numerical analysis of the eigenvalue problem with an absorption term. As a matter of fact, the theory only addresses case of no absorption.

Part III : Chapter 9 : Conclusion Knowing a priori error estimates of the method is a necessary first step in order to estimate the propagation of the uncertainty on the cross sections. Indeed, from the work of Charrier in [START_REF] Charrier | Strong and weak error estimates for the solutions of elliptic partial differential equations with random coefficients[END_REF], the error of the approximation is given by the deterministic discretization error and the stochastic error due to the coefficient uncertainties. Our study gives us the first term.

On the algorithmic side. We also propose a domain decomposition method (DDM), L 2jump, which can deal with globally non-conforming triangulations. This method consists in adding a Lagrange multiplier on each interface which ensures the conformity of the method and the projection between the local conforming triangulations. The numerical analysis of this DDM is done in this manuscript. Again, to our knowledge, it is the first time such analysis is done for a DDM with low-regularity solution. We implemented the L 2 -jump DDM in the solver MINOS of the APOLLO3 R platform. The solver MINOS includes also a non-overlapping optimized Schwarz DDM. This latter DDM cannot treat non-conforming triangulations, and we emphasize that the numerical analysis done for this method only considers the case of regular solution.

To decrease the number of outer iterations, we propose an implementation of an adaptive mesh refinement. In order to choose where to refine the triangulation, we use an a posteriori error estimate. Thus, the triangulation is better fitted to solve accurately the problem inside the power inverse iteration. Here we propose a proof of Lemma 3.4. In order to do that, we use the following proposition.

Proposition C.1. For any function f and g in L, for any α in ]0, 1], it stands:

2(f, g) 0,R ≤ α f 2 0,R + 1 α g 2 0,R .
Next we prove Lemma 3.4.

Proof. Let ψ be in L. The norm of Hψ in L reads: For all 1 ≤ g ≤ G, let (α g 1 , . . . , α g N -1 ) be a ( N -1)-uplet in ]0, 1] N -1 . Then, from Proposition C.1, it stands:

Hψ 2 L = G g=1   ψ g
Hψ 2 L ≥ G g=1   2 N h=2 ψ g h 2 0,R + ψ g 1 2 0,R - N -1 h=1 α g h ψ g h 2 0,R + 1 α g h ψ g h+1 2 0,R   ; ≥ G g=1   (1 -α g 1 ) ψ g 1 2 0,R + N -1 h=2 2 -α g h - 1 α g h-1 ψ g h 2 0,R + 2 - 1 α g N -1 ψ g N 2 0,R .
For all 1 ≤ g ≤ G, we denote:

• c g 1 = 1 -α g 1 . • ∀2 ≤ h ≤ N -1, c g h = 2 -α g h - 1 α g h-1 ; • c g N = 2 - 1 α g N -1
. Therefore, it stands:

Hψ 2 L ≥ α 2 H ψ 2 L ,
where α 2 H = min G g=1 min N h=1 (c h ). Now we have to choose the G family (α g h ) N -1 h=1 such that α H is strictly positive. A sufficient condition is that for all 1 ≤ g ≤ G, and for all 1 ≤ h ≤ N , c g h is strictly positive: This condition imply that, for all 1 ≤ g ≤ G:

c g 1 > 0 ⇒ 1 -α g 1 > 0 ⇒ α g 1 < 1 ∀2 ≤ h ≤ N -1, c g h > 0 ⇒ 2 -α g h - 1 α g h-1 > 0 ⇒ α g h-1 > 1 2 -α g h c g N -1 > 0 ⇒ 2 - 1 α g N -1 > 0 ⇒ α g N -1 > 1/2
For all 1 ≤ g ≤ G, one can remark that, for any

1 ≤ h ≤ N , if α g h-1 < 1 then 1 2-α g h < 1.
Thus, one can choose, for all 1 ≤ g ≤ G, α g N in 1 2 ; 1 and then α g h in 1 2α h+1

, 1 for h = 1, N -1 such that c g h are positive. From equation (C.1), thanks to the triangle inequality and Proposition C.1, one obtains that:

Hψ 2 L ≤ G g=1   ψ g N 2 0,R + 2 N -1 h=1 ψ g h 2 0,R + ψ g h+1 2 0,R   .
Therefore, it stands that:

Hψ L ≤ √ 2 ψ L .
C.2 Results for Chapters 5 and 6

Let (T h ) h be a given regular family of triangulations. We call K := [0, 1] d the reference element. Let h be given. For every K ∈ T h , we denote by x = F K (x) := A K x + b K ,

A K ∈ R d×d , b K ∈ R d
, the map from K to K. Introducing h K = diam(K) for all K ∈ T h , one may bound A K , (A K ) -1 , |det(A K )| with respect to h K . The change of variable formulas from K to K, and vice versa, can be found e.g. in [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF]§1].

Proof. (of Lemma 5.9) We follow [START_REF] Belgacem | Some Nonstandard Finite Element Estimates with Applications to 3D Poisson and Signorini Problems[END_REF]§2]. Given ψ h ∈ L k h , one has ψ h ∈ H µ (R), for all µ < 1/2. By the definition of the norm of H µ (R), we have the following equalities: h K . To estimate the remaining part, we recall that, for any K ∈ T h and any x ∈ K, it holds that, by going back the reference space, applying (cf. [Gris85, (1.3.2.12)]) on K and then going to the physical space: R\K 1 |x -y| d+2µ dy 1 ρ ∂K (x) 2µ , where ρ ∂K (x) = inf y∈∂K |x -y|. Thus we have:

ψ h 2 µ,R = ψ h 2 0,R + R R |ψ h (x) -ψ h (
K R\K |ψ h (x) -ψ h (y)| 2 |x -y| d+2µ dydx = K ∈T h K =K K K |ψ h (x) -ψ h (y)| 2 |x -y| d+2µ dydx K ∈T h K =K K K ψ h (x) 2 + ψ h (y) 2 |x -y| d+2µ dydx K ψ h (x) 2 ρ ∂K (x) 2µ dx. (C.4)
Going back to the reference element K and introducing ψ h|K (x) = ψ(x), it stands:

K ψ h (x) 2 ρ ∂K (x) 2µ dx h d-2µ K K ψ(x) 2 ρ ∂ K (x) 2µ dx.
Because µ < 1/2 (cf. [Gris85, Theorem 1.4.4.4]), ψ → ( K ψ(x) 2 ρ ∂ K (x) -2µ dx) 1/2 is a norm on Lk = Q k,k,k ( K). Thanks to the equivalence of the norms on finite dimensional vector spaces, one gets

K ψ h (x) 2 ρ ∂K (x) 2µ dx h d-2µ K ψ 2 0, K .
Finally, going back to element K, we know that ψ 2 0, K h -d K ψ h 2 0,K . Hence using (C.4) and the results that follow, we have:

K R\K |ψ h (x) -ψ h (y)| 2 |x -y| d+2µ dydx h -2µ K ψ h 2 0,K . (C.5)
Starting from (C.2) using (C.3) and (C.5), we obtain finally the global bound:

ψ h µ,R h -µ min ψ h 0,R .
As the family of triangulations is regular + , one has h -µ min h (θ-2)µ , which concludes the proof.

Proof. (of Lemma 6.2) For l = c, f , we introduce the operators from the normal trace spaces (H(div , R) ∩ H µ (R)) • n |Γ f c to the discrete spaces of normal traces T l,h on Γ f c :

         Π l,R : (H(div , R) ∩ H µ (R)) • n |Γ f c → T l,h|Γ f c q • n |Γ f c → q l,R • n |Γ f c .
With the help of a classical formula for the change of variables on Γ ([BoBF13, (2.1.62)]), one finds after a few elementary algebraic manipulations( * ) that

h d-1 K Γ (v |K • n) 2 dΓ Γ (v • n) 2 d Γ,
where Γ is equal to F -1 K (Γ ). On the reference element, it holds

v 2 H( d iv , K) Γ (v • n) 2 d Γ,
because the non-zero degrees of freedom are all located on Γ . Finally, one has the classical bounds [BoBF13, Lemma 2.1.7]:

v |K 2 0,K h 2-d K v 2 0, K div v |K 2 0,K h -d K d iv v 2 0, K , so that v |K 2 H(div ,K ) h -d K v 2 H( d iv , K) h -1 K Γ (v |K • n) 2 dΓ.
Adding up the contributions for ∈ N Γ , one finds:

δq f c H(div , R f ) h -1/2 f δq f c • n 0,Γ f c . (C.7)
By modifying the final computations in the proof of Lemma 6.2, one finds that for all 0 < < µ: δq f c • n 0,Γ f c h c|Γ f c q f,h 0,Γ f c [BGNR06, Lemma 4.9]

h c|Γ f c (h f |Γ f c ) -1/2 q f,h -1/2,Γ f c [Stei08, Lemma 10.10] h +1/2 f q f,h -1/2,Γ f c h +1/2 f Π f,R (q f • n |∂ R f ) -1/2,∂ R f h +1/2 f q f • n |∂ R f -1/2,∂ R f [BeBr01, Theorem 2.4-Remark 2.5] h +1/2 f q f , R f + div q f 0, R f .
Or, choosing = µη for η > 0 arbitrary small, that

δq f c • n 0,Γ f c h µ+1/2-η f q f µ, R f + div q f 0, R f .
Using (C.7), we conclude the proof. * Since the meshes are quasi-uniform on Γ f c , they are in particular regular. Keys words : numerical analysis, SP N equations, non-conforming.

Abstract : In this thesis, we investigate the resolution of the SPN neutron transport equations in pressurized water nuclear reactor. These equations are a generalized eigenvalue problem. In our study, we first considerate the associated source problem and after we concentrate on the eigenvalue problem. A nuclear reactor core is composed of different media: the fuel, the coolant, the neutron moderator... Due to these heterogeneities of the geometry, the solution flux can have a low-regularity. We propose the numerical analysis of its approximation with finite element method for the low regular case. For the eigenvalue problem under its mixed form, we can not rely on the theories already developed.

We propose here a new method for studying the convergence of the SPN neutron transport eigenvalue problem approximation with mixed finite element. When the solution has low-regularity, increasing the order of the method does not improve the approximation, the triangulation need to be refined near the singularities of the solution. Nuclear reactor cores are well-suited for Cartesian grids, but the refinement of these sort of triangulations increases rapidly their number of degrees of freedom. To avoid this drawback, we propose domain decomposition method which can handle globally non-conforming triangulations.

Figure 1 :

 1 Figure 1: Example of a pre-processed geometry.

Figure 2 :

 2 Figure 2: Two triangulations of the same geometry (middle), one conforming (left) and one non-conforming (right).

Figure 1

 1 Figure 1.2: Spherical coordinates of the direction vector.

  denote Y m n (ϑ, θ) := Y m n (Ω). For n, n ∈ N and m, m ∈ Z such that |m| ≤ n and |m | ≤ n , one can compute the term S 2 ΩY m n (Ω)Y m n (Ω)dΩ. We look for an expression of ΩY m n (Ω) such that the integral is a sum of scalar products of only two spherical harmonics. From the recursive relation (B.6), (B.7) and (B.8), it stands:

Figure 1 . 3 :

 13 Figure 1.3: The three main directions of the transport.

Problem 5. 30 .T

 30 For any d in L, find (y d , η d ) in Q × V such that: T o y d + grad x Hη d = 0; T Hdiv y d + T T e η d = d. (5.40) Theorem 5.31. Under condition 5.21, for any d in L, there exists a unique (y d , η d ) in X winch satisfy problem 5.29. Moreover, it stands:

Figure 6 . 1 :

 61 Figure 6.1: Two meshes of the same partition (middle), one nested non-matching (left) and one non-nested (right).

Figure 7 . 3 :

 73 Figure 7.3: Geometry for the checkboard testcase

Figure 7 . 5 :

 75 Figure 7.5: Large heavy steel reflector PWRlike Figure 7.6: PWR assembly.

Figure 7 . 7 :

 77 Figure 7.7: Subgrids use for the homogenization step.

Figure

  Figure 7.8: Large heavy steel reflector PWR-like

  introduce the bilinear form B(ζ, ξ) = B S (ζ, ξ) + B A (ζ, ξ) for ζ and ξ in X. Let ζ be the solution of Problem 8.1. One can remark, from the definition of B that for all ξ in X, it holds that:

Figure 8 . 1 :

 81 Figure 8.1: Reconstruction of the Flux

Figure 8

 8 Figure 8.4: Evolution of the number of elements during the power inverse iteration.

Figure 8 . 7 :

 87 Figure 8.7: Evolution of the number of elements during the power inverse iteration.

Figure 8 . 8 :

 88 Figure 8.8: Value of D on R.

Figure 8 . 9 :

 89 Figure 8.9: Solution of the source problem on a uniformly refined grid.

Figure 8 .

 8 Figure 8.10: Refined mesh for the source problem.

Titre:

  Analyse numérique d'une méthode de décomposition de domaine non-conforme pour les équations multigroupes SP N .Mots Clefs : analyse numérique, équations SP N , non-conforme.Résumé : Dans cette thèse, nous nous intéressons à la résolution des équations SPN du transport de neutrons au sein des coeurs de réacteurs nucléaires à eau pressurisée. Ces équations forment un problème aux valeurs propres généralisé. Dans notre étude nous commençons par le problème source associé et ensuite nous étudions le problème aux valeurs propres. Un coeur de réacteur est composé de différents milieux: le combustible, le fluide caloporteur, le modérateur... à cause de ces hétérogénéités de la géométrie, le flux solution du problème source peut être peu régulier. Nous proposons l'analyse numérique de l'approximation de la solution par la méthode des éléments finis du problème source dans le cas où la solution est peu régulière. Pour le problème aux valeurs propres, dans le cas mixte, les théories déjà développées ne s'appliquent pas. Nous proposons ici une nouvelle méthode pour étudier la convergence de la méthode des éléments finis mixtes pour les problèmes aux valeurs propres. Pour les solutions peu régulières, la montée en ordre de la méthode des éléments finis n'améliore pas l'approximation du problème, il faut raffiner le maillage aux alentours des singularités de la solution. La géométrie des coeurs de réacteur se prête bien aux maillages cartésiens, mais leur raffinement augmente vite leur nombre de degrés de liberté. Pour palier à cette augmentation, nous proposons ici une méthode de décomposition de domaine qui permet d'utiliser des maillages globalement non-conformes.Title : Numerical Analysis of a Non-Conforming Domain Decomposition for the Multigroup SP N Equations.

  G N Moments P N Couplings P N Moments SP N Couplings SP N

	1	1 3 5 7	4 16 36 64	9 69 185 357	4 8 12 16	5 11 17 23
	2	1 3 5 7	8 32 72 128	27 171 443 843	8 16 24 32	18 38 58 78
	Table 1.1: Angular moments and non null coupling terms for the SP N and P N approxi-mations in three spatial dimensions.

  From proposition 5.42 and the work of Osborn in [Osbo75, Theorem 2], one derives an estimate on the error on the eigenvalues. Corollary 5.43. Let ω be as in definition 5.41. Then for h small enough, the error on the eigenvalue is given by |νν h | h ω . Remark 5.44. If ν has an algebraic multiplicity m ν > 1, the previous analysis and the a priori estimate are still valid with ν h = 1 , where (ν h,i ) i=1,mν are the m discrete eigenvalues closest to ν, see again [Osbo75, Theorem 2].

	mν	mν i=1 ν h,i

  ) -1 C d is symmetric and positive definite, thus one can use the preconditioned conjugate gradient method (PCG) to inverse it [Shew04, Lath09], Algorithm 4. We recall that the matrices W g,g d,d are block diagonal, and each block corresponds to one subdomain. To inverse each block we use the LDL T factorization. ← R 0 |D 0 6: n pgc ← 0 7: while δ npgc > δ 0 do

	d R 0 4: D 0 ← P g 5: δ 0 8: Q npgc ← C T d (W g,g d,d ) -1 C d D npgc 9: δ npgc α npgc ← D npgc |Q npgc 10: Λ g n,ngs,nadi,npgc + αD npgc n,ngs,nadi,npgc+1 ← Λ g 11: R npgc+1 ← R npgc -αn pgc Q npgc 12:
	ngs,nadi,d -S g S,n,ngs,nadi,d ; n,ngs,nadi,d -C d Λ n,ngs,nadi ); d,d ) -1 (S g n,ngs,nadi,d = (W g P g	(7.4)
	The resolution over the domain decomposition is done in two step, The first step consists in finding the Lagrange multiplier Λ g n on the interfaces betwen the subdomains and then, the second step consists in solving the current P g n on each subdomain. The matrix C T d (W g d,d Algorithm 4 Preconditioned Conjugate Gradient Method
	1: Λ g n,ngs,nadi,0 ← Λ g n,ngs,nadi 2: S g S,n,ngs,nadi,d ← C T d (W g,g d,d ) -1 S g ,n,ngs,nadi,d -S g S,n,ngs,nadi,d 3: R 0 ← S g S,n,ngs,nadi,d -C T d (W g,g d,d ) -1 C d Λ g n,ngs,nadi,0

  table 7.3. The domain decomposition is show in Figure7.4. We refine the middle subdomain by a factor 2. The meshsize h used here is the greatest one.

	100/h 10 16 26 30 32 60 62 64 128 Rate	DoF 0.32e3 0.80e3 108 1.66e-05 9.63e-02 N out ε k eff ε φ 97 4.22e-05 1.53e-01 2.08e3 113 6.23e-06 5.95e-02 2.76e3 113 4.67e-06 5.16e-02 3.14e3 113 4.07e-06 4.84e-02 10.92e3 111 1.08e-06 2.59e-02 11.66e3 112 1.02e-06 2.51e-02 12.42e3 113 9.58e-07 2.43e-02 49.41e3 115 1.20e-07 1.25e-02 --h 2.11 h 0.99

Table 7 .

 7 1: Number of iterations and error on the criticity for the checkerboard test case with D = and RT 0 elements on a conforming mesh. Table7.2: Number of iterations and error on the criticity for the checkerboard test case with D = and RT 0 elements on a conforming mesh.

	100/h DoF N out 4 208 114 7.67e-06 7.15e-02 ε k eff ε φ 6 456 111 1.32e-06 3.21e-02 8 800 112 3.59e-07 1.84e-02 Rate --h 4.42 h 1.95

Table 7 .

 7

3: Number of iterations and error on the criticity for the checkerboard test case with D = and RT 0 elements on a non-conforming mesh.

Table 7 .

 7 4: Error on the criticity for the PWR-like test case with RT 0 elements on conforming mesh.

	410/h	DoF	SP 1	ε k eff	DoF	SP 3	ε k eff
	95	0.16e6 115.81 0.31e6 111.27
	190	0.63e6 28.91	1.24e6 26.31
	285	1.41e6 13.48	2.79e6 13.28
	380	2.49e6	8.01	4.95e6	8.39
	475	3.88e6	5.47	7.72e6	6.04
	570	5.58e6	4.10	11.11e6 4.72
	665	7.58e6	3.26	15.12e6 3.90
	760	9.90e6	2.72	19.73e6 3.38
	855	12.52e6 2.36	24.97e6	3.1
	950	15.44e6 2.09	30.82e6 2.75
	Rate	-		h 1.58	-		h 1.37

Table 7 .

 7 5: Error on the criticity for the PWR-like test case with RT 0 elements on nonconforming mesh.

  8.4. Let ζ and ζ h be respectively the solution of Problem 8.1 and Problem 8.2. Let ζh = (p h , φh ) be a reconstruction of ζ h in Q × V . For any K ∈ T h , we define the residual estimators by η

  Power Inverse Iteration With Adaptive Mesh Refinement1: initial state (P 0 , Φ 0 , Λ 0 , k eff,0 ) 2: S 0 ← M f Φ 0 3: n ← 1 4: until convergence do

	Algorithm 6 is applied to the resolution of the one-group diffusion model with RTN [0] finite elements in one dimension (d = 1) in § 8.1.4 and in two dimensions (d = 2) in § 8.1.5 ( 5:

* ).

Algorithm 6

  Table 8.1: Comparison of the power inverse iteration with and without AMR.

	0.08						1.1
	0.06						1.05
	0.04						1
	0.02						0.95
	0						0.9
	0	0.2	0.4	0.6	0.8	1	0		0.2	0.4	0.6	0.8	1
	Figure 8.2: Solution of the one-group diffusion.	Figure 8.3: Refined triangulation.
					I-AMR II-No AMR III-AMR
				N it	21		17		17
				ε k eff 7.00e -6	1.43e -4	7.00e -6

  Table 8.2: Comparison of the power inverse iteration with and without AMR.

	0.08						1.1
	0.06						1.05
	0.04						1
	0.02						0.95
	0						0.9
	0	0.2	0.4	0.6	0.8	1	0		0.2	0.4	0.6	0.8	1
	Figure 8.5: Solution of the one-group diffusion.	Figure 8.6: Refined triangulation.
	In table 8.2, we represent, as in table 8.1, the error and the number of outer iterations with AMR (I-AMR), without AMR (II-No AMR) and for the last column (III-AMR), the number of outer iteration is fixed to N it = 11, that is the number of outer iterations without AMR.
					I-AMR II-No AMR III-AMR
				N it	21		11		11
				ε k eff 2.00e -6	2.17e -4	5.1e -4

  Number of atoms of isotope i in a spatial volume around x at a given time tcm -3 σ t,i (E)Microscopic cross section of isotope i for a given energyE cm 2 σ s,i (Ω • Ω , E → E )Microscopic differential scattering cross sections of isotope i from energy E to E and direction Ω to Ω cm 2 .Mev -1 .sr -1 Total reaction rate cm -3 .Mev -1 .s -1 Microscopic scattering cross section of isotope i from energy E to E and direction Ω to Ω resulting in the formation of an isotope i cm 2 .Mev -1 .sr -1 σ * i→i ,x (E) Microscopic cross section of an event x of isotope i from energy E resulting in the formation of an isotope i

		Appendix C
		Technical Lemmas
	C.1 Results for Chapter 3
	Symbol	Meaning	Units
	C i (x, t)
	Σ t (x, E, t)	Macroscopic total cross section	cm -1
	Σ f (x, E, t)	Macroscopic fission cross section	cm -1
	Σ ν(E)	Fission yield	-
	χ(E)	Fission spectrum	Mev -1
	σ cm 2
	λ i→i	Probability of an isotope i changes into an isotope i by	s -1
		radioactive decay
	ζ i→i (E)	Microscopic reaction rate relating of the transformation	s -1
		of an isotope i into an isotope i

s (x, Ω•Ω , E → E , t) Macroscopic differential scattering cross section from energy E to E and direction Ω to Ω cm -1 τ t (x, Ω, E, t) * i→i ,s (Ω•Ω , E → E )

  Let us estimate first K∈T h ψ h 2 µ,K . According to Corollary 1.138 of[START_REF] Ern | Theory and Pratice of Finite Elements[END_REF], we know that

		ψ h	2 µ,K	h -2µ K	ψ h	2 0,K	h -2µ min ψ h	2 0,R ,	(C.3)
	where h min = min K∈T h	y)| 2 |x -y| d+2µ dydx K R |ψ h (x) -ψ h (y)| 2 |x -y| d+2µ dydx K∈T h 0,K + 2 K R\K µ,K + ψ h ψ h 2 K∈T h |ψ K∈T h = K∈T h = K∈T h

h (x)ψ h (y)| 2 |x -y| d+2µ dydx. (C.2)

† If ϕ W = 0, then νΣ f ϕ = 0. By definition of W , ϕ is solution of (2.2) with zero right-hand side. Thus, by uniqueness of the solution it follows that ϕ = 0.

* For non-nested meshes, numerical illustrations suggest that the convergence properties can be recovered in some situations (see[START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF] Table 2]). See also §6.1.2.

† f refers to fine discretization, while c refers to coarse discretization.

Part III : Chapter 7 : Numerical Applications

* The theoretical work done in section 8.1.1 and in section 8.1.2 is valid for 3 dimensions problems.

[Wang09] Y. Wang. Adaptive Mesh Refinement Solution Techniques for the Multigroup SN Transport Equation Using a Higher-Order Discontinuous Finite Element Method. Texas A&M University, 2009. [Wohl01] B. I. Wohlmuth. Discretization Methods and Iterative Solvers Based on Domain Decomposition, volume 17 of Lecture Notes in Computational Science and Engineering. Springer-Verlag, 2001.

follows that B µ is a compact operator. Denote by σ(B µ ) its spectrum. By construction, λ -1 ∈ σ(B µ ) if, and only if, λ is an eigenvalue of (2.2). Finally, we consider the discrete operator B h µ associated to the discrete source problem (5.2): given f ∈ H µ (R), we call B h µ f the second component of the couple (p h , φ h ) that solves (2.17) with source S f = νΣ f f . Under the assumptions of §5.1.1 and as noted at the beginning of §5.1.3, it holds for all f ∈ L, lim h→0 B 0 f -B h 0 f L = 0 . This property is the so-called pointwise convergence. However, for a mixed formulation, the fact that the family (B h 0 ) h converges pointwise towards the compact operator B 0 is not sufficient to guarantee that the family (B h 0 ) h converges in operator norm towards B 0 .

Convergence in Operator Norm

On the other hand, according to [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF], proving that lim h→0 B µ -B h µ L(H µ (R)) = 0 for discrete approximants (B h µ ) h is a sufficient condition to obtain convergence of the eigenvalues. In order to ensure the convergence in operator norm of the family (B h µ ) h towards the compact operator B µ , we need a technical assumption on the triangulations. Definition 5.8. A family of triangulations (T h ) h is regular + if it satisfies: ∃θ > 0, ∀h, h 2-θ min K∈T h diam(K).

(5.18)

In particular, a quasi-uniform family of triangulations is regular + (take θ = 1 in (5.18)).

For a regular + family, one has the following inverse inequality, whose proof is given in the Appendix C.

Lemma 5.9. Let µ ∈ [0, 1/2[. For a regular + family of triangulations, it holds: ∀h, ∀ψ h ∈ L k h , ψ h µ,R h -2µ+θµ ψ h 0,R .

(5.19)

We note that the hidden constant in equation (5.19) depends on k.

Theorem 5.10. Under the assumptions of Proposition 2.5 with r max < 1/2 plus νΣ f ∈ PW 1,∞ (R), let µ ∈ [0, r max [. Provided that the family of triangulations is regular + , one has:

h θµ .

(5.20)

Proof. According to (5.15), we know that

(5.21)

It remains to estimate (B µ -B h µ )f µ,R : for that, we use the triangle inequality

To bound the first term, we have according to Theorem 2.3 in [START_REF] Belgacem | Some Nonstandard Finite Element Estimates with Applications to 3D Poisson and Signorini Problems[END_REF] that ∀ψ ∈ PH 1 (R), ψπ 0 ψ µ,R h 1-µ ψ PH 1 (R) .

Applying the result to ψ = B µ f , we find B µ fπ 0 (B µ f ) µ,R h 1-µ f µ,R . To bound the second term, we use first the inverse inequality (5.19) on the discrete space L k h , valid for a regular + family of triangulations. Applying the result to ψ h = π 0 (B µ f ) -B h µ f and using again the triangle inequality, we now find that

where we have used (5.4) and (5.21) to derive the final estimate. Since µ < 1/2, we conclude by aggregating the results that (5.20) holds.

Thanks to [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF], convergence of the discrete eigenvalues to the exact ones is guaranteed, and so is the absence of spectral pollution (see § 4.3).

Optimal Convergence Rate

Let the assumptions of Theorem 5.10 hold. We determine now the rate of convergence of the eigenvalues in the spirit of [START_REF] Boffi | Optimal Convergence of Adaptive FEM for Eigenvalue Clusters in Mixed Form[END_REF]. Let ν = λ -1 be an eigenvalue of B µ . For simplicity, let us assume that ν is a simple eigenvalue, and denote by W the associated eigenspace. According to the absence of spectral pollution, for h small enough, the closest discrete eigenvalue, denoted by ν h , is also simple ; we denote by W h the associated eigenspace.

Definition 5.11. Let ω ν > 0 be the regularity exponent of the eigenfunction, i.e. either W ⊂ PH 1+s (R) for s < ω ν and W ⊂ PH 1+ων (R), or W ⊂ PH 1+ων (R) and W ⊂ PH 1+s (R) for s > ω ν . Let ω = min(ω ν , m + 1), where m ≥ 0 is the order of the RTN finite element.

Clearly, ω ν , and as a consequence ω, can be greater than r max . We shall prove that the approximation converges with a rate equal to twice the exponent ω defined above: this result is stated in Corollary 5.19 at the end of the subsection. Let µ ∈ [0, r max [ be given. As we defined B µ (resp. B h µ ), we define A µ (resp.

The following lemma introduces some equalities that we will use later on.

Lemma 5.12. Let ϕ and ϕ be given in W . Then, it holds:

(5.22)

(5.23)

Proof. The definitions of A µ , B µ imply that for all f ∈ H µ (R), for all (q, ψ) ∈ X:

whereas the definitions of A h µ , B h µ imply that for all f ∈ H µ (R), for all (q, ψ) ∈ X h :

(5.25)

The first equality (5.22) comes from (5.24) with:

The second one, (5.23), comes from the difference between (5.24) and (5.25) with: f = ϕ ; q = A h µ ϕ ; ψ = B h µ ϕ ; and with the symmetry of a(•, •) and t(•, •).

We remark that ϕ → ϕ W = (νΣ f ) 1 2 ϕ 0,R is a norm over W ( † ), and this norm is induced by the inner product

Proposition 5.13. Let ω be as in definition 5.11. For every ϕ in W , the following inequalities hold:

Proof. These two inequalities come from the first Strang's Lemma. The method is the same as for Theorem 5.3 (see remark 5.4 for the "smooth" case). Here, we use the equivalence of all norms on W to state the result.

Introducing δ(Z, Z ) = sup z∈Z, z 0 =1 inf z ∈Z zz 0,R for Z, Z closed subspaces of L, the gap between W and W h is defined by:

It allows us to evaluate the approximation of the continuous eigenfunctions by their discrete counterpart. Classically, this gap can be bounded with the help of Proposition 5.13, following [Osbo75, Theorem 1]:

(5.26)

Let us now define E h as the projector from L onto W h such that

(5.27)

Lemma 5.14. The operators

. This is equivalently expressed as

The penultimate equality stems from the fact that c(•, •) is symmetric, and the last one comes from the definition of φ and E h .

Let us bound each part separately below. One obtains from the difference between (5.22) and (5.23)

). Then, one can bound the first part:

The second part is bounded by:

In the first line we use Proposition 5.17 with f = (B µ -B h µ )ϕ and ϕ h = F h ϕ. In the third line we use the uniform continuity of S h in h, and in the last line we use the first inequality of Proposition 5.13 and the estimation (5.28). Therefore we have obtained (5.29).

From this estimation and the work of Osborn in [Osbo75, Theorem 2], one derives an optimal estimate on the error on the eigenvalues.

Corollary 5.19. Let ω be as in definition 5.11. Then for h small enough, the error on the eigenvalue is given by |ν -

Remark 5.20. If ν has an algebraic multiplicity m ν > 1, the previous analysis and the a priori estimate are still valid with ν h = 1 mν mν i=1 ν h,i , where (ν h,i ) i=1,mν are the m discrete eigenvalues closest to ν, see again [Osbo75, Theorem 2].

Extension to the Multigroup SP N Transport Equations

We extend here the results of § 5.1 on the approximation of the diffusion equations with mixed finite element method. We recall that the function spaces used for the SP N equations under the mixed setting are:

We also recall that H s stands for

We recall the condition 1.13 on page 27:

Numerical Analysis of the Generalized Eigenvalue Problem

Let us focus on the approximation of the generalized eigenvalue problem 3.2 in our lowregularity setting, under condition 5.21 and hypothesis 5.22, supplemented with condition 3.18. Thanks to condition 5.21 and condition 3.18, M f is a continuous operator from H µ to H µ , µ < 1/2. Let 0 ≤ µ < r max be given,we recall from § 3.3.2 that the operator B µ associated to the source problem 3.15, given f ∈ H µ , we call

Finally, we consider the discrete operator B h µ associated to the discrete source problem 5.24: given f ∈ H µ , we call B h µ f the second component of the couple (p h , φ h ) that solves (5.33) with source S f = M f f . Under 5.21 and hypothesis 5.22, theorem 5.37 holds so that for all f ∈ L, lim h→0 B 0 f -

This property is the so-called pointwise convergence. However, for a mixed formulation, the fact that the family (B h 0 ) h converges pointwise towards the compact operator B 0 is not sufficient to guarantee that the family (B h 0 ) h converges in operator norm towards B 0 . On the other hand, according to [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF], proving that lim h→0 B µ -B h µ L(H µ ) = 0 for discrete approximants (B h µ ) h is a sufficient condition to obtain convergence of the eigenvalues. In order to ensure the convergence in operator norm of the family (B h µ ) h towards the compact operator B µ , as in § 5.1.4, we need to suppose that the triangulation is regular + , see definition 5.18 on page 67. The inverse inequality in lemma 5.9 can be extended onto L k h .

Lemma 5.39. Let µ ∈ [0, 1/2[. For a regular + family of triangulations, it holds:

(5.49)

Theorem 5.40. Under condition 5.21, hypothesis 5.22 with r max < 1/2 and condition 3.18, for µ ∈ [0, r max [, provided that the family of triangulations is regular + , one has:

(5.50)

Proof. The proof is the same as the one for theorem 5.10.

Thanks to [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF], convergence of the discrete eigenvalues to the exact ones is guaranteed, and so is the absence of spectral pollution (see § 4.3). Moreover, we can derive an a priori error estimate on the eigenvalue. Let the assumptions of theorem 5.40 hold. We determine now a rate of convergence of the eigenvalues using the work of Osborn in [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF]. Let ν = λ -1 be an eigenvalue of B µ . For simplicity, let us assume that ν is a simple eigenvalue, and denote by W the associated eigenspace. According to the absence of spectral pollution, for h small enough, the closest discrete eigenvalue, denoted by ν h , is also simple ; we denote by W h the associated eigenspace.

Definition 5.41. Let ω ν > 0 be the regularity exponent of the eigenfunction, i.e. either W ⊂ PH 1+s for s < ω ν and W ⊂ PH 1+ων , or W ⊂ PH 1+ων and W ⊂ PH 1+s for s > ω ν .

Let ω = min(ω ν , m + 1), where m ≥ 0 is the order of the RTN finite element.

Again, this bound holds under the condition that the meshes are quasi-uniform on the part of the interface where they are non-matching. Due to (5.9), one has

As a consequence (follow § 5.1.3) we conclude that the estimate (6.18) holds.

Numerical Analysis of the Generalized Eigenvalue Problem

Let us focus on the approximation of the generalized eigenvalue problem (2.2) for lowregularity solutions with nested (matching or non-matching) meshes. We will follow the methodology of § 5.1.4.

Convergence in Operator Norm

Let 0 ≤ µ < r max be given, we introduce an operator B µ associated to the source problem (6.6):

the second component of the triple (p, φ, φ S ) that solves the source problem with S f = νΣ f f . For the same reason as in § 5.1.4, B µ is a bounded and compact operator. Next, let us consider the discrete operator B h µ associated to the discrete source problem: given f ∈ H µ (R), we call B h µ f the second component of the triple (p h , φ h , φ S,h ) that solves (6.10) with source S f = νΣ f f . Using estimate (6.18), we obtain, like in Chapter 5, the result below. Theorem 6.6. Under the assumptions of Theorem 6.3 with r max < 1/2 plus νΣ f ∈ PW 1,∞ (R), let µ ∈]0, r max [. Provided that the families of triangulations are regular + on every subdomain, one has for nested meshes:

where θ = min N i=1 θ i > 0, and for 1 ≤ i ≤ N , θ i is defined by (5.18) on R i . We conclude to the absence of spectral pollution.

Optimal Convergence Rate

Let the assumptions of Theorems 6.3 and 6.6 hold, and in particular the conditions for nested, non-matching meshes. We use the same notations as in § 5.1.4. In particular, let ω ν > 0 be the regularity exponent associated to ν with respect to ( PH 1+s (R)) s>0 , and introduce ω = min( ω ν , k + 1). Let µ ∈ [0, r max [ be given. As we defined B µ (resp. B h µ ), we define A µ and C µ (resp.

) the first and the third components of the triple (p, φ, φ S ) (resp. (p h , φ h , φ S,h )) that solves (6.6) (resp. (6.10)) with source S f = νΣ f f . For the DD+L 2 -jumps method, the transposition of Lemma 5.12 reads: Lemma 6.7. Let ϕ and ϕ be in W . Then, it holds:

(6.20)

(6.21)

Proof. The definitions of A µ , B µ and C µ imply that for all f ∈ H µ (R) and for all (q, ψ, ψ S ) ∈ X:

where the penultimate term

Whereas the definitions of A h µ , B h µ and C h µ imply that for all f ∈ H µ (R), for all (q, ψ, ψ S ) ∈ X h :

The first equality (6.20) comes from (6.22) with:

The second equality (6.21), comes from the difference between (6.22) and (6.23) with:

and with the symmetry of a and t. The formulas (6.20) and (5.22), resp. (6.21) and (5.23), are identical. As Strang's Lemma holds for the DD+L 2 -jumps method with nested meshes, we can also transpose Proposition 5.13. For that, we admit that the result of Lemma 6.2 can be improved for smooth functions q. As a matter of fact, in this case one may directly compare the discrete normal traces Π f,R (q•n |Γ f c ) and Π c,R (q•n |Γ f c ) to the exact normal trace q•n |Γ f c , and evaluate the difference in L 2 (Γ f c )-norm, because for smooth functions the exact normal trace always belongs to L 2 (Γ f c ). Proposition 6.8. For every ϕ in W , the following inequalities hold for the DD+L 2 -jumps method with nested meshes:

Estimate (5.26) on the gap between W and W h is still valid: δ(W, W h ) h ω . Let E h be the operator defined in (5.27). We recall that E h and B h µ commute (Lemma 5.14 holds). The restriction of E h to W , denoted by F h is a bijection that satisfies estimate (5.28), for h small enough. We will also make use of S h = F -1 h E h -I that satisfies Lemma 5.16 and Proposition 5.17. We recall that Bµ = B µ | W and Bh µ = F -1 h B h µ F h . The transposition of Theorem 5.18 is stated next. The proof is identical (replace ω by ω), so it is omitted. Theorem 6.9. For h small enough, one has for the DD+L 2 -jumps method with nested meshes:

Corollary 6.10. For h small enough, the error on the eigenvalue for the DD+L 2 -jumps method with nested meshes is given by:

Problem 6.14. Find u ∈ W such that ∀w ∈ W: c s,S (u, w) = f s,S (w). (6.31) Theorem 6.15. Under condition 6.11 and hypothesis 5.22, there exists a unique solution u ∈ W which satisfies problem 6.14.

Proof. Adapting the results of [START_REF] Ciarlet | Domain Decomposition Methods for the Diffusion Equation with Low-Regularity Solution[END_REF] for the neutron diffusion equations to the multigroup SP N transport equations, one can prove that c s,S satisfies an inf-sup condition. As the problem is not symmetric, c s,S has to satisfies the solvability condition. Using the same idea of the proof of theorem 3.16 with:

one can prove that ξ = 0. Then taking p a lifting of ψ S in Q, one can conclude the proof.

Therefore, problem 6.14 is well-posed.

In order to study abstract conforming discretization of the variational problem 6.14, we introduce the following finite dimension spaces:

We suppose that the mesh is the same for each energy group and for each moment. The discrete counterpart of the variational problem 6.14 reads: Problem 6.16. Find u h ∈ W h such that ∀w h ∈ W h : c s,S (u h , w h ) = f s,S (w h ). (6.32)

We consider problem 6.16 where the RTN finite element is used on each subdomain with a conforming mesh, or triangulation.

Theorem 6.17. Under condition 6.11 and hypothesis 5.22, there exists a unique solution u h ∈ W h which satisfies problem 6.16.

Proof. As (6.11)-(6.12) holds for each energy group and for each moment, one can prove that these conditions are satisfied on Q h and M h . Therefore c s,S satisfies a udisc.

Therefore, problem 6.16 is well-posed. Moreover, according to the first Strang's Lemma [START_REF] Ern | Theory and Pratice of Finite Elements[END_REF] and because c s,S verifies a udisc, the error reads:

As a consequence lim h→0 uu h W = 0. This result holds for nested and non-nested meshes. We study below how to improve the bound on the error, how to derive an Aubin-Nitsche estimate, and finally how to prove convergence for the generalized eigenvalue problem, for nested meshes. Those results hold under condition 6.11 on the coefficients and hypothesis 5.22 on the existence of a regularity exponent (plus νΣ g f ∈ PW 1,∞ (R), 1 ≤ g ≤ G, for the eigenproblem). We focus again on the low-regularity case.

of the discrete with DDM problem 6.16 is such that (p h , φ h ) ∈ X h , since p h ∈ Q h . Then restricting the test-fields in problem 6.16 to elements of X h × M h we observe that (p h , φ h ) satisfies the discrete without DDM problem 5.24 too, because all interface terms vanish. Hence, to estimate φφ h L in the DDM, we explicitly consider that the discrete fields (p h , φ h ) are also the solution to the variational formulation problem 5.24. Let us begin by a technical result. Lemma 6.21. Let the assumptions of lemma 6.19 hold. Let q ∈ Q∩H µ with 0 < µ < 1/2, and define

Proof. The proof is done in the same way as in the proof of proposition 5.26 using lemma 6.4.

Theorem 6.22. Under condition 6.11 and hypothesis 5.22 with r max < 1/2, one has for nested meshes:

Proof. Matching meshes. In this case, one can use the theory developed in § 5.2.1 to conclude that (6.35) holds. Nested meshes. From the proof of theorem 6.5, we know that for all energy groups g and all moments n, it stands:

where (p g n ) R is defined in (6.1.2). Defining the Q-interpolant p R of p by (p R ) g n = (p g n ) R , for all energy groups g and all moments n, one can prove that:

Then, by following the methodology of the proof of theorem 5.37, using proposition 5.26, lemma 6.21 and inequalities (5.35), one can conclude that estimate (6.35) holds.

Numerical Analysis of the Generalized Eigenvalue Problem

Let us focus on the approximation of the generalized eigenvalue problem 3.1 for lowregularity solutions with nested (matching or non-matching) meshes. We will follow the methodology of § 5.1.4. Let 0 ≤ µ < r max be given, we introduce an operator B µ associated to the source problem 6.14 given f ∈ H µ , we call B µ f = φ ∈ V the second component of the triple (p, φ, φ S ) that solves the source problem 6.14 with S f = M f f . For the same reason as in § 5.2.2, B µ is a bounded and compact operator. Next, let us consider the discrete operator B h µ associated to the discrete source problem 6.16: given f ∈ H µ , we call B h µ f the second component of the triple (p h , φ h , φ S,h ) that solves problem 6.16 with source S f = M f f . Using estimate (6.35), we obtain, like in section 5.2.2, the result below. Theorem 6.23. Under condition 6.11 and hypothesis 5.22 with r max < 1/2 plus condition 3.18, let µ ∈]0, r max [. Provided that the families of triangulations are regular + on every subdomain, one has for nested meshes:

where θ = min N i=1 θ i > 0, and for 1 ≤ i ≤ N , θ i is defined by (5.18) on R i . We conclude to the absence of spectral pollution. Moreover, we can derive an a priori error estimate on the eigenvalue. Let the assumptions of theorem 6.23 hold. We determine now a rate of convergence of the eigenvalues using the work of Osborn in [START_REF] Osborn | Spectral Approximation for Compact Operators[END_REF]. Let ν = λ -1 be an eigenvalue of B µ . For simplicity, let us assume that ν is a simple eigenvalue, and denote by W the associated eigenspace. According to the absence of spectral pollution, for h small enough, the closest discrete eigenvalue, denoted by ν h , is also simple ; we denote by W h the associated eigenspace. Definition 6.24. Let ω ν > 0 be the regularity exponent of the eigenfunction, i.e. either W ⊂ PH 1+s for s < ω ν and W ⊂ PH 1+ ων , or W ⊂ PH 1+ ων and W ⊂ PH 1+s for s > ω ν . Let ω = min(ω ν , m + 1), where m ≥ 0 is the order of the RTN finite element.

We remark that ϕ 2 W = (M f ϕ, ϕ) L is a norm over W (see the definition of • W in § 5.1). However, this norm does not derive from a inner product since M f is not symmetric. For this reason, we cannot use the same methodology as in § 5.1.4. Proposition 6.25. Let ω be as in definition 6.24. For every ϕ in W , the following inequalities hold:

Proof. The inequality comes from the first Strang's lemma. The method is the same as for theorem 5.3. Here, we use the equivalence of all norms on W to state the result.

From proposition 6.25 and the work of Osborn in [Osbo75, Theorem 2], one derives an estimate on the error on the eigenvalues.

Corollary 6.26. Let ω be as in definition 6.24. Then for h small enough, the error on the eigenvalue is given by |νν h | h ω .

Remark Here we give the list of all the physical quantities linked to the neutron density in the reactor core.

Symbol Meaning Units

N (x, Ω, E, t) Neutron density in a phase volume dxdΩdE at time t cm -3 .Mev -1 .sr -1 ψ(x, Ω, E, t)

A.3 List of Nuclear Data

Note that the microscopic cross sections are physical data whereas macroscopic cross sections are data resulting from the homogenization process.

Appendix B Spherical Harmonics Functions

We present here some properties on the Legendre polynomials, the associated Legendre polynomials and the spherical harmonics [START_REF] Hochstadt | The Functions of Mathematical Physics[END_REF][START_REF] Müller | Spherical Harmonics[END_REF]. As we show some recursive formulas, we denote for n ∈ Z:

B.1 Legendre Polynomials

The n th Legendre polynomial P n is defined as the solution on [-1; 1] of the Legendre's differential equation

Those polynomials are orthogonal and it stands for any n and m integers 1 -1

The Legendre polynomials satisfy also a recursive relation given by:

B.2 Associated Legendre Polynomials

To every single Legendre's polynomial P n , one can associated 2n + 1 polynomials defined as for m ∈ [-n; n] integer and x ∈ [-1; 1]

B.3 : Normalized Spherical Harmonics

Moreover, the following recursive formulas stand One can develop Legendre's polynomials over the spherical harmonics thank to the addition theorem

Another version of the addition theorem is given by:

From the third and the fourth recursive relations in (B.2), one can deduce that With a slight abuse of notations, we write Π l,R (q l • n |∂ R l ) = q l,R • n |∂ R l . We also introduce the operator Π 0 c,R on the vector space of normal traces of elements of Q c,h with lowest-order RTN finite element, i.e. the vector space T 0 c,h|Γ f c of piecewise constant functions on the interface mesh defined as the trace on Γ f c of the mesh used in R c . Note that because the meshes are nested, the restriction of Π f,R (resp., Π c,R and Π 0 c,R ) on T f,h|Γ f c (resp., on the subspaces T c,h|Γ f c and T 0 c,h|Γ f c where applicable) may also be considered as an orthogonal projection operator. Denoting q f,h = Π f,R (q • n |Γ f c ), we have:

As the meshes are quasi-uniform on the interface, one has h c|Γ f c h f |Γ f c . Then, starting from (C.6), thanks to the quasi-uniform mesh assumption for the inverse inequalities on Γ f c , cf. [Stei08, Lemma 10.10], we find

Above, we have used the continuity of the normal trace, resp. the stability of the RTN interpolant, to derive the last two inequalities.

Then, one can bound v |K H(div ,K ) by v |K • n 0,Γ for each index ∈ N Γ . To that aim, one goes back to the reference element K via the Piola transform, which reads [BoBF13, §2.1.3]: