J.-J Lautard † ‡

Y Maday § ¶ ‡

O Mula

MINARET OR THE QUEST TOWARDS THE USE OF TIME-DEPENDENT NEUTRON TRANSPORT SOLVERS FOR NUCLEAR CORE CALCULATIONS ON A REGULAR BASIS

Keywords: MINARET, neutronics, neutron transport, parareal in time algorithm, time domain decomposition, core calculations AMS subject classifications

 (developed at CEA in the framework of the APOLLO3 R project), we aim first of all at presenting a brief and comprehensive overview of the most widespread resolution techniques employed nowadays in neutron transport industrial codes. Given that the main obstacle in the use of this type of accurate solver on a regular basis relies in the long computing times, MINARET has been used in the present work as a support to rigorously quantify the efficiency of the most common sequential and parallel acceleration techniques that are currently used in this field. An important part of the paper will be devoted to study the performances of an acceleration method that has never been considered before in the resolution of this equation, which is the parallelization of the time variable. In this regard, the parareal in time algorithm (a domain decomposition method for the time variable, [20]) has been implemented to explore its potentialities in this particular application.

1. Introduction. In the framework of numerical simulations, the advances of computing architectures in the last decades have resulted in a progressive replacement of traditional coarse models by ever-increasing finer ones. In the field of nuclear core calculations and, more particularly, regarding the resolution of the time-dependent neutron transport equation, this development of the computing capabilities has led to significant advances. Indeed, codes like PARTISN [START_REF] Alcouffe | Time-dependent neutron transport on parallel architectures using PARTISN[END_REF], TORT-TD [START_REF] Saubert | The 3-D timedependent transport code TORT-TD and its coupling with the 3D thermal-hydraulic code ATTICA3D for HTGR applications[END_REF] and the recently developped time-dependent module in MINARET (that will be presented in this paper) confirm this fact and have proven that the traditional diffusion [START_REF] Chauvet | A multi-scale approach for the neutronic kinetics equation using the mixed dual solver minos[END_REF], improved quasi-static [START_REF] Girardi | Comparison of direct and quasistatic methods for neutron kinetic calculations with the EDF R&D COCAGNE code[END_REF] or point kinetics [START_REF] Reuss | Précis de neutronique[END_REF] traditional approximations can now start to be overcome. While this represents a significant progress for security assessment (approximations like diffusion are indeed not accurate enough in the study of fast transients or fast reactors), the long computing times of transport solvers in realistic 3D core geometries represent the main obstacle to face in order to definitely use this type of accurate solvers on a regular basis. In this framework, the present work is a contribution towards showing that these computing times can be very significantly reduced through acceleration methods that do not alter the main structure of the solvers. Towards this end, it will be explained how the traditional (sequential and parallel) acceleration methods employed in this field can be coupled efficiently and in a "non-intrusive" manner with other innovative acceleration methods that involve modern parallel architectures. For this reason, a non negligible part of the paper will be devoted to recall the current construction of these solvers as well as the classical acceleration strategies in order to better justify how new techniques can be added without modifying their main construction. A special effort has been done to provide tractable results about the performances of all the methods that will be discussed. These come from computations obtained with MINARET in a well-known benchmark test case of the literature [START_REF] Langenbuch | Coarse-mesh flux expansion method for the analysis of space-time effects in large light water reactor cores[END_REF]. The results will illustrate at the same time the work that has been done to optimize the implementation of MINARET.

Hence, the paper is organized as follows: we will start by recalling in section 2 the equation under consideration. Then, in section 3, we will present an overview of the most widespread discretization and resolution strategies that are employed in the field of nuclear safety calculations. This will be done through the particular example of the newly developed time-depend module of the discrete ordinates MINARET solver (we refer also to other codes like PARTISN and TORT-TD for similar approaches). Then, section 4 will be devoted to a simple validation study of MINARET in order to prove the validity of the results presented in the following sections.

We will after discuss the sequential and parallel acceleration techniques that are classically employed in these codes or that are currently the subject of active research and that do not require changes in the resolution scheme introduced in section 3. Section 5 will recall the two most traditional sequential accelerations (the Chebyshev extrapolation and the Diffusion Synthetic Acceleration). The combination of both strategies leads to speed-ups of a factor of around 100 as computations with MINARET will show. We will also explain how extra speed-ups can be obtained by the choice of an appropriate starting guess in the so-called outer iterations. In particular, the idea of involving diffusion computations in the starting guess had never been tested in the literature and seems to outperform the classical choices (a speed-up of a factor of about 3 has been noted).

The remaining part of the paper is devoted to parallel acceleration techniques. In section 6, we recall the usual method to parallelize the angular variable in discrete ordinates codes and illustrate its performances. The limitations of this strategy have motivated the exploration of other sources of concurrency in the litterature such as the parallelization of the spatial and energetic variables. Regarding the spatial variable, the so-called KBA algorithm was proposed in the late nineties (see [START_REF] Baker | An Sn algorithm for the massively parallel CM-200 computer[END_REF]) and seems to be the most efficient strategy known so far. Its efficient implementation on unstructured meshes represents one of the main difficulties to face as the works of [START_REF] Pautz | An algorithm for parallel S N sweeps on unstructured meshes[END_REF] show. Also, the parallelization of the energy variable seems to be a promising field to explore as outlined in [START_REF] Slaybaugh | Rayleigh quotient iteration in 3d, deterministic neutron transport[END_REF] and some recent contributions such as [START_REF] Adams | Provably optimal parallel transport sweeps on regular grids[END_REF] even aim at efficiently parallelizing jointly some of these three variables (energy, space and angle). However, among all these already existing approaches, the parallelization of the time variable seems to be missing and, for this reason, in section 7, a study about the additional speed-up that can be obtained with such a strategy will be presented. This task has been performed by a domain decomposition method called the parareal in time algorithm (see [START_REF] Lions | Résolution d'EDP par un schéma en temps pararéel[END_REF] [6] [START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF] for its theoretical foundations). The implementation is following some preliminary analysis made in [START_REF] Baudron | Parareal in time 3D numerical solver for the LWR Benchmark neutron diffusion transient model[END_REF] [START_REF] Baudron | The parareal in time algorithm applied to the kinetic neutron diffusion equation[END_REF] where the strategy was successfully applied to the time-dependent neutron diffusion equation. We will propose an extension of the classical theoretical expected speed-up formulae from the literature (see, e.g. [START_REF] Aubanel | Scheduling of tasks in the parareal algorithm[END_REF] for a summary on this topic) in order to take into account the fact that, in the present case, the parareal algorithm is being coupled with other iterative methods. The limitations regarding the maximum attainable speed-up will be discussed at this point. Despite this fact, a comparison between measured speed-up performances with MINARET and the theoretical optimal values will highlight the negligible impact of the communication time between processors and will illustrate the capacity of the method to accelerate long time transients. We will end up by outlining a strategy that is currently being developed to improve the performances of the parareal algorithm by coupling the method with the outer iterations in a more efficient manner. The results will be the published on a sequel of this paper.

2. The time-dependent neutron transport equation. The evolution of the angular flux ψ of neutrons in a reactor core R is governed by a linear Bolzman equation whose terms physically express a balance between the free neutrons that are created and that disappear in the core. We will consider here the three-dimensional case (R ⊂ R 3 ) where ψ depends on 7 variables, namely the time t ∈ [0, T ], the position within the reactor denoted as r ∈ R, the velocity of the neutrons v = 2E/m ω where E ∈ [E min , E max ] stands for the energy of the neutron, ω = v |v| ∈ S 2 stands for the direction of the velocity and m is the mass of the neutron. We will have v ∈ V, where V = S 2 ×[E min , E max ] is a compact subset of R 3 . The fission chain reaction that takes place inside the core leads to the presence of some radioactive isotopes that emit neutrons with a given delay (we refer to them as precursors of delayed neutrons). This phenomenon must be taken into account in our balance equation, hence the coupling of the Bolzman equation with a set of first order ODE's expressing the evolution in R of the precursors' concentration that will be denoted as C = {C ℓ } ℓ∈{1,...,L} .

The set (ψ, C) is thus the solution to the following initial value problem over the domain

D = {(t, r, ω, E) ∈ [0, T ] × R × S 2 × [E min , E max ]}: (2.1)            1 |v| ∂ t ψ(t, r, ω, E) + (L -H -F -Q)ψ(t, r, ω, E) = 0 ∂ t C ℓ (t, r) = -λ ℓ C ℓ (t, r) + Emax E ′ =Emin β ℓ (t, r, E ′ )(νσ f )(t, r, E ′ )φ(t, r, E ′ )dE ′ , ∀ℓ ∈ {1, . . . , L},
where φ(t, r, E) = S2 ψ(t, r, ω ′ , E)dω ′ is the scalar flux and the following operator notations have been used:

• Lψ(t, r, ω, E) = (ω.∇ + σ t (t, r, E)) ψ(t, r, ω, E) is the advection operator, • Hψ(t, r, ω, E) = S2 Emax E ′ =Emin σ s (t, r, ω ′ → ω, E ′ → E)ψ(t, r, ω ′ , E ′ )dE ′ dω ′
is the scattering operator,

• F ψ(t, r, ω, E) = χ p (t, r, E) 4π Emax E ′ =Emin (1 -β(t, r, E ′ )) (νσ f )(t, r, E ′ )φ(t, r, E ′ )dE ′
is the prompt fission operator,

• Qψ(t, r, ω, E) = L ℓ=1 λ ℓ χ d,ℓ (t, r, E)C ℓ (t, r
) is the delayed fission source. Note that this term depends on ψ because C ℓ depends on it in vertue of the ODE in time that they satisfy. In the enlisted terms:

• σ t (t, r, E) denotes the total cross-section,

• σ s (t, r, ω ′ → ω, E ′ → E)
is the scattering cross-section from energy E ′ and direction ω ′ to energy E and direction ω, • σ f (t, r, E) is the fission cross-section, • ν(t, r, E) is the average number of neutrons emitted, • χ p (t, r, E) and χ d,ℓ (t, r, E) are respectively the prompt spectrum and the delayed spectrum of precursor ℓ, • λ ℓ and β ℓ (t, r, E) are the decay constant and the delayed neutron fraction of precursor ℓ respectively,

• β(t, r, E) = L ℓ=1 β ℓ (t, r, E).
We will work with initial conditions ψ 0 and C ℓ,0 at t = 0 and vacuum boundary conditions over ∂R, i.e.

     ψ = 0, on [0, T ] × ∂R -× S 2 × R + ψ(0, .) = ψ 0 (.), on R × S 2 × R + C ℓ (0, .) = C ℓ,0 (.) on R,
where ∂R -:= {r ∈ ∂R | ω. -→ n < 0} denotes the part of the boundary where the angular flux is incoming. The knowledge of the initial conditions ψ 0 and C ℓ,0 is a complex issue in itself. In nuclear safety computations like the ones we are interested in, it is of special interest to analyze transients starting from a stable state of the core. The derivation of this state is related to an eigenvalue problem whose foundations are very well established. We refer to [START_REF] Reuss | Précis de neutronique[END_REF] and [START_REF] Dautray | Analyse mathématique et calcul numérique[END_REF] for physical and mathematical aspects of it and to [START_REF] Moller | Éléments finis courbes et accélération pour le transport de neutrons[END_REF] for numerical details about its computation in the MINARET solver. Remark 2.1 (A diffusion problem as an approximation to the transport equation (2.1)). Under some given physical hypothesis (see, e.g. [START_REF] Reuss | Précis de neutronique[END_REF] and [START_REF] Dautray | Analyse mathématique et calcul numérique[END_REF]), the angular mean value φ(t, r, E) = S2 ψ(t, r, ω ′ , E)dω ′ of the angular flux ψ(t, r, ω, E) satisfies a diffusion equation that has the advantage of being much less computationally expensive than the transport equation from the memory storage and from the computational time point of view. Although the present work deals with the resolution of the transport equation (2.1), the existence of such surrogate approximation will be used in our case in some acceleration techniques.

Classical discretization and implementation schemes in neutronics:

an example through the MINARET solver. With the exception of some simple cases (see [START_REF] Reuss | Précis de neutronique[END_REF] for further references) where problem (2.1) can exactly be solved, the resolution of (2.1) needs to be numerically addressed and requires discretizations and approximations of the involved variables. We will present here the classical discretizations and numerical schemes that are used in the field of core calculations. This will be done through the particular case of the MINARET solver, whose kinetic module has never been presented before, but that can be taken as a representative example due to its similarities with other solvers like PARTISN or TORT-TD.

We start by discretizing the energy variable and deriving the multigroup version of equation (2.1). The strategy is based on the division of the energy interval into

G subintervals: [E min , E max ] = [E G , E G-1 ] ∪ • • • ∪ [E 1 , E 0 ]. For 1 ≤ g ≤ G,
we denote by ψ g the approximation of ψ in the subinterval [E g , E g-1 ]. Let us introduce a series of time steps δt n > 0 and a serie of times t n such that t 0 = 0, t N = T and t n -t n-1 = δt n . An Euler-backward scheme for the time variable is then applied. Let ψ g,n (r, ω) be the approximation of ψ(t, r, ω, E) at time t = t n and for E ∈ [E g , E g-1 ]. Given {ψ g,n (r, ω)} G g=1 , the set of unknowns {ψ g,n+1 (r, ω)} G g=1 for the time t n+1 is the solution of the following set of coupled source problems:

(3.1)
Find over R × S 2 the angular flux ψ g,n+1 (r, ω) such that

L g -H g -F g -Qg ψ g,n+1 (r, ω) = Sg,n (r, ω), ∀g ∈ {1, . . . , G}.
The following notations have been used:

• Sg,n (r, ω) := ψ g,n (r, ω) V g δt n+1
, where V g is the average velocity of the neutrons whose energy belong to the interval [E g , E g-1 ]. Note that for the computation of ψ g,n+1 (r, ω), the term Sg,n (r, ω) is known and is a source for the equation.

• L g ψ g,n+1 (r, ω) = ω.∇ + σ g,n+1 t (r) + 1 V g δt n+1 ψ g,n+1 (r, ω) • H g ψ g,n+1 (r, ω) = G g ′ =1
H g ′ →g ψ g ′ ,n+1 (r, ω), with

H g ′ →g ψ g ′ ,n+1 (r, ω) = S2 σ g ′ →g,n+1 s (r, ω ′ → ω)ψ g ′ ,n+1 (r, ω ′ )dω ′ . • F g ψ g,n+1 (r, ω) = χ g,n+1 p (r) 4π G g ′ =1 1 -β g ′ ,n+1 (r) (νσ f ) g ′ ,n+1 (r)φ g ′ ,n+1 (r),
where φ g,n+1 (r) = S2 ψ g,n+1 (r, ω)dω ′ , ∀g ∈ {1, . . . , G}.

• Qg ψ g,n+1 (r, ω) = L ℓ=1 λ ℓ χ g,n+1 d,ℓ (r)C n+1 ℓ (r) The coefficients σ g,n+1 t (r), σ g ′ →g,n+1 s (r, ω ′ → ω), (νσ f ) g,n+1 , χ g,n+1 p (r) and χ g,n+1 d,ℓ (r) correspond to energy average values in [E g , E g-1 ] at time t = t n of the coeffi- cients σ t (t, r, E), σ s (t, r, ω ′ → ω, E ′ → E), σ f (t, r, E), ν(t, r, E), χ p (t, r, E) and
C n+1 ℓ (r) = 1 1 + λ ℓ δt n+1 C n ℓ (r) (3.2) + δt n+1 1 + λ ℓ δt n+1 G g ′ =1 β g ′ ,n+1 ℓ (r)(νσ f ) g ′ ,n+1 (r)φ g ′ ,n+1 (r).
The insertion of (3.2) into (3.1) finally yields the set of source problems:

(3.3)
Find over R × S 2 and ∀g the angular flux ψ g,n+1 (r, ω) such that: (L g -H g -F g ) ψ g,n+1 (r, ω) = S g,n (r, ω),

where:

• F g ψ g,n+1 (r, ω) = G g ′ =1 F g ′ ,g ψ g ′ ,n+1
and

F g ′ ,g ψ g ′ ,n+1 = χ g,n+1 p (r) 4π 1 -β g ′ ,n+1 (r) + L ℓ=1 λ ℓ β g ′ ℓ χ g d,ℓ δt n+1 1 + λ ℓ δt n+1 (νσ f ) g ′ ,n+1 (r)φ g ′ ,n+1 (r), • S g,n (r, ω) := ψ g,n (r, ω) V g δt n+1 + 1 1 + λ ℓ δt n+1
C n ℓ (r). Because of the coupling in the energy variable, system (3.3) is iterativelysolved over the index g with a numerical method that we will call "generalized Gauss-Seidel scheme" (these are the so called "outer iterations" in neutronics). The scheme goes as follows: let ψ g,n+1 (M) be the approximation of ψ g,n+1 at iteration number M . If we denote (3.4)

ψ n+1 =      ψ 1,n+1 ψ 2,n+1 . . . ψ G,n+1      ; φ n+1 =      φ 1,n+1 φ 2,n+1 . . . φ G,n+1      ; S n =      S 1,n S 2,n . . . S G,n      ,
then the scheme reads:

(3.5) M G,G ψ n+1 (M +1) = N G,G ψ n+1 (M ) + S n ψ n+1 (M =0) given,
where

A G,G = M G,G -N G,G , with (3.6) M G,G =      L 1 -H 1→1 0 • • • 0 -H 1→2 L 2 -H 2→2 • • • 0 . . . . . . . . . . . . -H 1→G -H 2→G • • • L G -H G→G      and (3.7) N G,G =      F 1,1 H 2→1 + F 2,1 • • • H G→1 + F G,1 F 1,2 F 2,2 • • • H G→2 + F G,2 . . . . . . . . . . . . F 1,G F 2,G • • • F G,G      .
System (3.5) is lower-triangular and is simply solved by forward substitution. Note that the difference between this scheme and a traditional Gauss-Seidel lies in the "explicit" treatment of the fission terms F g ′ ,g for g ′ ≤ g. This choice has been made for implementation purposes to minimize the number of computations: the fission terms are computed with the previous iterate so that their value remains fixed in the iterative resolution of (3.5) that is going to be introduced in (3.11).

It has been proven in chapter 1 of [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF] (theorem 1.6.1) that this scheme converges for small enough time steps. In MINARET, the iterations are performed until the average residual error in the scalar flux

(3.8) e n+1 outer (M + 1) :=      G g=1 φ g,n+1 (M+1) -φ g,n+1 (M) 2 L 2 (R) G g=1 φ g,n+1 (M+1) 2 L 2 (R)      1/2
goes below a given convergence threshold ε outer and ψ n+1 (∞) and φ n (∞) will denote the converged angular and scalar fluxes.

For g = 1, . . . , G and for a given iteration M , the problem to be solved in the forward substitutions to invert system (3.5) reads:

(L g -H g→g )ψ g,n+1 (M+1) (r, ω) = g ′ <g H g ′ →g ψ g ′ ,n+1 (M+1) + g ′ >g H g ′ →g ψ g ′ ,n+1 (
M) + G g ′ =1 F g ′ ,g ψ g ′ ,n+1 (M)
+ S g,n (r, ω), (3.9) which is a monoenergetic equation of the form

ω.∇ψ(r, ω) + σ t (r)ψ(r, ω) - S2 σ s (r, ω ′ → ω)ψ(r, ω ′ )dω ′ = q(r, ω), (3.10)
for all (r, ω) ∈ R × S 2 . The terms σ t , σ s , q must be understood as generic notations whose definition must be coherent with equation (3.9). Since equation (3.10) is integral in the angular variable and differential in space, a second numerical scheme is performed (called "inner or source iterations"). If ψ g,n+1 (M+1,m) is the approximation of ψ g,n+1

(M+1) at the m -th inner iteration, then ψ g,n+1 (M+1,m+1) is the solution of:

(3.11) L g ψ g,n+1 (M+1,m+1) (r, ω) = H g→g ψ g,n+1 (M+1,m) (r, ω) + S(r, ω), ∀g ∈ {1, . . . , G}, with S(r, ω) = g ′ <g H g ′ →g ψ g ′ ,n+1 (M+1) + g ′ >g H g ′ →g ψ g ′ ,n+1 (M)
+ G g ′ =1 F g ′ ,g ψ g ′ ,n+1 (M)
+ S g,n (r, ω).

It has been shown in [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF] (section 1.6.1.2 of chapter 1) that this strategy is equivalent to a Richardson scheme. The iterations are performed until the relative residual (3.12) e g,n+1 inner (m + 1) :=

φ g,n+1 (M,m+1) -φ g,n+1 (M,m) L 2 (R) φ g,n+1 (M,1) -φ g,n+1 (M,0) L 2 (R)
goes below a given convergence threshold ε inner .

The angular discretization of equation (3.11) has been performed with the discrete ordinates of order n technique (S n ), i.e., problem (3.11) is solved for a discrete number D of directions {ω d } D d=1 . The scattering terms H g ′ →g are computed in practice by assuming that the scatering cross sections are isotropic, i.e., that

σ s (r, ω ′ → ω) = σ s (r, ω ′ .ω),
This assumption allows to make an expansion of the terms H g ′ →g in Legendre polynomials (the indexes g, n, M and m have momentarily been removed to ease the reading). The expansion is truncated at a given order L so that

S2 σ s (r, ω ′ → ω)ψ(r, ω ′ )dω ′ ≈ L ℓ=0 2ℓ + 1 4π σ s,ℓ (r) S2 ψ(r, ω ′ )P ℓ (ω.ω ′ )dω ′ ,
where P ℓ is the Legendre polynomial of degree ℓ and σ s,ℓ = 2π

1 -1 σ s,ℓ (r, µ)P ℓ (µ)dµ. The use of the addition theorem P ℓ (ω ′ .ω) = 4π 2ℓ+1 ℓ m=-ℓ Y ℓ,m (ω ′ )Y ℓ,m (ω), where Y ℓ,m (ω) is the spherical harmonic of index (ℓ, m), eventually leads to S2 σ s (r, ω ′ → ω)ψ(r, ω ′ )dω ′ ≈ L ℓ=0 σ s,ℓ (r) ℓ m=-ℓ φ ℓ,m (r)Y ℓ,m (ω),
which is an expansion that depends only on the flux moments that are defined as

(3.13) φ ℓ,m (r) = S2 ψ(r, ω ′ )Y ℓ,m (ω ′ )dω ′ .
Note that φ ℓ,m (r) = φ(r) for ℓ = m = 0. More importantly, the fact that only the flux moments are required to compute the scattering term implies that the evaluation of the term N G,G ψ n+1 (M ) in system (3.5) does not need the knowledge of the angular fluxes ψ n+1 (M ) , but only the knowledge of flux moments φ n+1 ℓ,m (M ) (in the case where ℓ = m = 0, only φ n+1 (M ) is needed). We can therefore re-write system (3.5) as

(3.14) M G,G ψ n+1 (M +1) = N G,G φ n+1 ℓ,m,(M ) + S n φ n+1 (M =0) given,
where N G,G denotes the same operator as N G,G , but after having made the observation that the entries are flux moments.

The integrals in the angular variable of relation (3.13) are computed with a quadrature formula. In the case of MINARET, the level-symmetric quadrature rule has been used so that D = n(n + 2) in our case.

The treatment of the spatial variable is the point where the neutronic solvers most differ. In the case of MINARET, discontinuous Galerkin finite elements of arbitrary order have been implemented (however, the results presented in this paper involve only degree 1 polynomials -see remark 3.1). The three-dimensional spatial mesh is "partially unstructured" in the sense that it is built by extrusion of an initial twodimensional unstructured mesh. We refer to [START_REF] Moller | Minaret, a deterministic neutron transport solver for nuclear core calculations[END_REF] for further details on MINARET's mesh generator and also for a study made on the Jules Horowitz reactor, whose cylindrical geometry cannot be treated accurately on a cartesian grid and thus requires non cartesian grids. In addition to this, the finite element order can be locally adapted depending on the complexity of the phenomena involved in each region, which makes the solver be well-suited for the study of accidents.

Algorithm 1 summarizes the described two-stage nested iterative strategy implemented in MINARET.

Remark 3.1. The advantages of the hp adapativity of the MINARET solver will not be illustrated in the present work as our major goal is to provide tractable results about acceleration methods. We refer to [START_REF] Moller | Minaret, a deterministic neutron transport solver for nuclear core calculations[END_REF] as a first overview of these capabilities and to future papers in which MINARET will be used for accidental studies.

4. Definition of the numerical test cases and validation of the MINARET solver. We briefly explain in this section the two test cases that will be used to illustrate the numerical performances of the acceleration methods that are going to be discussed in the remaining of this paper.

The first test case (denoted below as "case A") represents a rod withdrawal and corresponds to a well-known benchmark in the literature of neutronics called TWIGL (see [START_REF] Langenbuch | Coarse-mesh flux expansion method for the analysis of space-time effects in large light water reactor cores[END_REF]). The geometry of the core is three-dimensional and the domain is 1b).

R = {(x, y, z) ∈ R 3 | 0 ≤ x ≤ 220 cm; 0 ≤ y ≤ 220 cm; 0 ≤ z ≤ 200 cm; }. A cross- sectional view at the height z = 180 cm is specified in table 1a
The results provided by the benchmark do not correspond to measured values, but to diffusion calculations, which is the reason why this test case was originally proposed for the validation of diffusion solvers. However, to the best of the authors knowledge, a specific kinetic benchmark for transport solvers does not exist in the literature.

Algorithm 1 The iterative strategy implemented in MINARET

1: for t n = ∆T to N ∆T do 2:

While(not converge) do (generalized GS iterations -see equation (3.5)) End While 16: end for It will therefore not be possible to study through this example the performances of MINARET in situations where diffusion is no longer accurate enough. Despite this fact, the benchmark will allow to

3:

Update fission operator

• do a simple validation of the solver's kinetic module by comparing our calculations with a diffusion solver (called MINOS; see, e.g., [START_REF] Baudron | MINOS: a simplified P N solver for core calculation[END_REF]) and with the original results presented in [START_REF] Langenbuch | Coarse-mesh flux expansion method for the analysis of space-time effects in large light water reactor cores[END_REF],

• provide tractable performance results. Regarding the validation, table 1b shows the power evolution computed with MINARET, MINOS and the original values from the diffusion benchmark given in [START_REF] Langenbuch | Coarse-mesh flux expansion method for the analysis of space-time effects in large light water reactor cores[END_REF]. As it can be noted, the MINARET solver presents the correct trend, which will be taken as an indication that the resolution has correctly been implemented.

The second test case (denoted below as "case B") uses the same geometry as the TWIGL benchmark but an oscillatory sequence of motion of the rods has been devised so that power fluctuations are produced. The simulated interval of time is [0, T ] with T = 250 s. (see table 2 for the details).

Both tests have been carried out with G = 2 energy groups, L = 6 precursors and vacuum boundary conditions. The degree of the discontinuous finite elements is P1 and the spatial mesh size is of order 20 cm. All the computations that will be presented hereafter have been obtained in a cluster of 38 nodes of 16 Gb memory, each one composed of 8 cores of 2814 MHz speed.

Remark 4.1. In the original TWIGL benchmark from the literature ( [START_REF] Langenbuch | Coarse-mesh flux expansion method for the analysis of space-time effects in large light water reactor cores[END_REF]), calculations are done in a quarter of a core with reflective boundary conditions in the inner parts of the core. In our case, the full geometry has been computed in order to be coherent with case B that has no spatial symmetries.

Remark 4.2. We insist on the fact that the test cases that have been considered will not illustrate the advantages of using transport rather than the classical diffusion approximation. They do not show the potentialities of MINARET's unstructured mesh either but we emphasize again that our choice has been determined with the aim of illustrating the algorithmic efficiency of our implementation and not the discretiza- tion per se; we thus provide results that allow to quantify the effect of the acceleration strategies. The other aspects about approximation rate and capabilities will be illustrated in subsequent works with the solver.

Sequential acceleration techniques.

The convergence of the iterative resolution of the multigroup problem given in (3.5) is often extremely slow and acceleration methods are required in order to obtain reasonable computing times.

We will recall in this section the two main traditional sequential accelerations. Their efficiency will be illustrated by some computations with MINARET. The first strategy is the Chebychev extrapolation in the outer iterations. It consists in adding a linear combination of the fluxes after each Gauss-Seidel iteration:

               M G,G ψ n+1 (M +1) = N G,G φ n+1 ℓ,m,(M ) + S n φ n+1 ℓ,m,(M +1) (r) = S2 ψ n+1 (M +1) (r, ω ′ )Y ℓ,m (ω ′ )dω ′ φ n+1 ℓ,m,(M +1) = α M+1 φ n+1 ℓ,m,(M +1) -φ n+1 ℓ,m,(M -1) + φ n+1 ℓ,m,(M -1) , M ≥ 1 φ n+1 ℓ,m,(M =0) given (0 ≤ ℓ ≤ L and -ℓ ≤ m ≤ ℓ).
(5.1)

We refer to [START_REF] Varga | Matrix iterative analysis[END_REF] for the exact form of the coefficients (α M ) and the theoretical foundations of this acceleration scheme.

The second acceleration scheme is the so-called diffusion synthetic acceleration (DSA) that has been added for the convergence of the inner iterations. It reads:

(5.2)

         L g ψ g,n+1 (M,m+1) (r, ω) = H g→g ψ g,n+1 (M,m) (r, ω) + S(r, ω), ∀g ∈ {1, . . . , G} φ g,n+1 (M,m+1) (r) = S2 ψ g,n+1 (M,m+1) (r, ω ′ )dω ′ ψ g,n+1 (M,m+1) (r, ω) = ψ g,n+1 (M,m+1) (r, ω) + e g,n+1 (M,m+1) (r).
The term e g,n+1 (M,m+1) is the solution of the diffusion problem

(5.3)                    -div   1 3 σ g t (r) + 1 V g ∆Tn+1 ∇e g,n+1 (M,m+1) (r)   + σ g t (r) + 1 V g ∆Tn+1 -σ g,g s (r) e g,n+1 (M,m+1) (r) = σ s (r) φ g,n+1 (M,m+1) (r) -φ g,n+1 (M,m) (r) , ∀r ∈ R e g (m+1) (r) = 0, ∀r ∈ ∂R,
where

φ g,n+1 (M,m) (r) = S2 ψ g,n+1
(M,m) (r, ω ′ )dω ′ . DSA is an acceleration scheme because it acts as a preconditioner of transport to solve equation (3.11). We refer to [START_REF] Adams | Fast iterative methods for discrete-ordinates particle transport calculations[END_REF] (sections I.D and II.B) and [START_REF] Moller | Éléments finis courbes et accélération pour le transport de neutrons[END_REF] for more theoretical details about this method. To preserve the total amount of neutrons, the spatial resolution of the DSA problem (5.3) is discretized with discontinuous Galerkin finite elements of the same degree as the one used for the transport problem (5.2). The discretized DSA problem is iteratively solved with a conjugate gradient method preconditioned by SSOR. If r i denotes the residual at the i-th iteration, the DSA iterations are performed until the ratio (5.4)

r i L 2 (R) r 0 L 2 (R)
goes below a given convergence threshold ε DSA .

We illustrate the performances of both acceleration methods through some numerical results obtained with MINARET for the test case A.

To begin with, table 3 lists the number of outer iterations M outer , inner iterations N inner and DSA iterations N DSA required to perform a propagation from time 0 to time 5/3 s. in an S 4 calculation (i.e, a calculation in which the S n discrete ordinates of order n technique is set to n = 4). We also provide the exact computing times obtained in our cluster. The convergence criteria associated to the errors (3.8), (3.12) and (5.4) have been fixed to ε outer = 10 -5 , ε inner = 10 -1 and ε DSA = 10 -2 . The product M outer N inner is also given as an estimation of the complexity of the resolution (the complexity added by the DSA can be neglected in a first approach in a sequential calculation). As the first case of table 3 shows, it is clear that the solver needs acceleration techniques in order to converge in a reasonable time. While the inclusion of the Chebyshev extrapolation (case 2) already represents a dramatic improvement in the computing time (by reducing about 10 times the number of outer iterations), this performance can still be improved by another factor of about 10 if the Chebyshev extrapolation is coupled with DSA in the inner iterations (case 3). This is achieved thanks to the reduction of the number of inner iterations. Another factor of about 3 can further be obtained if the initial guess φ n+1

(M =0,N =0)
of the outer iterations is well chosen (we will discuss here the case where L = 0, but note that it could very easily be generalized to any type of order L in the scatering expansion with Legendre polynomials). The classical choice is

(5.5) φ g,n+1 (M=0,N =0) = φ g,n (∞) , ∀g ∈ {1, . . . , G}.
This option is reasonable because for small time steps one can conjecture that the system does not change very much from t n to t n+1 . Other possibilities that exploit the information of the previous time steps have been explored (these are at the cost of storing additional information). One can first try a linear extrapolation of the flux:

(5.6)

φ g,n+1 (M=0,N =0) = φ g,n (∞) + t n+1 -t n t n -t n-1 φ g,n (∞) -φ g,n-1 (∞)
, ∀g ∈ {1, . . . , G}.

However, according to the point kinetics approximation, the behavior of the flux is rather exponential and another idea would be an exponential extrapolation:

(5.7)

φ g,n+1 (M=0,N =0) = φ g,n (∞) exp t n+1 -tn tn -t n-1 ln φ g,n (∞) φ g,n-1 (∞)
, ∀g ∈ {1, . . . , G}.

Another interesting option is to use the diffusion approximation to build a two-level propagation scheme. The idea goes as follows: the computation of the solution with the diffusion approximation can be obtained very quickly in comparison with the transport solution. For a given time t n+1 , we can therefore compute the solution coming from the diffusion (denoted here as φg,n+1 (∞) ) and use it as a predictor of the transport solution by setting:

(5.8) φ g,n+1 (M=0,N =0) = φg,n+1 (∞) , ∀g ∈ {1, . . . , G}..

As will be illustrated in the numerical results, this is a bad choice whose main problem is that the diffusion solution has a different trajectory than the transport one, hence the degraded computing times (the transport solver has to correct the trajectory and then converge to the transport solution). However, since the diffusion approximation seems to present the good trend, one can conjecture that

φ g,n+1 -φ g,n t n+1 -t n ≈ φg,n+1 -φg,n t n+1 -t n .
In this case, we can try:

(5.9)

φ g,n+1 (M=0,N =0) = φ g,n (∞) + φg,n+1 (∞) -φg,n (∞) , ∀g ∈ {1, . . . , G}.
The numerical results will show that this is a good starting guess. Furthermore, if we suppose that the trend is exponential as point kinetics suggests, an interesting initial guess could be:

(5.10)

φ g,n+1 (M=0,N =0) = φ g,n (∞) exp t n+1 -tn tn -t n-1 ln φg,n+1 (∞) φg,n (∞)
, ∀g ∈ {1, . . . , G}.

Starting guess Formula A (traditional)

φ g,n+1 (M=0,N =0) = φ g,n (∞)
B (linear extrapolation)

φ g,n+1 (M=0,N =0) = φ g,n (∞) + t n+1 -t n t n -t n-1 φ g,n (∞) -φ g,n-1 (∞) C (exponential extrapolation) φ g,n+1 (M=0,N =0) = φ g,n (∞) exp t n+1 -tn tn -t n-1 ln φ g,n (∞) φ g,n-1 (∞) D (plain multilevel) φ g,n+1 (M=0,N =0) = φg,n+1 (∞) E (multilevel linear) φ g,n+1 (M=0,N =0) = φ g,n (∞) + φg,n+1 (∞) -φg,n (∞) F (multilevel exponential) φ g,n+1 (M=0,N =0) = φ g,n (∞) exp t n+1 -tn tn -t n-1 ln φg,n+1 (∞) φg,n (∞)
Table 4: List of the explored starting guesses.

We summarize all the options in table 4. Their perfomances have been tested in the TWIGL benchmark case with a constant time-step of 5/3 s. In figure 1 we plot the computing times per time step as well as the cumulative ones. We also plot M outer and M outer N inner N DSA . From these figures, it seems thus clear that the use of a multilevel scheme outperforms the rest of the approaches provided that we do a linear or exponential extrapolation. The computing times are reduced by a factor of about 3 with this strategy. Options B and C provide a more moderate gain compared to the traditional starting guess A. As it can be observed from the figures, the speed up comes from the reduction of the number of outer iterations M outer , which results in a dramatic reduction of the total number of iterations M outer N inner N DSA . Note however that the starting guesses have been tried in a case in which diffusion is a good approximation to transport and it is therefore legitimate to wonder whether the starting guesses involving diffusion computations will still be efficient in other contexts. We believe that their performances should not be very degraded as long as the time step remains sufficently small, but computations in a more "transportoriented" test case are required to prove this claim.

6. Parallelization of the angular variable. We now turn to parallel acceleration techniques and start by recalling in this section the most classical and straightforward strategy employed in discrete ordinates codes:

The numerical scheme outlined in algorithm 1 puts in evidence that, for a given energy group g and a given inner iteration m, each angular flux ψ(r, ω d ) is computed independently (see lines 7 to 10 of algorithm 1) from the other unknowns ψ(r,

ω d ′ ) (d ′ = d).
The loop of lines 7 to 10 is therefore an embarrassingly parallel task that can be performed concurrently on several processors by uniformly distributing the set of angular fluxes to be treated among the different processors.

From an implementation point of view, the distribution of the tasks is performed in a master-slave fashion, that has the important advantage of alleviating the memory storage per processor in comparison with a sequential implementation. Indeed, each processor stores only the angular fluxes ψ(r, ω d ) of its assigned directions (and the moments of the flux are only stored by the master). Thanks to this fact, time-depend calculations involving a large number of directions and leading to HPC problems can be addressed. However, at the end of each inner iteration, the master gathers the angular fluxes {ψ(r, ω d )} D d=1 to compute the scalar flux φ(r) and then it performs the diffusion synthetic acceleration. As a consequence, as shown in, e.g., [START_REF] Jamelot | High performance 3D neutron transport on petascale and hybrid architectures within APOLLO3 code[END_REF] this part becomes a bottleneck for large numbers of processors because of the communication time and this isssue is one of the reasons that motivates to couple the parallelization of this variable with other strategies. Moreover, apart from this communication problem for large numbers of processors, the efficiency of this method is also a trade-off between:

• the number of directions assigned to each processor • the spatial complexity for the calculation of an unknown ψ(r, ω d ) (resolution of problem (3.11)) • the computing time required to perform the DSA step (that is run sequentially in MINARET) This trade-off is illustrated on table 5 and figure 2a where the numerical performances of this implementation in a strong scaling test are shown: the test case A has been performed for a fixed number of directions D = 24 with an increasing number of processors N that treat the loop over the directions. As it can be observed, for a reduced number of processors, the algorithm has excellent scalability properties (N ≤ 8). The behavior is degraded for larger values of N because the amount of work assigned to each processor decreases. The time to perform the loop on the angular directions is therefore reduced whereas the time to do the DSA remains constant because it is not parallelized: the DSA becomes a bottleneck. This issue could be overcome by its parallelization with spatial domain decomposition methods or multigrid techniques like in the works of [START_REF] Antonietti | Schwarz domain decomposition preconditioners for discontinuous galerkin approximations of elliptic problems: non-overlapping case[END_REF] and [START_REF] Napov | An algebraic multigrid method with guaranteed convergence rate[END_REF] respectively.

As a consequence of all these factors, in order not to lose much efficiency, there is a minimum number of directions ω d that need to be treated by each processor. In the present case, the most reasonable choice according to this criterion seems to be to assign N/D = 4 directions per processor (see table 5).

It is also desirable that the number of processors N is a divisor of the total number of directions D in order to have an uniform distribution of the tasks between processors. This is indeed a source of inefficiency as illustrated in table 5 for the case N = 10 (some processors will treat 3 directions and others only 2).

With the "optimal" number of N/D = 4 directions per processor being fixed, we perform a weak scaling test in which the angular S n approximation is increased (D increases) while incrementing the number N of processors. The results are summarized in table 6 and figure 2b where it can be noticed that, for a relatively reduced amount of processors, the efficiency is almost not degraded as N increases. As anticipated before, a degradation would be observed for large N because of significant communication times, but it is important to recall here that the S n angular accuracy should be chosen in coherence with the accuracy of the other of the variables (like, e.g., the spatial variable whose accuracy is given by the finite element polynomial approximation in our case). For this reason, and also for the time communication issues arising for large N , there is a limit in the total number of processors that is worth devoting to the parallelization of the angular variable. Thus, it is interesting to consider other parallelization strategies as a source of additional speed-up. As outlined in the introduction, significant efforts have been invested in previous works for the parallelization of the spatial variable with the KBA algorithm and promising results have also been obtained regarding the parallelization of the energy. As a complementary approach to all these methods, we propose to explore the additional speed-ups that the parallelization of the time variable could bring to the resolution of our equation. 7. Parallelization of the time variable. This section is organized as follows: after a brief recall of the basics of the parareal in time algorithm, an extension of the traditional theoretical speed-up formula will be proposed in order to properly take into account our particular case in which parareal is coupled with other iterative techniques at each time propagation. Finally, an analysis of the performances of the method for the resolution of transport transients with MINARET will be presented. The implemented results consider the parallelization of the time without coupling it with the parallelization of the rest of the variables. They are nevertheless representative enough of the accelerations that could be obtained in addition to the ones provided by the other methods.

D
7.1. The parareal in time algorithm. The unsteady problem (2.1) can be written in a more compact form:

(7.1) ∂y ∂t + A(t; y) = 0 , t ∈ [0, T ];
it is complemented with initial conditions: y(t = 0) = y 0 .

We assume that we have two propagators to solve (7.1): a fine one F τ1 τ0 (y(τ 0 )) that, starting from time τ 0 ∈ [0, T ] with the value y(τ 0 ), computes an approximation of the solution of (7.1) at time τ 1 ∈ [τ 0 , T ] accurately but slowly, and a coarse one G τ1 τ0 (y 0 ) that computes another approximation quickly but not so accurately (and not accurately enough). The fine propagator F can, e. g., perform the propagation of the phenomenon from τ 0 to τ 1 with small time steps δt with very accurate physics described by A. On the other hand, the coarse approximation G does not need to be as accurate as F and can be chosen much less expensive e.g. by the use of a scheme with a much larger time step ∆T ≫ δt or by treating "reduced physics" (i.e. by simplifying A into a less computer resources demanding operator).

In addition to these two propagators F and G, the parareal in time algorithm is based on the division of the full interval [0, T ] into N sub-intervals [0, T ] = N -1 n=0 [T n , T n+1 ] that will each be assigned to a processor P n , assuming that we have N processors at our disposal. The parareal in time algorithm applied to (7.1) is an iterative technique where, at each iteration k, the value y(T n ) is approximated by Y k n with an accuracy that tends to the one achieved by the fine solver when k increases. Y k n is obtained by the recurrence relation:

(7.2) Y k+1 n+1 = G Tn+1 Tn (Y k+1 n ) + F Tn+1 Tn (Y k n ) -G Tn+1 Tn (Y k n ), n = 0, ..., N -1 starting from Y 0 n+1 = G Tn+1 Tn (Y 0 n ). From formula (7.
2), it can first of all be seen by recursion that the method is exact after enough iterations. Indeed, for any n > 0, Y n n = F Tn 0 (y 0 ). However, convergence of Y k n to F Tn 0 (y 0 ) goes much faster than this as will be illustrated in our numerical example. Second, by the recurrence formula (7.2), the parareal in time algorithm can be cast into the category of predictor corrector algorithms, where the predictor

is G Tn+1 Tn (Y k+1 n ) while the corrector is F Tn+1 Tn (Y k n ) -G Tn+1
Tn (Y k n ) (we refer to [START_REF] Gander | Analysis of the parareal time-parallel time-integration method[END_REF] for a detailed discussion about the several possible interpretrations of the parareal method). 7.2. Algorithmics and theoretical speed-up. While the main results about the convergence properties of the method were studied in depth a decade ago (see, e.g. [START_REF] Lions | Résolution d'EDP par un schéma en temps pararéel[END_REF] [6] [START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF]), more recent efforts ( [START_REF] Minion | A hybrid parareal spectral deferred corrections method[END_REF] [5] [15] [START_REF] Berry | Event-based parareal: A data-flow based implementation of parareal[END_REF]) focus on the algorithmics to implement it in order to improve the speed-up provided by the original algorithm suggested in [START_REF] Lions | Résolution d'EDP par un schéma en temps pararéel[END_REF]. It consists of a master-slave type of implementation: the master carries out the coarse propagation in the whole time interval [0; T ], while each slave is in charge of fine propagations over its assigned time slice. Each slave sends F Tn+1 Tn (Y k n ) to the master so that the master computes the parareal corrections of equation (7.2), for all n. This original algorithm gives rise to two main computing drawbacks: the coarse propagation by the master is a bottleneck in the computation and the memory requirement in the master processor scales linearly with the number of slaves.

A remedy to both is a distributed algorithm that was suggested in [START_REF] Aubanel | Scheduling of tasks in the parareal algorithm[END_REF]: for each processor P n , the fine and the coarse solvers are propagated over [T n , T n+1 ] and the parareal correction Y k+1 n+1 is carried out. The process is repeated until convergence, i.e.

Y k+1 n -Y k n < η, for all n and where η is a given tolerance. A graphical description of the master-slave and distributed algorithms is shown in figures 3a and 3b in the ideal case where each processor is identical and the communication time is negligible. It is easy to realize that the distributed implementation does not change the number of iterations in order the parareal algorithm to converge but it provides better speed-ups than the original master-slave version (see formula (7.3) below). This is the reason why the distributed algorithm has been implemented in this study.

Remark 7.1. An improvement of the distributed algorithm from the processor utilization point of view is the so-called event-based parareal algorithm suggested by [START_REF] Berry | Event-based parareal: A data-flow based implementation of parareal[END_REF]. It exploits the fact that the coarse and fine propagations can be considered as a collection of tasks that can be treated by a processor as soon as their initial conditions are fulfilled. Once the task is performed, the processor treats the following task, if any, leading to an optimization of the processor utilization. However, since the present work focuses essentially on the feasibility and attainable speed-ups of parareal applied to equation (3.3), the distributed algorithm has been selected for its simpler implementation.

To the best of the authors knowledge, the theoretical analysis for the maximum attainable speed-up provided by the parareal algorithm in different types of algorithms has always been made under the assumptions that the computational cost of the fine (resp. coarse) solver does not vary from one processor to another. It is also assumed that the communication time is negligible. Under these two hypothesis, the maximum speed-up for the master-slave (S MS ) and distributed algorithms (S D ) are respectively (see [START_REF] Aubanel | Scheduling of tasks in the parareal algorithm[END_REF]):

(7.3) S MS = T seq T para,MS = N N r(1 + k * ) + k * ; S D = T seq T para,D = N N r + k * (1 + r) ,
where k * is the total number of parareal iterations needed in order to converge and

r = T G T F
, T G and T F are the computational costs of the coarse and fine propagators per processor. Note that S D > S MS for any k * , r and N > 1.

In the case that the fine and the coarse propagators solve each time step with an iterative numerical method (like in equation (2.1)), it is possible that the cost of the fine and the coarse solvers dramatically vary form one processor to another depending on the numerical complexity of the events that take place in each time slice ∆T (and this complexity cannot be predicted a priori). Figure 3c illustrates this fact. Formulae (7.3) need therefore to be extended to the broader case in which the computational costs T G = T G (k, p) and T F = T F (k, p) depend on the processor p and the parareal iteration k. It is easy to show that a more adequate formula for the speed-up in this case is (7.4)

                                 SD = T seq Tpara,D = T seq N -1 p=0 T G (0, p) + k * k=1 max p∈{0,...,N -1} (T G (k, p) + T F (k, p)) , SMS = T seq Tpara,MS = T seq N -1 p=0 T G (0, p) + k * k=1 N -1 p=0 T G (k, p) + max p∈{0,...,N -1} T F (k, p)
, where the communication time between processors has been neglected. Note that in the generalized formulae (7.4), we also find that SD > SMS since we have Tpara,MS -

Tpara,D ≥ k * k=1 p =p * (k)
T G (k, p) > 0, where T G (k, p * ) := max p∈{0,...,N -1}

T G (k, p). More importantly, there is an upper bound for these speed-ups given by N/k * so that

(7.5) SMS < SD < N k * .
Because of the upper bound N/k * , the algorithm turns out to be intrinsically less efficient than other parallelization techniques (like, e.g, spatial domain decomposition) in which one can theoretically expect to reduce by N the computing time with N processors. It could be stated that this lack of efficiency is due to the intrinsic difficulty that the parallelization of the time variable represents, because time is sequential by nature (see however section 7.4 for an current developement to improve this aspect). For this reason, the parareal in time algorithm should be employed to get additional speed-ups when other more efficient parallelization techniques reach saturation.

7.3. Numerical application. The parareal algorithm has been applied to the resolution of the test cases A and B. An S 4 transport propagator has been used as the fine solver whereas two coarse solvers have been tried out:

• an S 4 transport propagator (the only difference with the fine solver is the size of the time steps used: δt for F and ∆t = T n+1 -T n > δt for G), • a diffusion propagator. In all the computations, we will consider that convergence in the parareal scheme is achieved whenever the residual

r k = max n=0,...,N      G g=1 Φ g,n k -Φ g,n k-1 2 L 2 (R) G g=1 Φ g,n k 2 L 2 (R)      1/2 , k ≥ 1,
between the current parareal scalar flux solutions (φ g,n k ) 0≤n≤N 1≤g≤G and the previous ones

φ g,n k-1 0≤n≤N 1≤g≤G
go below a threshold η. Its value has been fixed to the precision of the numerical scheme that has been evaluated to be of order 10 -3 . Hence, η = 10 -3 . The tolerance in the convergence for the outer and inner iterations has been fixed to ε outer = 10 -5 and ε inner = 10 -1 . With this thresholds, parareal convergence has been achieved after only k * = 2, 3 or at most 4 iterations of the parareal in time algorithm.

In the following subsections, after giving a numerical proof of the convergence of the parareal algorithm in our case of study, some results about measured speed-ups will be presented.

7.3.1.

A numerical proof of the convergence. Figure 4 illustrates through one particular example that parareal effectively converges. The fine and coarse solvers are:

• Fine solver: S 4 transport with δt = 1/12 s.

• Coarse solver: Diffusion with ∆t = 60δt = 5 s. The points represent the errors

(7.6) e k (T n ) = G g=1 φ g,n k -φ g,n (∞) 2 L 2 (R) G g=1 φ g,n (∞) 2 L 2 (R) , ∀n ∈ {0, 1, ...N }, k ∈ {0, 1, 2}
between the parareal scalar fluxes Φ g,n k and the sequential fine solution φ g,n (∞) . 7.3.2. Speed-up performances. In the following strong and weak scaling tests, the fine solver has a fixed time step of δt = 1/12 s. Fig. 4: An example of the numerical convergence of the parareal algorithm in our neutron transport case.

Strong scaling results:. For the strong scaling analysis, the test case A has been solved with MINARET on an increasing number N of processors. The size of each sub-interval T n+1 -T n is constant for all n and equal to the time step of the coarse solver ∆t. In order to increase the number of processors solving the transient in the fixed time interval [0; 70 s.], the coarse time step has been reduced from ∆t = 60δt to ∆t = 20δt.

The measured speed-ups S measured are plotted in figure 5a. The evaluation of the maximum attainable speed-up formula SD can easily be derived thanks to the evaluation of the involved terms given in (7.4). When comparing the measured values with the theoretical ones, we obtain that S measured = SD so that the communication time between processors is negligible in our case and the obtained results are optimal.

Another interesting element to note is that one gets better speed-ups with a coarse diffusion propagator. This result seems reasonable because diffusion propagations are faster than transport ones.

We also observe that for a reduced number of processors, the speed-up increases until it reaches a plateau for more than 21 processors. This is due to the fact that, for large values of N , the size of the sub-intervals ∆t = T n+1 -T n decreases. As a result, the size of the problem addressed by each processor decreases and we reach a point in which the addition of more processors does not improve any longer the performances.

Weak scaling results. For this alternative evaluation of the scaling, we will focus on the test case B. We now consider the case in which the time step of the coarse solver ∆t is fixed to 60δt and the transient has a variable length T (N ) = N ∆t (i.e. the size of the problem linearly increases with the number N of processors). As an example, for N = 14, transient B will be solved in the time interval [0, 70 s.], whereas when N = 42 the time interval will be [0; 210 s.].

The measured speed-ups are plotted in figure 5b and, like in the strong scaling case, they are in perfect agreement with the theoretical formula SD . The most important result here is that the distributed algorithm can effectively speed-up long time calculations: the global trend for the speed-up is to increase with the number of processors. The discontinuity in the trend observed between a number of processors N = 41 and 42 comes from the fact that, due to the increasing size of the interval [0, T (N )], the number of parareal iterations k * raises from 3 to 4 at this stage.

7.4. Improvement of the efficiency through a degraded fine solver. In addition to the intrinsic lack of full efficiency of the plain parareal scheme (see relation (7.5) and the results obtained in section 7.3), the example of neutron transport also shows that, in general, the cost of the propagations of the fine solver (and also of the coarse one) can vary from one time t n to another time t n ′ and this imbalance leads to a loss in the performances of the parareal scheme (see figure 3c). In the present case, this is due to the fact that the resolution of each time step is performed by iterative techniques (our so-called outer and inner iterations). As a result, for a given convergence criterion in the outer and inner iterations, the cost of the fine propagations F Tn+1 Tn can dramatically vary from one processor to another depending on the numerical complexity that takes place in each time slice [T n ; T n+1 ] to which processor P n is assigned (notice, by the way, that this situation is very general and can occur in the resolution of other PDEs).

A strategy that is currently being studied is to degrade the fine solver F and build a solver F J that performs a reduced number J of outer iteration at each time step δt instead of making M n,εouter ≫ J outer iterations until convergence to ε outer of the residual error (see formula (3.8)). As a consequence, instead of converging "locally" at every parareal iteration, the convergence is expected to be achieved "globally" across the parareal iterations through the modified parareal scheme:

Y k+1 n+1 = G Tn+1 Tn (Y k+1 n ) + F Tn+1 J,Tn (Y k n ) -G Tn+1 Tn (Y k n ), n = 0, ..., N -1.
This development is in the same spirit as the works of Maday and Turinici in [START_REF] Maday | The Parareal in Time Iterative Solver: a Further Direction to Parallel Implementation[END_REF] or Minion in [START_REF] Minion | A hybrid parareal spectral deferred corrections method[END_REF], [START_REF] Emmet | Toward an efficient parallel in time method for partial differential equations[END_REF], where the idea of degrading the fine propagator and coupling it with other iterative methods such as spatial domain decomposition and spectral deferred corrections has already been tested. The main novelty of our proposed approach is that the degradation is purely algebraic and the main difficulty is of course to prove the convergence of the degraded parareal scheme. Preliminary results in this respect can be found in chapter 3 of [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF] in which the neutron diffusion equation has been adressed as an example. The study of this idea will be the topic of two subsequent papers: in the first one we will provide theoretical results of convergence and provide numerical results in a simple case. On a second stage, we plan to implement the scheme in the MINARET solver to carry a study about the real gain in efficiency that the method could provide.

7.5. A parareal in space and energy algorithm?. The parareal method was originally suggested for the time variable but it is quite straightforward to realize that the variable t of equation (7.1) is "dummy" in the sense that it could also represent a spatial variable: parareal provides also a method to parallelize 1D advection equations. The extension to 3D spatial advective problems like the current one (see equation (3.11)) seems therefore theoretically possible: for each angular unknown flux ψ(r, ω d ), the spatial mesh could be divided in a manner that is coherent with the direction of propagation ω d . Each part of the mesh could be assigned to a different processor that would perform the fine propagation (i.e. the transport propagation of ψ(r, ω d )). The coarse solver could consist in a diffusion approximation of the original equation (2.1). Furthermore, if we now observe the multigroup problem (equation (3.3)) or the outline of the resolution of a transient in algorithm 1, it can be seen that, for a given time step, the energy groups are solved through a loop that could in turn be also parallelized by the parareal in time algorithm: the coarse solver would propagate a reduced number of energy groups while the fine solver would propagate the problem for all the energy groups.

Conclusion. An analysis of traditional and innovative acceleration techniques that can be included in time-dependent neutron transport solvers without altering their main construction has been presented. A particular effort has been made to provide tractable numerical results about the performances of the methods. This has been achieved by the presentation computational results derived in a classical benchmark test case from the literature and obtained with the newly developed timedepend module of the MINARET solver.

Hence, it has been shown that the inclusion of the classical Chebyshev extrapolation and diffusion synthetic acceleration reduce the computing times by a factor of about 100. It has further been noted that one can still reduce the computing time by a factor of about 3 with the use of an appropriate starting guess for the outer iterations that involves diffusion propagations.

Regarding parallel accelerations, it has first been explained how the parallelization of the angular directions can efficiently speed-up calculations under several conditions (the efficiency of an S 12 calculation involving 168 directions is 90%). However, the parallelization of other variables is justified because there are limitations in the amount of processors that one should devote to this strategy for discretization and communication reasons. For this reason, it is interesting to explore complementary strategies such as the parallelization of the time variable by the parareal in time algorithm. The efficiency of this method is much lower than the performances provided by other parallelization techniques, but this is due to the difficult task of parallelizing a variable that is sequential by nature. It has nevertheless been illustrated that the method can provide additional speed-ups for the computation of -long time-neutron transport transients (e.g., a transient of 210 s. can be accelerated by a factor of 6 with 50 processors). In an ongoing work, a parareal scheme that involves a degraded fine solver with truncated outer iterations is being analyzed in an attempt to improve the efficiency of the method (see chapter 3 of [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF] for preliminary results).

Acknowledgments. This work was supported in part by the joint research program MANON between CEA-Saclay and University Pierre et Marie Curie-Paris 6.

  χ d,ℓ (t, r, E). We also have β g,n+1 (r) = L ℓ=1 β g,n+1 ℓ (r). The Euler backward scheme applied to the precursors' equation provides C n+1 ℓ (r) for any ℓ ∈ {1, . . . , L}:

  . The first group of rods (blue) is withdrawn from t = 0 (z = 100 cm measured starting from below) until t = 26.6 s. (z = 180 cm) at a constant speed. The second group of rods (red) is inserted from t = 7.5 s. (z = 180 cm) until t = 47.7 s. (z = 60 cm) and the simulated interval of time is [0, T ] with T = 70 s. (see table

4 : 5 : 6 : 8 : 9 :

 45689 for g = 1 to G do Update scattering (except self-scattering) While(not converge) do (source iterations) 7: for ω = ω 1 to ω D do Update self-scattering Solve spatial problem (3.11) for ω

  Cross-sectional view of the core (z = 180 cm).

  ) TWIGL (diffusion) (b) Rods' vertical motion and power evolution.

  Cross-sectional view of the core (z = 180 cm).

  Rods' vertical motion and power evolution.

Fig. 1 :

 1 Fig. 1: Performances of the initial guesses.

Fig. 2 :

 2 Fig. 2: Efficiency in the parallelization of the angular variable (case A)

Fig. 3 :

 3 Fig. 3: Two different algorithms to implement the parareal in time method (a-b) and an illustration of the imbalance in the tasks (c) observed when the parareal algorithm is coupled with other iterative schemes for each time step propagation (in the example, k * = 3 and N = 7 processors and the red/green colors represent the coarse/fine propagations).

  Strong scaling results (test case A).

Fig. 5 :

 5 Fig. 5: Scaling results with the pararel in time algorithm.

Table 1 :

 1 Case A (TWIGL benchmark).

Table 2 :

 2 Case B.

Table 3 :

 3 Case Chebyshev DSA M outer N inner N DSA M outer N inner Computing An illustration of the impact on the speed-up performances of the Chebyshev extrapolation and the DSA.

	time (s)

Table 5 :

 5 Efficiency in the strong scaling test for the angular variable (case A)

		24 24	24	24	24	24	24	24
	N	1	2	4	6	8	10	12	24
	D/N proc	24 12	6	4	3	2 or 3	2	1
	Efficiency 1	1 0.99 0.96 0.933	0.77	0.87 0.76
	S n approx 2	4	6	8	10	12
	N dir		8	24	48	80	120	168
	N proc	2	6	12	20	30	42
	N dir /proc 4	4	4	4	4	4
	Efficiency	1 0.96 0.95 0.90 0.922 0.90

Table 6 :

 6 Efficiency in the weak scaling test for the angular variable (case A).

The authors are also indebted to AREVA-NP and EDF for their financial support. The authors would like to thank A.M. Baudron for fruitful discussions.


John R Tucker

Alireza Mowla

Jeremy Herbert

Miguel A Fuentes

Craig S Freakley

Karl Bertling

Yah Leng Lim

Ranveer S Matharu

Julien Perchoux

Thomas Taimre

Stephen J Wilson

Aleksandar D Rakić
email: rakic@itee.uq.edu.au

Self-mixing sensing system based on uncooled vertical-cavity surface-emitting laser array: linking multichannel operation and enhanced performance

Many systems utilise a monolithic array of semiconductor lasers to increase the capability and functionality of single channel systems [1][2][3][4]. Such systems are able to dramatically increase the amount of data that can be acquired within a certain time frame or the amount of time it takes to obtain a certain amount of data. This is particularly important in the context of biomedical and industrial applications that have unstable or dynamically changing flow velocities.

In recent times, there has also been increasing interest in acquiring images with a monolithic array of lasers using the self-mixing (SM) effect [3,4]. This sensing technique is based on the interference of light reflected from an external target with the light inside the laser cavity, which leads to measurable changes in the output power and the laser junction voltage. This is an attractive setup as the laser is used as both the source and detector, which greatly reduces the complexity of the optical alignment and makes the system more compact.

Initially, it may seem that simultaneous multichannel operation of SM sensors may lead to problems common in other systems based on a monolithic array of lasers such as optical, electrical and thermal crosstalk (XT). A fraction of light emitted from a vertical-cavity surfaceemitting laser (VCSEL) in a dense monolithic array will, upon reflection from the external target, inevitably fall on the neighbouring elements in the array [5]. In a SM system based on such geometries optical XT may occur through the injection locking between different channels (VCSELs) of the array [6,7]. However, this is extremely unlikely to occur in most applications as injection locking requires a number of conditions to be met including mutual coherence and mode matching. This may be further avoided by a small change in the operating temperature or injection current of the laser [8]. Electrical XT may also arise through electrical connections to the laser array but other systems have been developed where this form of XT is insignificant at frequencies below 1 GHz [1,9,10]. Electrical XT may also be caused by the sharing of carriers via resistive paths in the laser array itself although this can be circumvented with appropriate semiconductor isolation structures [11].

Based on the results from other laser array based systems [2], it may also be expected that the heat generated by multiple lasers may decrease the performance of individual SM sensors. Surprisingly, we find that the signalto-noise ratio (SNR) of a monolithically integrated SM VCSEL sensor arrays is actually enhanced rather than reduced with multichannel operation. This is a consequence of the following phenomenon: maximum SNR achievable (over a range of operating injection currents) is increasing with the laser temperature, a result akin to that observed by Matharu et al. in a single channel system [12]. The remainder of this letter is therefore focused on describing an experiment that proves this phenomenon and discussing how the performance may be further enhanced through the array design.

The effect of multichannel operation in SM sensors was investigated in an uncooled 24×1 VCSEL array custom manufactured by Litrax Technology Co., Ltd with the system shown in Fig. 1. The VCSEL array was mounted in a surface mount ceramic package and individual VCSELs were wire-bonded to pins of the package. The VCSELs are based on proton implanted technology and have a pitch of 220 µm. A FLIR i7 thermal imaging camera was used to examine the temperature profile of the array for different operational settings.

A single macro lens (f =105 mm) was employed to focus all 24 laser beams onto a rotating disk. The target was placed 35 cm from the laser array, which resulted in a beam spacing of 220 µm and a magnification factor of 1. The laser array was oriented vertically. The aluminium disk, 10 cm in diameter, was sandblasted to provide a rotating diffusive target. The disk was driven by a DC servo motor with a 43:1 gear reduction controlled by an Elmo Whistle motor controller. The stability of the disk's angular velocity for the duration of the experiments was better than 1%, as established by monitoring the frequency of the motor's rotary encoder. A custom built 24-channel laser driver-receiver was used to individually bias each laser with a constant injection current. The SM signals were obtained through terminal voltage variations across individual VCSELs. The terminal voltage fluctuations were then individually amplified using 24 ac-coupled, single stage, low noise preamplifiers with gain G = 100. Additional stages of amplification with variable gain (G = 64 to 6400) were used to bring the signal to a level suitable for processing by a 16-bit data acquisition card.

Each channel had a sampling rate of 62.5 kHz. The power spectral density of the time domain signal was calculated in the LabVIEW programming environment to obtain the SM signal spectrum and 32 consecutive spectra were averaged. The Doppler frequency was extracted by fitting the power spectral density spectrum to a Gaussian distribution of Doppler components [13].

In this work, we are interested in uncooled operation of the laser array as Peltier devices consume high amounts of power and make a sensing system less compact and portable. For uncooled operation, the temperature of a single laser device may be increased by increasing the injection current or the amount of power dissipated as heat [14]. In array based operation, the heat generated by a single laser spreads to other lasers in the array. The resulting temperature rise from a neighbouring laser is inversely proportional to the pitch [15], which leads to the central lasers showing the highest rise in temperature when all lasers are operating simultaneously. To examine the effect of multichannel operation on the performance of a SM sensor we first measured the SNR of one of the central lasers (Channel 12) for a wide range of bias currents from threshold to 10 mA. We then remeasured the SNR for the same bias currents while turning on other lasers in the array. To maximise the temperature rise we simultaneously turned on the remaining 23 lasers in the array. These lasers, which we will subsequently referred to as the neighbouring lasers, were biased at the same injection current while the bias current of the VCSEL under investigation was changed. Figure 2 shows the SNR of the central (sensing) VCSEL over a range of bias currents.

Each curve corresponds to a different operating current in all neighbouring VCSELs. The well-known increase and subsequent decrease of SNR with bias current can be observed for all operating conditions [16,17]. For single channel operation (neighbouring VCSELs not operating), the SNR increases rapidly above threshold before reaching a maximum. The SNR then gradually decreases as the injection current is increased beyond the current for the maximum SNR. When the currents in the neighbouring lasers were increased to 2 mA, which is below the threshold current of each laser, the SNR in the sensing VCSEL increased significantly. This indicates that the increase in SNR is not due to optical XT as the lasers have no observable signal level below the threshold current. Subsequent measurements on a probing station showed that the electrical isolation between adjacent elements was of the order of 5 MΩ, which suggests that the increase in SNR was due solely to the temperature rise in the sensing laser. The injection current for the optimum SNR also slightly increases when the neighbouring lasers are turned on, which suggests that the initial temperature of the laser array is on the hot side of the parabolic threshold current temperature profile [12,18]. Further increase in operating current of the neighbouring VCSELs leads to a monotonic increase in maximum achievable SNR.

To investigate this further, we measured the SNR of the central laser (Channel 12) while simultaneously operating six other lasers on the array. In each experiment, all seven lasers were equally spaced. Three different sets of equal spacing were investigated. The injection current in all of the neighbouring lasers was set to 10 mA. Therefore, in each experiment, the same electrical power was delivered to the laser chip. To examine this issue in more detail, we used a thermal IR camera to obtain the temperature profile of the array. Figure 3 shows a close-up photo of the array as well as the thermal images for the different separation distances used in the experiment.

Figure 4 shows the SNR of the central VCSEL as a function of injection current. Three curves correspond to three different separation configurations, with corresponding temperature distributions shown in Fig. 3. Clearly, the SNR increases as the separation distance between lasers decreases, which can be attributed to the higher temperature of the central (sensing) VCSEL due to the heating sources being concentrated nearer the sensing VCSEL.

A consequence of this result is that the SNR may be further improved by manufacturing an array with elements spaced more closely together. This would not only increase the spatial resolution of an imaging system based on this array, but also allow for a smaller lens to be used due to the smaller array size, making the whole optical system more compact.The temperature rise in the array may also be increased by increasing the element diameter as suggested by Osinski et al [15]. However, larger diameter VCSELs may exhibit multiple transverse modes, which can affect the accuracy and performance of SM sensors [19].

While a change in operating current and temperature strongly affects the SNR of the system it also the emission wavelength of the We show here that this substantial change in lasing wavelength has predictable, and almost negligible effect, on the Doppler frequency and calculated target velocity. The dependence of the lasing wavelength on current and heatsink temperature was measured by placing the laser chip in a temperature-controlled mount. Figure 5 (1) where T is the heatsink temperature in • C, and I is the laser current in mA, and λ is wavelength in nm. The change in lasing wavelength with current is mainly caused by heating of the laser cavity [20].

At each temperature-current pair (represented by black dots in Fig. 5) a SM Doppler spectrum was acquired and the Doppler frequency was extracted. The total shift in Doppler frequency caused by changes in current and temperature was 16 Hz (from 2.784 kHz to 2.768 kHz) which corresponds to a fractional change of approximately 0.6%. This systematic error may be removed by calculating the wavelength at any given point in the current-temperature space using (1). Figure 3 shows that temperature variation under different operating regimes of our array is much smaller than the current-temperature range used in this experiment. Therefore, we can safely conclude that the effect of current and temperature change in an uncooled VCSEL array will not only be predictable and easily removable but shows negligible effect on the performance of such arrays. This means that there is no significant reduction in the accuracy of the velocity determination with multichannel operation of a SM array based sensor, which allows it to be used for a wide range of sensing applications.

In summary, we have compared the performance of single channel and multichannel operation of an uncooled single mode monolithic 24×1 VCSEL array in a SM sensing system. We found that the SNR of monolithically integrated SM sensor arrays is inherently enhanced through multichannel operation, as a result of the increased operational temperature of the sensor. This effect can be augmented by reducing the spacing between VCSELs on the array or increasing the diameter of each laser. Finally, the temperature rise due to multichannel operation had no significant effect on the measured Doppler frequency, which corresponds to no significant change in the accuracy of the measured velocity. These findings suggest that SM imaging systems based on VC-SEL arrays have the potential to produce images with high spatial resolutions at high frame rates. These features make parallel SM sensing systems especially attractive for biomedical and industrial applications where there are unstable or dynamically changing flow velocities.

This research was supported under the Australian Research Councils Discovery Projects funding scheme (DP 120 103703). Y.L.L. acknowledges support under the Queensland Governments Smart Futures Fellowships programme.

Fig. 1 .

 1 Fig. 1. Experimental setup for measuring velocity on a rotating disc. The disc is tilted around the vertical axis by 5 • to produce a small velocity component in the direction of the laser beam.

Fig. 2 .

 2 Fig. 2. SNR of central laser (Channel 12) versus injection current with different levels of current in neighbouring lasers. The current in the neighbouring lasers are all biased at the same current while the current of the central laser is varied.( circles -Single Channel, squares -Multichannel 2 mA, up triangles -Multichannel 5 mA, down triangles -Multichannel 8 mA, diamonds -Multichannel 10 mA).

Fig. 3 .

 3 Fig. 3. Close up photo of the array and thermal images for the different separation distances used for the experiments in Fig. 4. a) Photo, b) 660 µm spacing, c) 440 µm spacing, d) 660 µm spacing. Channel 1 is located at the left hand side of the images.

Fig. 4 .

 4 Fig. 4. SNR of central laser (Channel 12) versus injection current with simultaneous operation of six neighbouring lasers biased at 10 mA. There is an equal separation distance between lasers with the separation distance varying for different experiments. (circles -660 µm spacing, squares-440 µm spacing, triangles -220 µm spacing).

  shows the VC-SEL wavelength dependence on current and temperature measured on a stand-alone VCSEL of the same design as the VCSEL array used in this work. Measurements were executed for a range of temperatures between 10 and 60 • C and currents up to 10 mA. Polynomial fit to these experiments can be conveniently represented by the following function: λ(T, I) = 844.1+0.063T +0.066I +0.00029T I +0.0128I 2 ,

Fig. 5 .

 5 Fig. 5. Surface fit of the variation in wavelength with current and temperature to experimental data from stand-alone Litrax VCSEL.