Bertrand

Christian, Emile, Eléonore Bruno Bruno

François, Frank, Gérald Fannie Flora

Jean-Philippe Laurent

Mark, Olivier Marie-Agnès Lucie

Romain Séthy

Sofiane

Thierry Thibault

Didier, Eric David

Massively Parallel Cartesian Discrete Ordinates Method for Neutron Transport

Keywords: Parallelism, distributed computing, HPC, multi-core, vectorization, taskbased runtime system, Cartesian S N, Sweep

High-fidelity nuclear reactor core simulations require a precise knowledge of the neutron flux inside the reactor core. This flux is modeled by the linear Boltzmann equation also called neutron transport equation. In this thesis, we focus on solving this equation using the discrete ordinates method (S N ) on Cartesian mesh. This method involves a source iteration scheme including a sweep over the spatial mesh and gathering the vast majority of computations in the S N method. Due to the large amount of computations performed in the resolution of the Boltzmann equation, numerous research works were focused on the optimization of the time to solution by developing parallel algorithms for solving the transport equation. However, these algorithms were designed by considering a super-computer as a collection of independent cores, and therefore do not explicitly take into account the memory hierarchy and multi-level parallelism available inside modern super-computers. Therefore, we first proposed a strategy for designing an efficient parallel implementation of the sweep operation on modern architectures by combining the use of the SIMD paradigm thanks to C++ generic programming techniques and an emerging task-based runtime system: PaRSEC. We demonstrated the need for such an approach using theoretical performance models predicting optimal partitionings. Then we studied the challenge of converging the source iterations scheme in highly diffusive media such as the PWR cores. We have implemented and studied the convergence of a new acceleration scheme (PDSA) that naturally suits our Hybrid parallel implementation. The combination of all these techniques have enabled us to develop a massively parallel version of the S N Domino solver. It is capable of tackling the challenges posed by the neutron transport simulations and compares favorably with state-of-the-art solvers such as Denovo. The performance of the PaRSEC implementation of the sweep operation reaches 6.1 Tflop/s on 768 cores corresponding to 33.9% of the theoretical peak performance of this set of computational resources. For a typical 26-group PWR calculations involving 1.02 × 10 12 DoFs, the time to solution required by the Domino solver is 46 min using 1536 cores.

iii S N cartésien massivement parallèle pour la simulation neutronique Résumé : La simulation haute-fidélité des coeurs de réacteurs nucléaires nécessite une évaluation précise du flux neutronique dans le coeur du réacteur. Ce flux est modélisé par l'équation de Boltzmann ou équation du transport neutronique. Dans cette thèse, on s'intéresse à la résolution de cette équation par la méthode des ordonnées discrètes (S N ) sur des géométries cartésiennes. Cette méthode fait intervenir un schéma d'itérations à source, incluant un algorithme de balayage sur le domaine spatial qui regroupe l'essentiel des calculs effectués. Compte tenu du très grand volume de calcul requis par la résolution de l'équation de Boltzmann, de nombreux travaux antérieurs ont été consacrés à l'utilisation du calcul parallèle pour la résolution de cette équation. Jusqu'ici, ces algorithmes de résolution parallèles de l'équation du transport neutronique ont été conçus en considérant la machine cible comme une collection de processeurs mono-coeurs indépendants, et ne tirent donc pas explicitement profit de la hiérarchie mémoire et du parallélisme multi-niveaux présents sur les super-calculateurs modernes. Ainsi, la première contribution de cette thèse concerne l'étude et la mise en oeuvre de l'algorithme de balayage sur les super-calculateurs massivement parallèles modernes. Notre approche combine à la fois la vectorisation par des techniques de la programmation générique en C++, et la programmation hybride par l'utilisation d'un support d'exécution à base de tâches: PaRSEC. Nous avons démontré l'intérêt de cette approche grâce à des modèles de performances théoriques, permettant également de prédire le partitionnement optimal. Par ailleurs, dans le cas de la simulation des milieux très diffusifs tels que le coeur d'un REP, la convergence du schéma d'itérations à source est très lente. Afin d'accélérer sa convergence, nous avons implémenté un nouvel algorithme (PDSA), adapté à notre implémentation hybride. La combinaison de ces techniques nous a permis de concevoir une version massivement parallèle du solveur S N Domino. Les performances de la partie Sweep du solveur atteignent 33.9% de la performance crête théorique d'un super-calculateur à 768 cores. De plus, un calcul critique d'un réacteur de type REP 900MW à 26 groupes d'énergie mettant en jeu 10 12 DDLs a été résolu en 46 minutes sur 1536 coeurs.

Mots clés : Parallélisme, calcul distribué, HPC, multi-coeur, vectorisation, ordonnanceur à base de tâche, S N cartésien, Sweep.

Discipline : Informatique

LaBRI (UMR CNRS 5800) Université de Bordeaux, 351, cours de la libération 33405 Talence Cedex, FRANCE Equipe Projet Inria HiePACS1 Inria Bordeaux -Sud-Ouest-200818243Z, 200, avenue de la vielle tour, 33405 Talence Cedex, FRANCE

Introduction

In today's world, the energy demand is constantly growing. Among the various sources of energy, the nuclear represents 11% of the overall total energy production 1 . Nuclear energy (actually electricity) is produced by the mean of nuclear power plants, using fissile fuel placed in reactor cores. Since the beginning of the nuclear industry in the early 1950s, nuclear reactor core operators, such as Electricité de France (EDF), conduct engineering studies to enhance safety and efficiency of the nuclear power plants. Indeed, before operating the reactor, safety reports require to precisely evaluate the location and the magnitude of the maximal power peak within the reactor core (pin peak power). This is mandatory in order to fulfill the regulatory hurdles. Efficiency concerns will require the evaluation of the reactor longest cycle length as possible, while maximizing the nominal power (see [START_REF] Sanchez | Prospects in deterministic three-dimensional whole-core transport calculations[END_REF] for more details). In order to meet these goals, nuclear engineering studies require to achieve high-fidelity predictive nuclear reactor core simulations. These simulations involve coupled multi-physics calculations that encompass thermal-hydraulics and neutronic studies. In particular, the neutronic studies, or neutron transport calculations, on which we focus in this thesis, consist in describing precisely the neutron flux distribution inside the reactor core. This flux, which corresponds to the neutron phase-space density, depends on seven variables: three in space ( #" r = (x, y, z) t ), one in energy (E), two in direction ( #" Ω ≡ (θ, φ)) and one in time (t). Thus, the precise simulation of the neutron flux distribution in the reactor core would require a tremendous computational effort.

Indeed, let us consider the core of a Pressurized Water Reactor (PWR) 900 MW. From the presentation in [START_REF] Reuss | Précis de neutronique[END_REF], we make the following observations.

• There are 157 fuel assemblies in the core, each of which being formed of 289 pin-cells.

Each pin-cell is composed of O [START_REF] Barrault | Efficient Parallel Resolution of the Simplified Transport Equations in Mixed-Dual Formulation[END_REF] radial zones and of O(50) axial zones.

• The large variations of the energetic spectrum of neutrons requires considering O(2 • 10 4 ) spatial mesh points.

• We must consider O( 102 ) angular directions for describing the traveling directions of the neutrons.

It follows that O (10 12 ) values of the flux should be evaluated for each timestep, which represents a very large amount of data. Thereby, to cope with daily industrial calculations, a two-step homogenization approach [START_REF] Sanchez | Prospects in deterministic three-dimensional whole-core transport calculations[END_REF] is generally employed. The first step of this approach leads to the evaluation of homogenized cross-sections using a lattice code, while the second step consists in evaluating the neutron flux in the whole reactor core using a core solver. The core solver is generally based on a simplified model (e.g. diffusion) of the exact neutron transport problem, and thus in fact comprises modelization errors. In order to quantify these modelization errors, existing in core industrial solvers, it is necessary to solve the Boltzmann Transport Equation (BTE) with a reference solver. Two main families of methods are used to perform reference calculations: probabilistic and deterministic methods. The first ones consist in using statistical tools to simulate the history of a large number of neutrons in the core, taking into account all the interactions that may occur between neutrons and the matter in the core. Hence, probabilistic methods enable to avoid the phase-space discretization problems, and can handle complex geometries, which make them very attractive for reactor physics analysis [START_REF] Bielajew | PRESTA: the Parameter Reduced Electron-Step Transport Algorithm for Electron Monte Carlo Transport[END_REF][START_REF] Briesmeister | MCNP TM -A General Monte Carlo N-Particle Transport Code[END_REF][START_REF] Salvat | PENELOPE-2006: A code System for Monte Carlo Simulation of Electron and Photon Transport[END_REF]. Unfortunately, the convergence of probabilistic methods allows to estimate the phase-space density with an accuracy that converges rather slowly. The other class of methods, reference deterministic methods, are computationally very demanding because they require a full discretization of the BTE to achieve acceptable levels of accuracy. For this reason, until a few years ago these methods were impracticable for 3D cases because of the limitations on the computing power.

Computer resources have since grown in capability, and several research works are therefore being conducted in developing reference 3D deterministic solvers [START_REF] Gregory G Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF][START_REF] Kaushik | Enabling High-Fidelity Neutron Transport Simulations on Petascale Architectures[END_REF][START_REF] Sanchez | Prospects in deterministic three-dimensional whole-core transport calculations[END_REF][START_REF] James S Warsa | Krylov subspace iterations for deterministic k-eigenvalue calculations[END_REF]. These are based on either discrete ordinates (S N ), or spherical harmonic (P N ), or Method of Characteristics (MOC) which is a special case of S N . All these methods share the same energetic discretization, and differ on angular and spatial discretizations. The S N method consists in considering only a finite set of angular flux components ψ( #" r , E, #" Ω i ) that corresponds to a finite set of carefully selected angles #" Ω i . In the P N method, the transport equation is projected onto a set of spherical harmonics, allowing to mitigate the fundamental shortcoming present in the S N method: the "ray-effects". The MOC uses the same energetic and angular discretizations as the S N method, but its spatial discretization is based on the long characteristics. This discretization enables the MOC to deal with the heterogeneity and complexity of the reactor core like probabilistic methods. Although the MOC has been proven to be very efficient in 2D, large-scale 3D cases are still very computationally demanding [START_REF] Sanchez | Prospects in deterministic three-dimensional whole-core transport calculations[END_REF].

In this thesis, we consider the problem of solving accurately and efficiently the steady-state BTE using the S N method on Cartesian mesh. As previously mentioned, this reference method enables the validation of approximate core industrial solvers. In particular, during the reactor core refueling process, it must be guaranteed that the optimization of a new fuel-loading pattern (see [START_REF] Yamamoto | A Quantitative Comparison of Loading Pattern Optimization Methods for in-core Fuel Management of PWR[END_REF]) can be done in a short period of time [START_REF] Sanchez | Prospects in deterministic three-dimensional whole-core transport calculations[END_REF]. Consequently, a reference calculation, involving several coupled multi-physics iterations, has to be completed in a limited amount of time. It is therefore necessary to manage to develop highly-optimized algorithms and numerical methods to tackle the large amount of calculations required by the S N method.

In front of this problem, the landscape of today's parallel computers, on which the simulations have to be executed, has dramatically shifted since the end of frequency scaling of monolithic processors, which has motivated the advent of the multicore processor in the early 2000s [START_REF] Fruehe | Multicore Processor Technology[END_REF]. Modern clusters are composed of heterogeneous computing nodes. These computing nodes are equipped with processors having tens of CPU cores, capable of issuing vector instructions on wide registers. A key point to note here is that the off-chip bandwidth, determining the latency of memory accesses between the computer main memory and the CPU registers, is not growing as fast as the computing power of these CPUs [START_REF] Mckee | Reflections on the Memory Wall[END_REF], leading to a "memory wall". The consequence of this is to put a dramatic emphasis on the sustainable peak-performance per core, since the gap between this metric and the theoretical peak performance trends to increase for memory bound applications. Furthermore, computing nodes may comprise accelerators like Graphics Processing Units (GPUs) and manycore devices. This represents a huge amount of computing power that can be used to tackle large numerical simulations in science and engineering such as predictive nuclear reactor core simulations. As a consequence of this hardware evolution, several projects are dedicated in the scientific community to the improvement of the performance of numerical simulation codes. Hence, one solution for maximizing the sustainable peak performance of a computational kernel is to rely on highly optimized external numerical libraries. For instance, the Trilinos project [START_REF] Michael | An Overview of the Trilinos Project[END_REF] aims at providing optimized packages, developed with advanced objectoriented techniques and state-of-the-art parallel programming paradigms, which can therefore be directly used as building blocks for solving large-scale and complex multi-physics problems. For some specialized applications, as the neutron transport simulation problem, the solution for achieving high-performance on emerging architectures is oriented towards computer code modernization as justified by the large number of recent researches on this subject [START_REF] Paul | Sustained Systems Performance Monitoring at the US Department of Defense High Performance Computing Modernization Program[END_REF][START_REF] Heroux | ASC Co-design Proxy App Strategy[END_REF][START_REF] Jundt | Optimizing codes on the Xeon Phi: a case-study with LAMMPS[END_REF][START_REF] Mielikainen | Intel Xeon Phi accelerated Weather Research and Forecasting (WRF) Goddard microphysics scheme[END_REF][START_REF] Möller | A Case Study on Using a Proto-Application as a Proxy for Code Modernization[END_REF][START_REF] Radhakrishnan | Using Coarrays to Parallelize Legacy Fortran Applications: Strategy and Case Study[END_REF][START_REF] Reinders | High Performance Parallelism Pearls: Multicore and Many-core Programming Approaches[END_REF][START_REF] Taylor | Code modernization and modularization of APEX and SWAT watershed simulation models[END_REF]. This initiative is further promoted by hardware vendors such as the "Intel code modernization enablement program" 1 which provides tools and guidelines for the development of state-of-the-art and cutting-edge numerical simulation tools, capable of scaling on exascale machines. Moreover, the Exa2CT 2 project aims at facilitating the development of highly optimized scientific codes, capable to scale on exascale machines, through the development of new algorithms and programming techniques, validated on proto-applications, and that can be directly re-integrated into parent computational codes.

In the particular case of neutron transport simulations, a lot of efforts have been dedicated for improving the efficiency of the discrete ordinates method. These efforts can be classified in two fields: numerical methods and parallelization strategies that were developed accordingly. The most computationally demanding portion in the S N method is the space-angle problem called sweep operation, for each energy group. This operation acts like a wave front propagating throughout the spatial domain, according to the angular direction. Therefore, the successive parallel algorithms developed for the S N method focus on the parallelization of this sweep operation. The first parallel sweep algorithm KBA [START_REF] Baker | An SN Algorithm for the Massively Parallel CM200 Computer[END_REF] is implemented in numerous early S N solvers [START_REF] Gregory G Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF][START_REF] Humbert | Parallelization of PANDA Discrete Ordinates Code Using Spatial Decomposition[END_REF]. KBA splits the 3D spatial grid on a 2D process grid, enabling each process to perform the sweep on a local subdomain in a classical fork-join mode. The efficiency of this algorithm has been extensively studied in the literature through performance models [START_REF] James | Comparison via Parallel Performance Models of Angular and Spatial Domain Decompositions for Solving Neutral Particle Transport Problems[END_REF][START_REF] Hoisie | Performance and Scalability Analysis of Teraflop-Scale Parallel Architectures using Multidimensional Wavefront Applications[END_REF][START_REF] Kerbyson | Performance Modeling of Deterministic Transport Computations[END_REF][START_REF] Silber-Chaussumier | Recouvrement des communications et des calculs, du matériel au logiciel[END_REF]. Moreover, the advent of modern massively parallel computers has motivated extensions of this algorithm to provide sufficient concurrency for an improved scalability on large number of cores. In the unic code system [START_REF] Kaushik | Enabling high-fidelity neutron transport simulations on petascale architectures[END_REF], the authors implemented a parallel decomposition over space, energy and angle, enabling to extract more parallelism. A similar approach has been recently introduced in the Denovo code [START_REF] Thomas M Evans | Three-Dimensional Full Core Power Calculations for Pressurized Water Reactors[END_REF], where a new multilevel parallel decomposition allows concurrency over energy in addition to the space-angle. Furthermore, the pdt code [1] implements an extension of the KBA algorithm that enables a decomposition of the 3D spatial grid over a 3D process grid. These approaches enable to tackle neutron transport problems featuring larger number of energy groups. However, the parallelism in these codes adopts a uniform view of a supercomputer as collection of distributed computing cores (Flat approach), without explicitly addressing all the hierarchical parallelism provided by the modern architectures.

Recently, in the neutron transport community, some research initiatives have emerged to cope explicitly with the hierarchical topology (cluster of multiprocessor with multiple cores) of modern architectures. These initiatives are conducted through the development of protoapplications. For instance, the snap code, as a proxy application to the Partisn [START_REF] Baker | PARTISN on Advanced/Heterogeneous Processing Systems[END_REF] code focuses on exploring a Hybrid programming model for the sweep operation and auto-vectorization capabilities. Similarly, kripke is a proxy application to the ardra [START_REF] Aj Kunen | KRIPKE-A Massively Parallel Transport Mini-App[END_REF] radiation transport code, that is being developed to study the implications of new programming models and data layouts on its parent code performances.

As we have seen, the neutron transport simulation is a large and complex research area, in an international competitive environment. We record ourselves in this area and, in this thesis, we give some original contributions to the field. Indeed, we considered all the issues associated with the emerging architectures, that encompass the vectorization (SIMD), multithreading and message-passing paradigms. We make a key point on maximizing the sustainable peak performance of our implementation on distributed multicore-based machines using emerging software tools and frameworks (task-based runtime systems) that enable to reach higher performance, on today platforms, and the performance portability on future exascale architectures. More precisely we followed the trends on code modernization, and we developed highly-scalable parallel algorithms for the Cartesian S N neutron transport simulations over distributed multicore-based architectures. We especially tailored the discrete ordinates method to explicitly use all levels of parallelism available on these architectures. Furthermore, our approach relies at each level on theoretical performance models. In Chapter 1, we present background on reactor physics, the neutron transport equation, and we describe the sweep operation. Then, in Chapter 2, we discuss the vectorization strategies for maximizing the single-core peak performance of the sweep operation, and we justify our approach with a theoretical performance model. In Chapter 3, we discuss the design and theoretical performances of a new efficient parallel Hybrid sweep algorithm, targeting multicore-based architectures, using a new performance model which extends existing ones. In Chapter 4, we present an implementation of our Hybrid sweep algorithm using emerging task-based models on top of generic runtime systems. The integration of this sweep algorithm into our S N solver Domino is then presented in Chapter 5. A key point in this chapter is the implementation of a new efficient piece-wise acceleration method required to In this chapter, we introduce the neutron transport simulation as a requirement for enhancing the safety, efficiency and design of nuclear reactor cores. We start by presenting how neutrons inside a nuclear reactor core are modeled and the fundamental equation governing this process: the Boltzmann transport equation. Then, we derive the eigenvalue form of this equation that is used to determine the criticality of a reactor core. After that, we present the multigroup approximation that is used for the discretization of the energy variable and we introduce the discrete ordinates method, and the sweep operation that is used for solving the space-angle problem. Finally, we give an introduction to the modern parallel architectures on which the simulations will have to be executed.

Analysis of nuclear reactor cores

This section presents a general view on the reactor physics engineering and the mathematical model describing the neutron transport in the reactor core. This presentation is inspired from the works in [START_REF] Duderstadt | Nuclear Reactor Analysis[END_REF][START_REF] Lewis | Computational Methods of Neutron Transport[END_REF].

Basic concepts of nuclear reactor physics

A nuclear power plant is an industrial facility dedicated to electricity production, from the energy generated by the fissions of heavy nuclei (e.g. 235 U or 239 Pu) taking place in the nuclear reactor core. Each fission is induced by a neutron and releases an average energy of about 200 MeV [START_REF] Duderstadt | Nuclear Reactor Analysis[END_REF] in the form of heat, and some additional neutrons (2 or 3 on average) which can in turn induce other fissions, hence leading to a chain reaction (Figure 1.1). Therefore, the Figure 1.1: The fission chain reaction (illustration inspired by the CEA web site.) fission energy increases the temperature of the water1 circulating in the reactor core, also called a coolant. This high temperature water, by means of a heat exchanger, generates high pressure and high temperature steam from low temperature water. Finally, the steam is used to run turbogenerators producing the electricity. Hence, the more neutron fission, the more energy is released by the reactor. However, the fission probability of a nucleus is higher if the speed of the incident neutron is slow. Thus, in order to increase this fission probability, a moderator is generally used to slow down the neutrons2 . For a Pressurized Water Reactor (PWR), on which we will focus for the remainder of this dissertation, the coolant also acts as a moderator.

As we have seen, the neutrons play an important role in the analysis of nuclear reactor cores, because they determine the fission process and thus the power generated by a nuclear reactor [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]. In addition, the neutron population of the core is a fundamental information. This information is used to conduct safety studies and optimization of both fuel reload and reactor core design. It is therefore of great interest to evaluate precisely the neutron population (or neutron flux distribution) in each region of the reactor core #" r , at a given time t, as a function of the neutron's energy E, and according their propagation direction #" Ω. This represents therefore a complicated problem in which the phase-space has 7 dimensions.

The distribution of neutrons inside the core is influenced by the different types of nuclear reactions that neutrons can undergo. As already mentioned, a neutron when captured by a heavy nucleus can induce a fission of the nucleus into lighter nuclei. This fission results in more neutrons, energy release and some fission products. In addition, a neutron can simply be captured (pure absorption) by the nucleus without inducing a fission; or it can scatter off (bounce off) the nucleus. Finally the neutron can leak outside the reactor. Hence, evaluating the neutron distribution inside the core requires taking into account all these reactions all together. There are essentially two main families of methods used to simulate the neutron population: the probabilistic methods (Monte-Carlo) and the deterministic ones. The former consist of using statistical tools to characterize the life of a neutron from its "birth" to "death" (absorption by a nucleus or leaking out of the reactor), through a simulation of the history of a large number of particles (neutrons). The main advantage of the Monte-Carlo methods is that they avoid the phase-space mesh problems, because no discretization is required. Unfortunately, probabilistic methods allow phase-space density estimation with an accuracy that converges rather slowly with the number N of particles (∝ 1/ √ N ). The latter, deterministic methods, rely on finding a numerical solution of a mathematical equation describing the flow of neutrons inside the core. The fundamental equation governing the flow of neutrons is the Boltzmann Transport Equation (BTE), also called the neutron transport equation.

Boltzmann transport equation

In a nuclear reactor core, there are basically two different situations with regard to the interaction of a neutron with the nuclei present in the core.

• The neutron interacts with no nuclei; it therefore moves spatially at the same speed, without shifting from its initial direction: transport without collision.

• The neutron interacts with a nucleus. In this case, it can either be scattered by the nucleus, that is by changing its direction and its energy, or cause the fission of the target nucleus.

Therefore, the neutron transport equation or linear Boltzmann Transport Equation (BTE) is obtained by establishing a balance between arrivals and migrations of neutrons at the spatial position #" r , traveling with an energy E (or a speed v), toward direction #" Ω at a given time t. It is presented in equation (1.1), where: Σ t ( #" r , E, t) is the total cross-section. It characterizes the probability that a neutron of energy E interacts with a nucleus at the position #" r .

Σ s ( #" r , E ′ → E, #" Ω ′ • #" Ω, t)
is the scattering cross-section. It characterizes the probability that a neutron of energy E and direction #" Ω will be scattered by a nucleus and result in an outgoing neutron of energy E ′ and direction #" Ω ′ .

Σ f ( #" r , E, t) is the fission cross-section. It characterizes the probability that an interaction between a neutron of energy E with a nucleus at position #" r results in a fission.

ν is the average number of neutrons produced per fission.

χ(E) is the fission spectrum. It defines the density of neutrons of energy E produced from fission.

S ext ( #" r , E, #" Ω, t) is an external source of neutrons. It is generally used to start the chain reaction.

S 2 is the unit sphere.

1 v ∂ψ ∂t ( #" r , E, #" Ω, t) = - Transport #" Ω • - → ∇ψ( #" r , E, #" Ω, t) - Collision Σ t ( #" r , E, t)ψ( #" r , E, #" Ω, t) + Scattering ∫ ∞ 0 dE ′ ∫ S 2 d #" Ω ′ Σ s ( #" r , E ′ → E, #" Ω ′ • #" Ω, t)ψ( #" r , E ′ , #" Ω ′ , t) + Fission χ(E) 4π ∫ ∞ 0 dE ′ ∫ S 2 d #" Ω ′ νΣ f ( #" r , E ′ , t)ψ( #" r , E ′ , #" Ω ′ , t) + External source S ext ( #" r , E, #" Ω, t), (1.1)
Two different boundary conditions (BC) are generally used in reactor physics applications: vacuum and reflective boundary conditions. They are respectively defined in equation (1.2) and equation (1.3).

Vacuum BC:

ψ( #" r , E, #"

Ω, t) = 0 when #" Ω • #" n < 0, (1.2)
Reflective BC:

ψ( #" r , E, #" Ω, t) = ψ( #" r , E, #" Ω ′ , t) when #" Ω • #" n < 0 and #" Ω • #" n = - #" Ω ′ • #" n , (1.3)
where Γ is the boundary of the computational domain (the reactor core in our case), and #" n an outward normal to Γ. Vacuum BC represents the case where there is no incoming neutron from the outside of the core and is used when the full description of the core is given. Reflective BC is used to take advantage of the symmetries presented by the physical problem. In this dissertation, we will consider only the vacuum BC.

Let H and F be the transport and fission operators as defined by the following equations (1.4) and (1.5):

Hψ( #" r , E, #" Ω, t) = #" Ω • -→ ∇ψ( #" r , E, #" Ω, t) + Σ t ( #" r , E, t)ψ( #" r , E, #" Ω, t) - ∫ ∞ 0 dE ′ ∫ S 2 d #" Ω ′ Σ s ( #" r , E ′ → E, #" Ω ′ • #" Ω, t)ψ( #" r , E ′ , #" Ω ′ , t), (1.4) Fψ( #" r , E, t) = χ(E) 4π ∫ ∞ 0 dE ′ ∫ S 2 d #" Ω ′ νΣ f ( #" r , E ′ , t)ψ( #" r , E ′ , #" Ω ′ , t). (1.5)
Then, the equation (1.1) becomes:

1 v ∂ψ ∂t ( #" r , E, #" Ω, t) = -Hψ( #" r , E, #" Ω, t) + Fψ( #" r , E, t) + S ext ( #" r , E, #" Ω, t). (1.6)
For a nuclear reactor core under normal operating conditions, we seek to have a self sustained chain reaction in the absence of external sources of neutrons; this corresponds to an equilibrium state between fission neutron production and migrations of neutrons by transport or collision. In order to determine the state of the core in such a situation, it is necessary to find a nonnegative stationary solution of the Boltzmann equation (1.6) without the external source S:

Hψ( #" r , E, #" Ω) = Fψ( #" r , E).
However, for any given set of cross-sections, there generally exists no stationary solution to the BTE [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]. One way of transforming the previous equation so that it accepts a nonnegative solution is to adjust the average number of neutron procuced per fission, ν, so that a global time-independant balance can be preserved. Hence, we replace ν by ν/k and we obtain the following generalized eigenvalue problem:

Hψ( #" r , E, #" Ω) = 1 k Fψ( #" r , E, #" Ω), (1.7)
for which there will be a largest value of k such that a nonnegative solution exists. The largest such eigenvalue, which is also the spectral radius of operator H -1 F, is called the effective multiplication factor, and denoted k ef f . Physically, this coefficient allows to determine the criticality of the reactor core.

• If k eff < 1, (Hψ( #" r , E, #" Ω) > Fψ( #" r , E, #" Ω))
, the neutron production from fissions is less than migrations: the chain reaction turns off, and the reactor is said to be subcritical.

• If k eff > 1, (Hψ( #" r , E, #" Ω) < Fψ( #" r , E, #" Ω))
, the neutron production from fissions is larger than migrations: the reactor is said to be supercritical.

• If k eff = 1,
there is a strict equilibrium between production and migration of neutrons: the reactor is said to be critical.

problem (1.7) is solved using an inverse power algorithm, which leads to the computation of the neutron flux ψ and the eigenvalue k, by iterating on the fission term as presented in equation (1.8).

Hψ n+1 = 1 k n Fψ n , k n+1 = k n √ < Fψ n+1 , Fψ n+1 > < Fψ n , F ψ n > . (1.8)
Algorithm 1 describes the continuous form of the power iterations, also called external iterations.

The power algorithm converges slowly near the criticality. We therefore use the Chebyshev polynomials to accelerate its convergence. This is done by evaluating the neutron flux at iteration n + 1 as a linear combination of the solutions obtained in the last 3 iterations [START_REF] Saad | Chebyshev Acceleration Techniques for Solving Nonsymmetric Eigenvalue Problems[END_REF]. Each iteration of the power algorithm involves an inversion of the transport operator (Line 3 of Algorithm 1) and a source computation (Line 5 of Algorithm 1) which requires a discretization of the transport and fission operators. The next section describes the energetic discretization of these operators.

Multigroup formulation of the transport equation

The discretization of the energy variable is realized according to the multigroup formalism as described in [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]. In this formalism, the energy domain is split into a finite set of energy intervals, called energy groups and delimited by a decreasing sequence of carefully selected energy values Algorithm 1: Inverse power algorithm Input :

ψ 0 , νΣ f , Σ s , Σ t Output: k eff , ψ 1 S = Fψ; 2 while |k-k old | k old ≥ ϵ k or ||S-S old || ||S old || ≥ ϵ ψ do
      H 11 H 12 • • • H 1G H 21 H 22 • • • H 2G . . . . . . . . . . . . H G1 H G2 • • • H GG             ψ 1 ψ 2 . . . ψ G       =       S 1 S 2 . . . S G       , ( 1.9)
where:

ψ g ( #" r , #" Ω) = ∫ E g-1 Eg ψ( #" r , E, #" Ω)dE, (1.10)
is the multigroup angular flux of group g, and

H gg ψ g ( #" r , #" Ω) = #" Ω • #" ∇ψ g + Σ g t ( #" r , #" Ω)ψ g ( #" r , #" Ω) - ∫ S 2 d #" Ω ′ Σ g→g s ( #" r , #" Ω ′ )ψ g ( #" r , #" Ω ′ • #" Ω), H gg ′ ψ g ′ ( #" r , #" Ω) = - ∫ S 2 d #" Ω ′ Σ g ′ →g s ( #" r , #" Ω ′ • #" Ω)ψ g ′ ( #" r , #" Ω ′ ), g ̸ = g ′ S g ( #" r ) = 1 k n eff χ g 4π G ∑ g ′ =1 ∫ S 2 νΣ g ′ f ( #" r )ψ n g ′ ( #" r , #" Ω ′ )d #" Ω ′ , ( 1.11)
define the multigroup equations. The definition of the multigroup flux ψ g in equation (1.10) is justified by the energy separability hypothesis: the angular flux within the group g is approximated as a product of a function f and ψ g as defined by equation (1.12).

ψ( #" r , E, #" Ω) ≈ f (E)ψ g ( #" r , #" Ω), such that ∫ E g-1 Eg f (E)dE = 1.
(1.12)

In the multigroup formulation (equation (1.11)), the neutron parameters (cross-sections and fission spectrum), are group-wise defined. A comprehensive description of the evaluation of multigroup cross-sections can be found in [START_REF] Carlson | International Evaluation of Neutron Cross Section Standards[END_REF][START_REF] Said | Neutron Cross Sections: Neutron Resonance Parameters and Thermal Cross Sections Part B[END_REF]. We recall that the multigroup scattering cross-sections Σ g ′ →g s are expanded on the Legendre polynomials basis (see [START_REF] Hébert | Applied Reactor Physics[END_REF]). In the following, for the sake of clarity, all the numerical algorithms will be presented only for the order 0, corresponding to the case of an isotropic collision. In this case, the scattered neutron has a uniform probability to go along all directions in S 2 , and the scattering cross-section is given by:

Σ g ′ →g s ( #" r , #" Ω • #" Ω ′ ) = 1 4π Σ s0 ( #" r ). (1.13)
The resolution of the multigroup problem (1.9) is done using a block Gauss-Seidel (GS) algorithm as presented in Algorithm 2. It is worth noting here that the convergence of this Algorithm 2: Gauss-Seidel algorithm Input : ψ, S Output: ψ 1 while Non convergence do

2 for g ∈ 1, G do 3 Q = S g - ∑ g ′ ̸ =g H gg ′ ψ g ′ ( #" r , #" Ω); 4 H gg ψ g ( #" r , #" Ω) = Q;
GS algorithm depends on the sparsity profile of the transport matrix, which is determined by the scattering of the problem (presence of the coefficient Σ g ′ →g s in equation (1.11)). As the upscattering1 is only possible for some thermal groups2 (see [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]), almost all the non-zero elements of the transport matrix are located in its lower triangular part, including the diagonal.

Hence, it is not necessary for the GS algorithm to iterate for all energy groups: the resolution is direct for all fast groups, using a forward substitution, and the GS iteration apply only for the thermal groups. Each one of the GS iterations involves resolutions of G one-group space-angle problems:

H gg ψ m+1 g ( #" r , #" Ω) = - ∑ g ′ ̸ =g H gg ′ ψ m g ′ ( #" r , #" Ω) + S g ( #" r ), g = 1, . . . , G. (1.14)
We refer to these one-group problems as monokinetic equations. One should note here that, for each Gauss-Seidel iteration, the resolutions of all the monokinetic equations are necessarily done sequentially. This is an intrinsic property of the Gauss-Seidel algorithm. In order, for example, to parallelize the resolutions of these equations, one could use for example the Block-Jacobi algorithm as presented in [START_REF] Gregory G Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF]. But as our target applications do not feature many energy groups (less than 26), and because the Jacobi iterations converge two times slower than compared to Gauss-Seidel [START_REF] Randall | Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems[END_REF], we would rather keep using the Gauss-Seidel algorithm.

There are two major classes of resolution methods of monokinetic equations: P N and S N methods. The P N method consists of expanding the angular flux and the scattering crosssection on a truncated Legendre polynomials basis, whereas the discrete ordinates method (S N ) consists of discretizing the angular variable on a finite number of directions. We will focus on the latter method in the rest of this dissertation.

Discrete ordinates method on Cartesian meshes

In this section we present the discretizations of the angular and spatial variables of the monokinetic equations.

Angular discretization

Let us consider the monokinetic transport equation (1.14), on which both group indices and iteration indices are removed to simplify the notations. The angular dependency of this equation is resolved by looking for solutions on a discrete set of carefully selected angular directions { #"

Ω i ∈ S 2 , i = 1, 2, • • • , N dir }, called discrete ordinates: #" Ω i • #" ∇ψ( #" r , #" Ω i ) + Σ t ( #" r , #" Ω i )ψ( #" r , #" Ω i ) Lψ( #" r , #" Ω i ) - Rψ( #" r , #" Ω i ) ∫ S 2 d #" Ω ′ Σ s ( #" r , #" Ω ′ • #" Ω i )ψ( #" r , #" Ω ′ ) = Q( #" r , #" Ω i ) ∀i, (1.15)
where Q( #" r , #" Ω i ) gathers monogroup fission and inter-group scattering sources. Basically, the choice of the discrete ordinates is the same as to uniformly distribute a finite number of points on the unit sphere; which is not a trivial task because of the curvature of the sphere. In general, these discrete ordinates are determined thanks to a numerical quadrature formula. In Domino, we use the Level Symmetric quadrature formula, which leads to N dir = N (N + 2) angular directions, where N stands for the order of the Level Symmetric quadrature formula. Each angular direction is associated to a weight w j for integral calculation on the unit sphere S 2 , such that:

∫ S 2 g( #" Ω)d #" Ω ≃ N dir ∑ j=1 w j g( #" Ω j ), (1.16)
for any function g summable over S 2 . Hence, the scattering term Rψ( #" r , #" Ω i ) in equation (1.15) becomes:

Rψ( #" r , #" Ω i ) ≃ N dir ∑ j=1 w j Σ s ( #" r , #" Ω j • #" Ω i )ψ( #" r , #" Ω j ).
Using equation (1.13), we obtain:

Rψ( #" r , #" Ω i ) = 1 4π Σ s0 ( #" r )ϕ 00 ( #" r ),
where ϕ 00 ( #" r ) is the zeroth angular flux moment defined as:

ϕ 00 ( #" r ) = ∫ S 2 d #" Ω ′ ψ( #" r , #" Ω ′ ) ≃ N dir ∑ j=1 w j ψ( #" r , #" Ω j ).
One should note that the major drawback of the discrete ordinates method is the problem of so called "ray effects". Indeed, as shown in [START_REF] Kaye | Ray Effects in Discrete Ordinates Equations[END_REF], for some problems, such as those featuring a localized fixed-source in a pure absorber media, the discrete ordinates method can give a flux only for regions "surrounded" by the directions of the ordinates, taking their origin in the fixedsource. To remedy this issue, one solution is to increase the order of the angular quadrature used. However, in a nuclear reactor core, which is our focus in this work, the neutron sources are uniformly distributed in the whole core, and thus the probability of having non-covered regions is smaller. Therefore, our goals are not oriented towards numerical methods for eliminating these effects. The resolution of (1.15) is obtained by iterating on the in-scattering term R as presented in Algorithm 3. In highly diffusive media, the convergence of this algorithm is very slow, and therefore an acceleration scheme must be combined with this algorithm in order to speed-up its convergence. Section 5.1 presents a new acceleration scheme used in Domino to improve the convergence of scattering iterations. As presented in Algorithm 3, each scattering iteration involves the resolution of a fixed-source problem (Line 3), for every angular direction. This is done by discretizing the streaming operator which is presented in following section.

Algorithm 3: Scattering iterations

Input : ψ k Output: ψ k+ 1 2 1 while Non convergence do 2 Rψ k ( #" r , #" Ω) = ∫ S 2 d #" Ω ′ Σ s ( #" r , #" Ω ′ • #" Ω i )ψ k ( #" r , #" Ω ′ ); 3 3 Lψ k+ 1 2 ( #" r , #" Ω) = Rψ k ( #" r , #" Ω) + Q( #" r , #" Ω);

Spatial discretization

Let us consider the fixed-source monokinetic equation along the angular direction

#" Ω = (Ω x , Ω y , Ω z ) (Line 3 of the Algorithm 3): Ω x ∂ψ ∂x (x, y, z) + Ω y ∂ψ ∂y (x, y, z) + Ω z ∂ψ ∂z (x, y, z) + Σ t ψ(x, y, z) = B(x, y, z), (x, y, z) ∈ V ijk
(1.17) where B gathers all neutron sources including scattering and fission; the angular variable is omitted to lighten the notations. In this work, we focus on a 3D reactor core model, represented by a 3D Cartesian domain D. The spatial variable of this equation is discretized using a diamond difference scheme (DD), as presented by A. Hébert in [START_REF] Hébert | High order diamond differencing schemes[END_REF]. In particular, the DD0 scheme, as implemented in Domino, is derived by combining the moment of order 0 of the transport equation with some closure relations.

Let us first define a map from the cell

V ijk = [x i , x i+1 ] × [y j , y j+1 ] × [z k , z k+1 ], of size ∆x i × ∆y j × ∆z k , to the reference mesh V ref = [-1, 1] × [-1, 1] × [-1, 1] as follows: M : V ijk → V ref    x y z    →     x = 2x-(x i +x i+1 ) ∆x i ŷ = 2y-(y j +y j+1 ) ∆y j ẑ = 2z-(z k +z k+1 ) ∆z k     .
Using this map, the equation (1.17) is rewritten in:

2Ω x ∆x i ∂ψ ∂x (x, y, z) + 2Ω y ∆y j ∂ψ ∂y (x, y, z) + 2Ω z ∆z k ∂ψ ∂z (x, y, z) + Σ t ψ(x, y, z) = B(x, y, z) (x, y, z) ∈ V ref .
(1.18)

Moment of order 0 of the transport equation. The moment of order 0 of the transport equation is obtained by integrating the equation (1.18) on the cell V ref :

Σ t ∫ 1 -1 ∫ 1 -1 ∫ 1 -1 ψ(x, y, z)dxdydz + 2Ω x ∆x i ∫ 1 -1 ∫ 1 -1 (ψ(1, y, z) -ψ(-1, y, z)) dydz + 2Ω y ∆y j ∫ 1 -1 ∫ 1 -1 (ψ(x, 1, z) -ψ(x, -1, z)) dxdz + 2Ω z ∆z k ∫ 1 -1 ∫ 1 -1 (ψ(x, y, 1) -ψ(x, y, -1)) dxdy = ∫ 1 -1 ∫ 1 -1 ∫ 1 -1 B(x, y, z)dxdydz. (1.19)
To lighten the notations, we define 2 volumetric moments,

ψ 000 = 1 8 ∫ 1 -1 ∫ 1 -1 ∫ 1 -1 ψ(x, y, z)dxdydz, B 000 = 1 8 ∫ 1 -1 ∫ 1 -1 ∫ 1 -1 B(x, y, z)dxdydz,
and 6 surface flux moments (3 incoming and 3 outgoing),

ψ X±,00 = 1 4 ∫ 1 -1 ∫ 1 -1 ψ(±1, y, z)dydz, ψ Y ±,00 = 1 4 ∫ 1 -1 ∫ 1 -1 ψ(x, ±1, z)dxdz, ψ Z±,00 = 1 4 ∫ 1 -1 ∫ 1 -1 ψ(x, y, ±1)dxdy.
Then equation (1.19) becomes:

Ω x ∆x i ( ψ X+,00 -ψ X-,00 ) + Ω y ∆y j ( ψ Y +,00 -ψ Y -,00 ) + Ω z ∆z k ( ψ Z+,00 -ψ Z-,00
)

+ Σ t ψ 000 = B 000 ,
(1.20) where ψ X-,00 , ψ Y -,00 , ψ Z-,00 are known thanks to the boundary conditions. However, 4 unknowns have to be determined according to this equation: the volumic angular flux ψ 000 and the 3 outgoing surface moments ψ X+,00 , ψ Y +,00 , ψ Z+,00 . Indeed, since it is not possible to determine 4 unknowns from a single equation, we must combine equation (1.20) with 3 other equations. These equations are the closure relations provided by the diamond differencing scheme.

Closing relations These relations are obtained by cancelling the first term appearing in the moment of order 1 of the transport equation. The closing relations for the DD0 scheme are given as follows:

       ψ X+,00 = 2ψ 000 -ψ X-,00 ψ Y +,00 = 2ψ 000 -ψ Y -,00
ψ Z+,00 = 2ψ 000ψ Z-,00 (1.21) The equations (1.20) and (1.21) are solved by "walking" step by step throughout the whole spatial domain and to progressively compute angular fluxes in the spatial cells. In the literature, this process is known as the sweep operation. The whole S N algorithm as implemented in Domino is presented in Algorithm 4. In this algorithm, the sweep operation (Line 15) gathers the vast majority of computations performed. This operation is the focus of this dissertation, and we are going to describe it in detail in the following section.

Cartesian transport sweep operation

Without loss of generality, we consider the example of a 2D spatial domain discretized into 6 × 6 cells. The DD0 scheme in this case defines 3 DoFs per spatial cell: 1 for the neutron flux and 2 for the neutron current. The sweep operation is used to solve the space-angle problem defined by equations (1.20) and (1.21). It computes the angular neutron flux inside all cells of the spatial domain, for a set of angular directions. These directions are grouped into four quadrants in 2D (or eight octants in 3D). In the following, we focus on the first quadrant (labeled I in Algorithm 4: Discretized algorithm of the S N method as implemented in Domino

Input : νΣ f , Σ s , Σ t Output: k eff , ψ 1 ▷ Initialization of external iterations 2 ϕ = ( 1, 1, • • • , 1 ) t ; 3 C = G ∑ g=1 νΣ f,
S =       χ 1 χ 2 . . . χ G       • C 7 ▷
⃗ Ω k . ⃗ ∇ψ k + Σψ k = Q ∀ ⃗ Ω k ∈ S N ; 16 ϕ g = ∑ Ω k ω k ψ k ; 17 C old = C; 18 C = ∑ g νΣ f,g .ϕ g ; ▷ Fission sources update 19 k = √ <C,C> <C old ,C old > ; 20 C = 1 k C; Figure 1.2a).
As shown in Figure 1.2b, each cell has two incoming dependencies ψ L and ψ B for each angular direction. At the beginning, incoming fluxes on all left and bottom faces are known as indicated in Figure 1.2c. Hence, the cell (0, 0) located at the bottom-left corner is the first to be processed. The treatment of this cell allows the updating of outgoing fluxes ψ R and ψ T , that satisfy dependencies of the cells (0, 1) and (1, 0). These dependencies on the processing of cells define a sequential nature throughout the progression of the sweep operation: two adjacent cells belonging to successive diagonals cannot be processed simultaneously. Otherwise, treatment of a single cell for all directions of the same quadrant can be done in parallel. Furthermore, all cells belonging to a same diagonal can be processed in parallel. Hence, step by step, fluxes are evaluated in all cells of the spatial domain, for all angular directions belonging to the same quadrant. The same operation is repeated for all the four quadrants. When using vacuum boundary conditions, there is no incoming neutron to the computational domain and therefore processing of the four quadrants can be done concurrently. This sweep operation is subject to numerous studies regarding design and parallelism to reach highest efficiency on parallel architectures. In Chapter 4, we will present our task-based approach that enables us to leverage the full computing power of such architectures. boundary conditions (c) Domain decomposition and boundary conditions. The corner cell (0, 0) is the first to be processed for a quadrant, and its processing releases computations of its neighbors (Here (0, 1) and (1, 0)).

Modern parallel computers and performance evaluation

The landscape of today's high-performance computing capability includes a large number of different and powerful computing devices. These are essentially multicore processors and Graphics Processing Units (GPUs). While the former can be considered as an evolution of the classical processors, the latter are of a new kind because of the programming model shift required.

Although the focus of this dissertation is on multicore-based architectures, one should keep in mind that this work provides an opening towards GPU-based architectures.

Architecture and design of parallel computers

The advent of the multicore processor In the area of computer architecture, Moore's Law [START_REF] Falk | Law and Social Change: The Semi-Autonomous Social Field as an Appropriate Subject of Study[END_REF] 1 had allowed for the production of powerful monolithic processors until the beginning of the 2000s. The increase of the computing power of the processors was thus mainly due to the increase of the clock rate. However, as the transistors get smaller, their power density2 does not scale down below a threshold as it was previously predicted by the Dennard scaling [START_REF] Robert H Dennard | Design of Ion-Implanted MOSFET's with Very Small Physical Dimensions[END_REF][START_REF] Hadi Esmaeilzadeh | Dark Silicon and the End of Multicore Scaling[END_REF]. The consequence of this "power wall" (limits on temperature of the chip and its power consumption) has led processor designers to shift toward multicore processors to keep the growth of the processor performance. Furthermore, to increase the power efficiency (Performance/Watt), processor cores integrate new functional units that are capable of performing a Single Instruction on Multiple Data (SIMD) that will be discussed in Chapter 2. These units, also called vector units, are now common on all modern processors and will continue to be present on processors henceforth.

Modern distributed memory clusters

Building ever more powerful computers can be done by aggregating either several multicore processors (sockets) in a single computing node, sharing the main memory of the computer, or several computing nodes interconnected by a high-speed network. In the first case, we obtain either a Symmetric Multi-Processor (SMP) if all the CPU cores on the processor chip have the same access latency to the main memory, or a Non Uniform Memory Access (NUMA) computing node, providing hierarchical access to the memory. Currently, modern distributed computers are built from an interconnection of NUMA nodes, and therefore the parallelism on these computers is decoupled into four parts:

• SIMD units in each CPU core;

• CPU cores in each socket;

• sockets in each node;

• nodes of a supercomputer.

Consequently, as already noted "the free lunch is over"3 , and taking advantage of the increasing computing power brought by multicore-based architectures, requires a deep shift on the traditional programming models used until now. Indeed, these architectures are generally addressed by means of either multithreading or message-passing techniques. While the first solution is limited to shared-memory systems, the second one can be used both on shared and distributed memory systems. A combination of both approaches can also be used. In Chapter 4, we will discuss these solutions in detail.

Metrics for performance evaluation

We primarily use parallel computers to speed-up the computation time required for executing a given workload. To evaluate the performance of the workload on a given parallel computer, the main metric is the elapsed time which allows to determine the parallel efficiency of the workload.

Parallel efficiency This metric is related to the speed-up that can be achieved when running an application on a parallel computer. We consider the following notations:

• T s the serial computation time;

• T p the parallel computation;

• P the number of parallel processing units in use.

Then, the speed-up (S) and parallel efficiency (E) are defined by:

S = T s T p , E = S P . (1.22)
However, even though the parallel efficiency metric determines the scalability of an application, it does not indicate if the considered application is efficiently using the full computing power of the computer. To achieve this, we rather use another metric.

Flop/s This metric determines the number of floating point operations (Flop) executed per second (s) on a given machine, which can be considered as the "speed of execution" of a computational kernel on a given machine. Each computational kernel is associated with a number of Flop which evaluation can be done either by counting all the floating point operations that exist in the kernel, or using performance monitoring tools that collect hardware performance counter events during the execution of the computational kernel. In general, the Flop/s metric characterizes the performance of an application on a given machine and it should be compared to the theoretical peak performance of the target computer, which is the maximum number of Flop that a computer can perform per second.

▶ Having defined these metrics, we will use them throughout this dissertation to quantitatively evaluate the performances of our implementation on the target architectures. As mentioned, the point we made in this thesis is to tailor the discrete ordinates method for emerging architectures. In order to meet this goal, the first step in our process is the vectorization of the sweep operation as presented in Chapter 2.

Chapter 2

On the Vectorization of the Sweep Operation Contents As we mentioned in Chapter 1, on the way to designing a highly efficient massively parallel neutron transport solver targeting modern multicore-based supercomputers, we should take care of all the microarchitectural improvements brought by these computers. Since the early 2000s and the advent of the multicore processor, the main processor manufacturers have designed and introduced new vector execution units into the CPU cores, capable to apply a Single Instruction on Multiple Data (SIMD) residing in dedicated registers. The consequence of adding vector units to the CPU cores is a large boost of the corresponding processors peak performance. This SIMD execution model, also called vectorization, can therefore be used to reduce the run time of computational workloads, and thus to maximize the sustainable peak performance [START_REF] Victor W Lee | Debunking the 100X GPU vs. CPU Myth: an Evaluation of Throughput Computing on CPU and GPU[END_REF] and the energetic efficiency of these workloads [START_REF] Firasta | Intel AVX: New Frontiers in Performance Improvements and Energy Efficiency[END_REF]. In light of this observation, there is more and more research on designing new algorithms and re-factoring legacy scientific codes to enhance their compatibility to vectorization. We can cite the works in [START_REF] Kreutzer | A Unified Sparse Matrix Data Format for Efficient General Sparse Matrix-Vector Multiplication on Modern Processors with Wide SIMD Units[END_REF] where the authors describe some optimizations to help efficiently use vector units for the classic operation of sparse matrix-vector multiply consisting in designing a vector-friendly storage format. Moreover, in [START_REF] Chhugani | Efficient Implementation of Sorting on Multi-Core SIMD CPU Architecture[END_REF], authors presented a design of a vectorized implementation of the MergeSort algorithm. In [START_REF] Henretty | Data Layout Transformation for Stencil Computations on Short-Vector SIMD Architectures[END_REF] the authors presented an efficient vectorized implementation of a stencil computation. For this reason, we propose in this chapter to study the vectorization capabilities offered by the most computational demanding operation in the discrete ordinates method: the sweep operation, in order to maximize its sustainable single-core performance.

The remainder of this chapter is organized as follows: in section 2.1, we recall the vector (or SIMD) programming model and how it can be used on modern multicore processors. Then, we give in section 2.2 a theoretical analysis of the sweep operation by the means of its computational intensity. In section 2.3 and section 2.4, we present two different strategies that can be used to vectorize the sweep operation.

Review of the SIMD paradigm

This section recalls the SIMD programming model from early vector processors to modern multicore processors.

General presentation of the SIMD paradigm

Historical trends

In Flynn's taxonomy of computer architectures [START_REF] Michael | Some Computer Organizations and Their Effectiveness[END_REF], Single Instruction Multiple Data (SIMD) corresponds to a class of computers that can perform a single instruction on a set of data stream. Such computers are also called vector processors and the process of applying a SIMD instruction is commonly called vectorization. According to the study in [START_REF] Espasa | Vector architectures: past, present and future[END_REF], vector processors can be classified into two categories: memory-to-memory and memory-to-register. In the first category, corresponding to early vector processors, the vector functional units directly retrieves data from the memory, process it and write back the result to memory. memory-to-register vector processors have vector registers that are used to store data retrieved from the main memory, and on which vector instructions operate. In practice, early vector processors used to exhibit a high ratio of bandwidth over peak performance and the speed of vectorized algorithms was optimal for basic BLAS1 like operations on long vectors. On modern chips, the situation is totally different and we will see that the performance of most vectorized algorithms is bounded to the memory-to-register bandwidth.

SIMD on modern CPUs

Increasing the processor clock rate, or frequency, has been, until the end of the last century, the major factor for improving the floating point performance of processors. However, since the early 2000s this frequency scaling stalled in favor of multicore chips allowing to have powerful processors while minimizing their power consumption. Furthermore, to increase the computing power of a single core, processor designers tend to add larger vector units into the processor cores. These units are capable to apply SIMD instructions, which are extensions of the processor Instruction Set Architecture (ISA), on fixed-size vector registers. Thereby, these vector units can be used to speed-up computations by a factor theoretically equal to the SIMD width of the considered processor. Table 2 stood as the maximum theoretical speed-up attainable on that architecture. Vector instructions comprises in addition of the classic floating point arithmetics (add/sub/mul/div), other types of instructions such as gather and scatter, or logical operations. Each of these instructions is characterized by a latency and a throughput1 , determined by the considered architecture. On Westmere, the load and store instructions, from CPU L1 cache to SIMD registers, have a throughput of 16 Bytes/cycle. However, as noted in [START_REF] Gwennap | Sandy Bridge Spans Generations[END_REF], on Sandy Bridge, which has doubled the SIMD width from 128 bits to 256 bits, the load throughput has doubled to 32 Bytes/cycle compared to Westmere, whereas the store throughput remains the same at 16 Bytes/cycle. Furthermore, in [START_REF] Vladimirov | Arithmetics on Intel's Sandy Bridge and Westmere CPUs: not all FLOPS are created equal[END_REF], authors show that the Intel Sandy Bridge microarchitecture gives no performance gain for the division as compared to the Intel Westmere microarchitecture, on a class of benchmarks. Thereby, for a code bottlenecked by division, the theoretical 2-fold speed-up when moving from Westmere to Sandy Bridge, to which one may expect, can not be observed. On the Skylake microarchitecture2 , it is possible to perform up to 16 floating point instructions using a single vector instruction. Consequently, without the usage of SIMD units, we can only take advantage of 1/16 of the CPU peak performance. For this reason, it is essential that scientific codes manage to make usage of these units in order to maximize their sustainable floating point performance. However, the performance of a perfectly vectorized algorithm will generally depend on the arithmetic intensity of this algorithm.

The critical arithmetic intensity issue

As mentioned previously, successive supercomputer generations brought a regular and impressive improvement of their peak performance. Surprisingly enough, the computer bandwidth that measures the maximal data flow between the computer RAM (off-chip bandwidth) and the floating point units (FPUs) did not increase as fast as the peak performance. The consequence of the broadening gap between the off-chip bandwidth and the chip peak performance is to put a dramatic emphasis on the arithmetic intensity of computational kernels and defined by:

I a =
Number of floating point operations Main memory traffic (Read+Write) .

It measures the average number of instructions executed by the CPU per byte read from or written to the main memory, and should be compared to the processor-dependent constant critical arithmetic intensity, defined by:

I c = Peak floating point operations Off-chip Bandwidth .
If the arithmetic intensity I a of a given kernel is lower than the critical value I c , then its performance does not depend on the computational power of the target processor, but mainly on the system memory bandwidth: the algorithm is then said to be memory bound. Consequently, the impact of the vectorization on such a kernel, with a low arithmetic intensity is negligible. Contrariwise, if I a is higher than I c , then the computational kernel is said to be cpu bound, and it therefore can benefit from the speed-up of computations enabled by the vectorization. The critical arithmetic intensity of processors increases with each generation causing an ever larger fraction of algorithms to be memory bound. Indeed, on Table 2 single-socket and the whole node, as the peak sustainable DRAM bandwidth varies with the number of cores in use. We used the Stream benchmark [START_REF] Mccalpin | STREAM: Sustainable Memory Bandwidth in High Performance Computers[END_REF][START_REF] Mccalpin | Memory Bandwidth and Machine Balance in Current High Performance Computers[END_REF] to evaluate the bandwidth in three cases1 . We found that the single-core I c of the considered processors are 3.6 Flop/Byte for bigmem and 3.4 Flop/Byte for a node of the athos platform. For a single-socket these values are of 8.3 Flop/Byte on the bigmem node compared to 15.0 Flop/Byte, on the more recent athos platform highlighting the increase of the critical arithmetic intensity for successive processor generations. As a consequence of this processor evolution, the whole design of our sweep implementation aimed at maximizing the arithmetic intensity in order to exploit the full power of modern multicore processors, thanks to the vectorization.

On the way for vectorization

Exploitation of vector units can be done by relying on the compiler auto-vectorizer. However, while this may work for simple loops, it is generally not the case for complex kernels. In [START_REF] Maleki | An Evaluation of Vectorizing Compilers[END_REF], the authors give several reasons on why the compiler fails in generating vector code. For instance, compilers generally lack accurate interprocedural analysis, which would help to enable important transformations such as changing the memory layout of data structures needed to vectorize the code. For example, such a memory layout change is required when data are accessed with nonunit stride accesses. Meanwhile, by explicitly using assembly instructions or compiler intrinsics corresponding to the target architecture, one can exploit these vector units with the expense of some hardware constraints specific to vector instructions. In fact, vector instructions operate on packed data loaded inside specialized registers of fixed size. To perform fast load and store operations, data items need to be well aligned on cache boundary: 16 bytes for sse, 32 bytes for avx and 64 bytes for avx-512. Sometimes we have to resort to padding to satisfy this requirement. As an example, when loading 256 bits packet data with Intel avx in a contiguous memory region of size 256+32×3 = 352 bits, we need to extend this region with 512-352 = 160 bits at the memory allocation stage. Attention should be paid to padding as it increases global memory consumption and useless computations. The drawback of inlining compiler intrinsics into the code is that it degrades the readability of the code and is error prone. To overcome this, one can rely on C++ generic programming concepts by overloading the arithmetic operators to call the corresponding compiler intrinsics.

Generic programming and tools for the vectorization

One solution to get a vectorized code is to use external libraries offering optimized instructions for a given architectures. Examples of these libraries include Intel MKL [START_REF] Wang | Intel Math Kernel Library[END_REF] or ACML [3]. However, these libraries are especially designed and tuned for linear algebra routines and thus can not be easily used in other contexts. It is also possible to rely on the compiler auto-vectorizer, which can enforce the vectorization of loops via a set of hints dictated by the programmer. Thus, the pragma simd extension, firstly introduced by Intel Cilk Plus [START_REF] Michael D Mccool | Structured Parallel Programming: Patterns for Efficient Computation[END_REF][START_REF] Arch D Robison | Cilk plus: Language Dsupport for Thread and Vector Parallelism[END_REF] and integrated in gcc and icc compilers, tells the compiler that a given loop can be vectorized. The performance that can be obtained from such an approach depends on the compiler in use [START_REF] Estérie | Boost.SIMD: Generic Programming for portable SIMDization[END_REF]. There are also some SIMD-enabled languages which feature some extensions for writing explicit vectorized code. This is the case with the Intel Cilk array notations [START_REF] Michael D Mccool | Structured Parallel Programming: Patterns for Efficient Computation[END_REF][START_REF] Arch D Robison | Cilk plus: Language Dsupport for Thread and Vector Parallelism[END_REF] where A[:] states that an operation on the array A is a vector instruction. The Intel SPMD Compiler (ISPC) [START_REF] Pharr | ISPC: A SPMD Compiler for High-Performance CPU Programming[END_REF], an LLVM-based language and compiler similar to CUDA/OpenCL, allows to automatically vectorize a C-based code, by mapping several SPMD program instances to SIMD units. On the other hand, one can use an external library that calls to the compiler intrinsics of the target architecture, and providing a classical API to the application developer that enhances the portability of the code. For instance, Boost.SIMD [START_REF] Estérie | Boost.SIMD: Generic Programming for portable SIMDization[END_REF] is a C++ template library which provides a class holding packed data residing in a SIMD register, along with a set operations overloaded to call the corresponding compiler intrinsics. Eigen [48] is a similar library, but unlike Boost.SIMD, it provides a Map object that can be used to interface with raw buffers. Thereby, it is possible to easily and elegantly take advantage of the vectorization of operations offered by Eigen, in a code using other containers than those of Eigen. This is the case for our Domino code, which uses the Legolas++ [START_REF] Kirschenmann | Vers des noyaux de calcul intensif pérennes[END_REF], a generic C++ library internally developed at EDF. This library provides basic constructing blocs to build Linear Algebra Solvers. In particular, it allows to write algorithms that apply indifferently on arrays of scalar and on arrays of packs of scalars.

Legolas++ is built on top of Eigen for SIMD issues and Intel TBB [START_REF] Pheatt | Intel® Threading Building Blocks[END_REF] for multi-threading.

Hence, in the following, we rely on the Eigen library to vectorize the sweep kernel, allowing us to benefit from new vector instructions of the next generation of processors without having to modify the code. As an example, the same source code will automatically be compiled for sse or avx machine instructions depending on the C++ compiler options. Note that Eigen, internally invokes the SIMD instructions explicitly and that the vectorized binary performance do not depends on the auto-vectorization capability of the C++ compiler.

Arithmetic intensity of the sweep kernel

In this section, we evaluate the arithmetic intensity of the sweep operation, as a function of the spatial and angular discretizations, in order to characterize the performance of the sweep operation. We first consider the case of a single cell with a single direction.

Memory traffic and flops per cell and per direction

We recall that in the sweep operation, the processing of a single spatial cell c ijk for a single angular direction #" Ω, consists in updating: the scalar flux at the cell-center and the neutron current on outgoing faces. This processing is depicted on Figure 2.1 (in DD0) and described on Algorithm 5, where ψ 0 u and ψ 1 u represent incoming and outgoing neutron current along the u dimension (u = x, y, z). To evaluate the arithmetic intensity of Algorithm 5, we need to determine the flops and the number of memory accesses per cell, per angular direction.

Memory accesses

According to the Algorithm 5, on each spatial cell, for each angular direction, the update of the outgoing neutron current requires to load:

• 12 input scalar values according to the quadrature formula used: inverse of the mesh steps δ u , u = x, y, z; source term S; total cross-section Σ t ; outgoing neutron current ψ 0 u , u = x, y, z; coordinates and weight of the angular direction (ν, η, ξ), ω;

• 1 input/output value: scalar flux ϕ;

Algorithm 5: Sweep of a single spatial cell (in DD0), for a single angular direction

In : δx = 2 ∆x ; δy = 2 ∆y ; δz = 2 ∆z ; S; Σ t ; ϕ; {ψ 0 u |u = x, y, z}; #" Ω ≡ {(ν, µ, ξ), ω} Out: ϕ; {ψ 1 u |u = x, y, z}; 1 ϵ x = ν δx; ϵ y = η δy; ϵ z = ξ δz; 2 ψ = ϵxψ 0 x +ϵyψ 0 y +ϵzψ 0 z +S ϵx+ϵy+ϵz+Σt ; 3 ψ 1 x = 2ψ -ψ 0 x ; 4 ψ 1 y = 2ψ -ψ 0 y ; 5 ψ 1 z = 2ψ -ψ 0 z ; 6 ϕ = ϕ + ψ • ω;
• and 3 output values: outgoing neutron current ψ1 u , u = x, y, z.

This amounts to a total of 16 read/write memory accesses. As a matter of fact, the storage of the neutron current is not necessary for a stationary computation. This reduces the drop of the number of memory accesses to 13. In this case, the neutron currents will simply be denoted by ψ u , u = x, y, z.

Flops

The floating point operations performed on a cell consist of 7 additions, 3 subtractions, 10 multiplications and 1 division. Note that the quantities2 ∆u , u = x, y, z are pre-evaluated, allowing to save 3 divisions and 3 multiplications per cell. On modern processors, each of add/sub/mul operation account for 1 flop. The question of how many flops to count for one floating point division is a debatable one as the answer depends on how this operation is implemented on the target architecture. There exist however some tools that can help in evaluating the exact number of flops of a given kernel for a given architecture. Some of these tools, as likwid 1 or papi 2 , monitor events occurring in hardware performance counters of the target processor, in order to determine actual flops executed by the processor in a laps of time. Some other tools evaluate the flops with extended benchmarks. For instance, Lightspeed3 is a library featuring a set of routines for accurate flops counting of operations other than add/sub/mul. With this library, the number of flops of such an operation is equal to average value of the ratio between the execution time of the considered operation and the execution time of a multiplication (which counts as 1 flop), on the same architecture. According to our experiments with Lightspeed, we found that the division operation costs 4.8 flops (resp. 5.8 flops) on bigmem (resp. athos). However, for our studies, we want an architecture-independent measure of the flops of our sweep kernel, in order to compare its performance on the different architectures. Consequently, we have conventionally set the number of flops of a division to be 5. Therefore, the processing of a single spatial cell for a single angular direction cost 25 flops.

Thus, in single precision, the arithmetic intensity of the sweep of a single cell with one angular direction is 25 4×16 ≈ 0.39 if the neutron currents are stored, and 25 4×13 ≈ 0.48 otherwise. Therefore, it is found that these values are below I c of our target machines (see Table 2.2), indicating that it is not efficient to sweep a single cell with a single direction. We are going to increase it by grouping the cells into a MacroCell object.

Definition 1 (MacroCell).

A MacroCell is a structure that represents a 3D block of contiguous spatial cells, and some physical data associated to these cells as described on Listing 2.1. The number of cells aggregated into the x, y and z dimensions is respectively denoted by n x , n y , and n z .

General formula of the arithmetic intensity of the sweep

We consider the sweep of a MacroCell of size n x × n y × n z , using a quadrature formula having D angular directions. In the previous section, we have shown that there are 25 Flop/cell/direction. Therefore, in a general case, the number of flops is 25Dn x n y n z , and the associated number of memory accesses is given on Table 2.3. The sizes of variables are detailed as follows.

Storage of the currents

Variable Size (×4 Bytes)

{δ u |u = x, y, z} n x + n y + n z S n x n y n z Σ t n x n y n z ϕ n x n y n z (ν, µ, ξ) 3D ω D No {ψ u |u = x, y, z} D(n x n y + n x n z + n y n z ) Yes {ψ l u |u = x, y, z and l = 0, 1} D(3n x n y n z +n x n y +n x n z +n y n z )
Table 2.3: Memory accesses required to perform a sweep over a MacroCell of size n x × n y × n z using a quadrature formula comprising D angular directions.

• δ u : There are n u cells along the dimension u = x, y, z; each cell being associated to a mesh step. Thus, the total number of mesh steps are the sum of the number of cells for the three dimensions (x, y, z): n x + n y + n z .

• S, Σ t , ϕ: For each spatial cell, there is one of each; thus, the size of each of these variables is equal to the total number of cells in the MacroCell: n x n y n z .

• For each angular direction, there are:

-3 + 1 values which represent the coordinates and the weight associated with the direction.

-1 value per incoming face of the MacroCell, for the neutron currents, if they are not stored. The size of each MacroCell face is either n x n y , or n y n z or n x n z . If the neutron currents are stored, we must use 3 additional values per cell.

Hence, the general formula evaluating the arithmetic intensity for a sweep of a MacroCell, if all the neutron currents are stored, I s a is given by equation (2.1).

I s a = 25Dn x n y n z 4((n x + n y + n z ) + 3n x n y n z + D(3n x n y n z + n x n y + n x n z + n y n z ) + 4D) (2.1)
However, as mentioned in section 2.1.1, flops are cheaper than memory accesses on modern processors because of the gap between the processor peak performance and the off-chip memory bandwidth. Thus, to take advantage of all the power of modern processors, it is essential to save the memory bandwidth by reducing the unnecessary number of memory accesses. We consider this requirement and we do not store the neutron currents. In this case, the arithmetic intensity is given by equation (2.2).

I a = 25Dn x n y n z 4((n x + n y + n z ) + 3n x n y n z + D(n x n y + n x n z + n y n z ) + 4D) (2.2)
A numerical evaluation of these formulas is presented in Figure 2.4 on page 33.

Algorithm 5 shows that the computations performed on a single cell is similar for any two angular directions and for any two cells. The full sweep algorithm in 3D is presented on Algorithm 6, where the angular directions are grouped in octants (Line 1). Each octant has M Algorithm 6: The general sweep algorithm

1 forall o ∈ {1, . . . , 8} do 2 forall c ∈ Cells do 3 ▷ c = (i, j, k) 4 δx = 2 ∆x ; δy = 2 ∆y ; δz = 2 ∆z ; 5 forall #" Ω ∈ { #" Ω o d , d = 1, . . . , M } do 6 ▷ #" Ω ≡ {(ν, µ, ξ), ω} 7 ϵ x = ν δx; ϵ y = η δy; ϵ z = ξ δz; 8 ψ = ϵxψx+ϵyψy+ϵzψz+S ϵx+ϵy+ϵz+Σt ; 9 ψ x = 2ψ -ψ x ; 10 ψ y = 2ψ -ψ y ; 11 ψ z = 2ψ -ψ z ; 12 ϕ = ϕ + ψ • ω;
angular directions such that D = 8M . While the processing of a cell for all angular directions belonging to a single octant can be done once, the processing of two cells in a single step is only possible when the cells belong to a same diagonal plane. Therefore, two different strategies can be used to parallelize the sweep operation. The first strategy leads to parallelize the computations over angular directions (Line 5), while the second one leads to parallelize the computations over the cells on the same diagonal (Line 2). In the following, we are going to illustrate how to efficiently exploit both strategies in order to maximize the single-core performance of the sweep kernel, by vectorizing the computations.

Vectorization over spatial domain

In this section, we discuss the design and implementation issues of the spatial vectorization through performance modeling. We will use the following definitions.

Definition 2 (Front). For a fixed octant, we define a front of order f as the list of all cells (a, b, c) such that a + b + c = f .

Definition 3 (Step). A step processes a set of cells with a set of angular directions in a single

instruction. Each of these sets may be reduced to a singleton.

In the remaining of this section, without loss of generality, we focus on sse for which the SIMD width is packSize = 4. For the sake of clarity, we first consider a 2D spatial domain, which size is defined by n x = 5 and n y = packSize, as depicted in Figure 2.2. , where all the cells on the same front can be processed at once. Our aim here is to use vector instructions to parallelize the processing of all the cells belonging to the same front. However, the number of cells on the fronts is not constant. Hence, we cannot directly use vector instructions as they operate on fixed-size registers. For this reason, we must add padding cells to the fronts 0, 1, 2 and 5, 6, 7 in order to have fixed-size vectors that fit the SIMD width of the target processor. This strategy is illustrated on from the bottom-left to the upper-right corner, each step consists in the process of a block of packSize cells in vector mode. The neutron currents along the x dimension ψ x have the same size as in the scalar algorithm. Thus, they can be read once and kept in the SIMD registers until the sweep is finished. One should note that, as we do not store the neutron currents, ψ x is overwritten after each step, and thus it must be reset to the correct value before the processing of fronts of order 1, 2 and 3:

The algorithm

ψ x [packSize -1 -f ] = 0, f = 0, . . . , packSize -1,
where f is the front order. The situation is a little different along the y dimension. Indeed, the processing of the block of cells at step 0 requires to load data indexed by 0, 1, 2 and 3 from the ψ y buffer. This can be done using an aligned load provided that the buffer is correctly aligned to cache boundary. However, at the step 1, we have to load data indexed by 1, 2, 3 and 4: we cannot use anymore an aligned load to retrieve that data. This is a first limitation of this algorithm as an unaligned load operation has an overhead relative to an aligned load operation1 .

Maximum theoretical speed-up

The considered example requires n x n y = 20 scalar steps to process the sweep for one quadrant, in scalar mode. In vector mode, the number of steps is equal to the number of fronts, that is n x + n y -1 = 8. Thus, the maximum theoretical speed-up brought by this spatial SIMD is s max = 20/8 = 2.5. In general, the formula giving the speed-up is given by equation (2.3).

s max = packSize n x n x + packSize -1 , ( 2.3)
If n x is large, then the theoretical speed-up reaches its optimal value of packSize. In general, if n y > packSize, then the domain is split in slices of width equal to packSize along the y dimension. In this case, the efficiency formula becomes:

s max = n x n y (n x + packSize -1) ⌈ ny packSize ⌉ , ( 2.4)
where the ceiling takes into account the padding that should be used along the y dimension.

The generalization of the spatial SIMD algorithm in 3D is relatively straightforward. Indeed, we just need to loop over all the planes along the z dimension, and to vectorize the computations on each of these planes, which is exactly a 2D sweep as presented in section 2.3.1. Thus the theoretical efficiency in 3D remains also the same as in 2D (equation (2.4)).

To summarize, the spatial sweep algorithm as presented in this section, has only two limitations that are limiting its efficiency: unaligned load/store and additional padding cells. The performance results from a preliminary implementation were disappointing and have therefore confirmed the identified limitations. In the following section, we will present the angular vectorization strategy that allows us avoiding these limitations.

Vectorization over the angular variable 2.4.1 The algorithm

We consider that the neutron currents are not stored for maximizing the arithmetic intensity of the sweep. The vectorization of a cell processing during the sweep, along the angular directions, requires to move the loop over directions into the innermost level as presented Algorithm 6. Inside each spatial cell, we compute simultaneously several angular directions belonging to the octant currently being swept (the forall loop on Line 5 of Algorithm 6), using SIMD instructions which operate on packs of directions. Handling any angular quadrature order requires us to set up a padding system: for example, S 16 Level Symmetric angular quadrature formula gives 16(16 + 2) = 288 angular directions or 288/8 = 36 directions per octant. When using single precision with Intel avx, as 36 is not a multiple of 8, we perform 40 angular directions processing per spatial cell corresponding to an efficiency of 36/40 = 0.9 (see Table 2 of S 2 , where the vectorization does not give any improvement, the padding performance penalty is acceptable even for a quadrature formula featuring a small number of directions. For instance, using a S 4 quadrature, it is theoretically possible to speed-up the computations by a factor of 3 with sse. In addition, using product quadrature formulas such as Gauss-Legendre, we have more flexibility to reduce the efficiency loss due to the padding. Indeed, one can choose a combination of the number of azimuthal and polar directions that gives a total number of directions divisible by packSize.

As mentioned in section 2.1.2, we use vectorized instructions provided by the Eigen library [48] to enforce the vectorization of the sweep kernel while maintaining the readability of the code. Listing 2.2 shows a snapshot of this vectorization, which features some constructs of Eigen. Eigen::Array is a class representing a matrix. The type of the matrix elements is defined by the first template parameter. In our case, we use this matrix class for representing a block of contiguous data, which size is equal to the SIMD width packSize of the target architecture. As mentioned previously, Eigen::Map is a class mapping an existing array of data which can be declared as well-aligned.

In the following, we discuss the performance the sweep operation vectorized over angular directions on Intel Westmere and Ivy Bridge microarchitectures. The performance measurements were carried out on the bigmem and athos platforms (see Appendix A).

Study of the arithmetic intensity of the sweep

To assess the advantage of not storing the neutron currents, we have numerically evaluated the arithmetic intensity of the sweep operation, defined by equations (2.1) and (2.2), as a function of the MacroCell size and for a different angular quadratures. The results are presented on Figure 2.4. It shows the variation of the arithmetic intensity of the sweep as a function of the MacroCell size and for different angular quadratures. The first observation that should be made here is that when the neutron currents are stored (Figure 2.4a), then independently of the angular quadrature used, the arithmetic intensity of the sweep operation is below the critical threshold I c of both platforms. Thus, is this case, the sweep kernel is memory bound: its performance depends only on the sustainable bandwith of the target platform, and the vectorization will provide no performance gain. Contrariwise, when the neutron currents are not stored (Figure 2.4b), the arithmetic intensity of the sweep dramatically improves. Starting from a MacroCell of size 3×3×3, I a is over the critical treshold I c for both processors: the sweep kernel is then cpu bound, and its performance depends on the theoretical peak performance of the considered CPU. In this case, it makes more sense to vectorize the computations over angular directions. In the following, we are going to study the actual performance obtained on the test machines, showing the impact of vectorization.

Actual performances vs Roofline model

In this section, we first discuss the speed-up brought by the vectorized implementation of the sweep operation on Intel Westmere and Intel Ivy Bridge microarchitectures. Then, we compare the performances of the vectorized implementation against a theoretical performance model determined by the processor architecture: the Roofline model, as presented in [START_REF] Williams | Roofline: an Insightful Visual Performance Model for Multicore Architectures[END_REF].

Performances of the angular vectorization SSE Table 2.5 shows the performances of both scalar and angular-vectorized single-core implementations of a sweep, using the sse-enabled computing node bigmem. The performance measurements were obtained by averaging the computation time of ten successive runs. We observe that: first, when the spatial computational domain has a single cell, then the speed-up brought by the vectorization is only 1.14 using a S 4 quadrature, which is below the sustainable speed-up of 3 that is expected in single precision (see We observe that the computation time in scalar mode on the Intel Xeon E52697 V2 is slower than that was obtained on the Intel Xeon E78837. For instance, the sweep in scalar mode an a MacroCell of size 64 × 64 × 64 with a S 4 angular quadrature takes 82.21 ms on the Intel Xeon E78837 processor compared to 59.28 ms on Intel Xeon E52697 V2 processor, which represents an improvement of a factor of 27.8%. To justify this factor, one should consider the improvements brought by Ivy Bridge microarchitecture, including: higher DRAM bandwidth, larger cache sizes and the new µop cache for storing the decoded instructions1 . These new capabilities of the Ivy Bridge could therefore boost the throughput of instructions and thus the performances of a computational kernel. Although we have not quantified the individual effects of each of these features, the improvement we observed complies with the study in [START_REF] Jarp | Evaluation of the Intel Sandy Bridge-EP server processor[END_REF], where the authors observed an average speed-up of 33% over a series of benchmarks, between Westmere and Sandy Bridge architectures.

S 4 S 8 S 12 S 16 n x T scalar T sse (×) T scalar T sse (×) T scalar T sse (×) T scalar T sse (×) (ms) (ms) (ms) (ms) (ms) (ms) (ms) (ms
The second observation is that, as on the Westmere node, the speed-up brought by the vectorization increases with the MacroCell size and the angular quadrature. However, this speed-up is lower than expected. For instance, let us consider the performance of a sweep of a MacroCell of size 64 × 64 × 64 with a S 16 quadrature. Using avx we got a speed-up of 4.32, relative to a scalar execution, which is less than the expected speed-up of 7.2 (see Table 2.4). One should keep in mind that, as explained in [START_REF] Vladimirov | Arithmetics on Intel's Sandy Bridge and Westmere CPUs: not all FLOPS are created equal[END_REF], a floating point division performs slightly worse on Sandy Bridge than on a Westmere. Thus, if a avx-vectorized code uses floating point divisions, then its performance may not correspond to the optimal speed-up that could be expected.

Roofline model

Here, we compare the actual single-core performance of the vectorized implementation of the sweep operation against the Roofline model [START_REF] Williams | Roofline: an Insightful Visual Performance Model for Multicore Architectures[END_REF] which defines the maximum performance attainable on a given multicore architecture. This Roofline model uses two parameters (the sustainable DRAM bandwidth B, and the theoretical peak performance P of the target processor) to evaluate the maximum attainable performance by a kernel, as a function of its arithmetic intensity I a . The evaluation of this Roofline model is defined by equation (2.5). Performance = min (P, I a * B) .

(2.5)

In the following, as we are focusing on the single-core performances, we consider a particular case of the Roofline model by evaluating its parameters for a single core. Figure 2.6 shows the performance variation of the vectorized implementation of the sweep on the Intel Xeon E78837 processor, as a function of the arithmetic intensity I a , and for different quadrature order. As mentioned previously, we observe that the performance of the sweep in- creases with I a and the order of angular quadrature used. Furthermore, the actual performances follow the same trend as the predicted performance from the Roofline model, and reach a plateau the value of which depends on the quadrature: for a S 16 , the plateau is 13.6 Gflop/s, corresponding to 62.9% of the theoretical peak performance per core of the bigmem node. One should note that the saturation of actual performances occurs a little later after the ridge point indicated by the Roofline model. In addition, the discrepancies between actual performances and the upper bound given by the Roofline model, are larger for lower values of I a (< 20 Flop/Byte). For instance, using a S 16 quadrature, the ratio between the prediction of the Roofline model and the actual performance of the sweep is 20.9 for I a = 0.89 and of 1.5 when I a = 154.6. Indeed, our evaluation of the flops for the sweep operation (25 flops per cell per angular direction) takes into account only floating point operations performed in the computational kernel, and thus ignores all the remaining operations performed in the code. This hypothesis is valid, provided that the MacroCell in use is large enough to provide a large amount of computations in the kernel to counterbalance the remaining operations in the code. Consequently, for low values of I a the flops are underestimated which explains the large discrepancy between the Roofline model and actual measurements.

The same study was performed on the Intel Xeon E52697 V2 processor (Figure 2.6), and yields similar conclusions. One should note that the saturation of the performances on the Intel Xeon E52697 V2 processor occurs nearly for the same value of I a as on the Intel Xeon E78837 processor. Indeed, the critical arithmetic intensities of both processors are approximately the same (respectively 3.4 and 3.6 Flop/Byte), which justify the similarity of the observed trends.

▶ In this chapter, we presented a first parallel implementation of the sweep operation targeting modern SIMD-enabled multicore processors. We first studied the arithmetic intensity of the sweep operation depending on whether the neutron currents are stored or not, for different order of the angular quadrature and for different spatial mesh resolution. We found out that when the neutron currents are not stored, then the arithmetic intensity is dramatically improved which allowed us to efficiently vectorize the sweep operation. While the spatial vectorization is a promising algorithmic strategy, its efficiency is theoretically limited by the required padding to enforce data alignment. Meanwhile, the angular vectorization gives good performances on our test machines. For instance, on the SSE-enabled Intel Xeon Westmere processor, the angular-vectorized implementation of the sweep achieves 13.6 Gflop/s when using a S 16 angular quadrature and MacroCell of size 100 × 100 × 100, corresponding to 62.9% of the theoretical peak performance per core of that processor. In practice, we intend to use vectorized kernels in parallel. We will see in following sections that this additional level of parallelism is applied to other MacroCells that can be processed in parallel. Indeed, the smallest the MacroCells are, the biggest the parallelism potential. Hence, the optimal macrocell size is a trade-off between the kernel performance and the parallel efficiency.

Comforted with this high performance per core, we are going to present in Chapter 3, a theoretical study of the sustainable performance of the sweep operation on clusters of multicore processors. Our goal is to develop a highly efficient parallel neutron transport solver on massively parallel architectures. To achieve this goal, we rely on the discrete ordinates method (S N ) for solving the neutron transport equation. The vast majority of computations, within the S N method, is spent in the sweep operation. Therefore, it is essential to first design a highly efficient parallel sweep algorithm. Considering this requirement, we adopted a bottom-up process by first answering the question on how to maximize the floating point performance of the sweep operation on a single computing core. In Chapter 2, we found that for today and future architectures, to maximize the performance per core, it is essential to make use of the current SIMD vector units. In this chapter, we move to the upper-level to study the parallel performance of the sweep operation on distributed multicore-based architectures.

A substantial amount of work on the neutron transport theory, for the last decade, was focused on the development of sweep algorithms that are capable of maintaining a high efficiency on large number of cores. In section 3.1, we will introduce basic notations that are going to be used to develop the theoretical performance models. Then, we present some previous work on the modelization of the sweep operation (section 3.2) before presenting the two major contributions we made on this performance modelization (section 3.3 and section 3. A comparative study of the theoretical models is given in section 3.5, where all performance measurements were carried-out using vector instructions corresponding to target architectures.

Preliminary definitions

We consider a 3D spatial domain D, discretized into N x , N y , N z cells along the x, y and z axis:

D = 1, N x × 1, N y × 1, N z ,
and a 3D angular discretization comprising N dir = 8 × M directions, where M is the number of directions per octant. The sweep problem is defined as following: Problem 1 (Parallel Sweep). Find the minimal execution time of a parallel implementation of the sweep, over a spatial domain D, along N dir angular directions, on distributed multicore supercomputers.

The vast majority of computations performed during the resolution of the neutron transport problem, according to the S N method, is concentrated in the sweep operation. Thereby, for the last decades, a substantial part of the researches on neutron transport simulations was dedicated to the development of scalable and efficient sweep algorithms. Successive works proposed different approaches to overcome the sweep problem, and focused mainly on the development of new strategies aiming to maintain high parallel efficiency on large number of cores. However, as presented in section 1.4, the architecture of today and future parallel computers has shifted from traditional machines. Thus, algorithms and traditional parallel programming models should be adapted accordingly. From a performance point of view, the efficiency metric must be used in addition with a complementary metric: the "speed of execution" of a given kernel, or Flop/s. Indeed, this latter metric characterizes how much an application uses the power of a given computer (see section 1.4.2).

Let us first give some definitions and present notations that are going to be used to describe the theoretical models. For convenience, all the notations are reported in Table 3.1, where the following variables are used:

• A x , A y and A z are the number of cells along the x, y, z dimensions of a single MacroCell (see Definition 1);

• S x = N x /A x , S y = N y /A y and S z = N z /A z , the global number of MacroCells along each axis of the spatial domain;

• P , Q and R, the number of processes along x, y and z axis;

• L x = S x /P , L y = S y /Q, and L z = S z /R the local number of MacroCells on each process, along each axis of the spatial domain.

Definition 4 (Task). A task represents the computation of all neutron fluxes on a single

MacroCell, of Cartesian coordinates (a, b, c), for a subset of angular directions, A m (belonging to a single octant o), and for a subset of energy groups A g . We denote this task as T o,a,b,c,g .

In the remaining of this chapter, we consider only a sweep for one group. Thus, for the sake of clarity a task will be simply denoted by T o,a,b,c .

Definition 5 (Front). A front of order f , for a given octant o denoted U o

f , is the set defined as following:

U o f = { T o,a,b,c | a + b + c = f } .
∀o, F = S x + S y + S z -2 is the number of diagonal planes (front) of MacroCells per octant, on the whole domain. This allows us to define the list of tasks belonging to the front f, for all the octants,

U f = 8 ∪ o=1 U o f ,
and the total number of tasks available in the sweep problem,

U = F ∪ f =1 U f . Proposition 1. Assume that S x ≥ S y ≥ S z . Then: ∀o, ∀f, U o f ≤ S y × S z .
Proof. The map defined as:

F : U o f → 1, S y × 1, S z T o,a,b,c → ( b c
) , is an injection. Therefore, it follows that the number of elements in the set U o f is less than that of 1, S y × 1, S z ■

In the remaining of this chapter we will assume that S x ≥ S y ≥ S z .

Literature review on performance models of the sweep

The BTE resolution represents a significant portion in the main applications targeted by the DOE's Advanced Simulation Computing1 , formerly known as ASCI [START_REF] Alexander R Larzelere | Creating Simulation Capabilities[END_REF]. This initiative led to numerous research activities on the resolution of BTE. An important part of these researches concerns the development of efficient parallel sweep algorithms, and can be categorized in two different classes: Flat and Hybrid models. The first class corresponds to models implemented by assuming a uniform view of a parallel computer as a collection of independent computing cores. Flat models are usually implemented by using the standardized message-passing interface MPI [START_REF] Gropp | A High-Performance, Portable Implementation of the MPI Message Passing Interface Standard[END_REF]. The second class integrates the hierarchy of modern distributed memory clusters as described in section 1.4.1, and are implemented using several parallel programming models together (message-passing and threading for instance). This section reviews some previous works on performance modelization of the sweep operation using Flat and Hybrid models.

Flat models

In [START_REF] Baker | An SN Algorithm for the Massively Parallel CM200 Computer[END_REF], Koch, Baker, and Alcouffe have proposed the KBA algorithm which decomposes the 3D Cartesian grid onto a 2D process grid. This reference algorithm, in the field of S N Cartesian neutron transport simulation, is also used in the neutron transport code Denovo [START_REF] Gregory G Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF], developed at ORNL 2 . Moreover, the codes Pentran [START_REF] Courau | 3D Neutron Transport and HPC: A PWR Full Core Calculation Using PENTRAN SN Code and IBM BLUEGENE/P Computers[END_REF] and unic [START_REF] Kaushik | Enabling high-fidelity neutron transport simulations on petascale architectures[END_REF] partition the global problem onto a 3D virtual grid of S × A × G processes, where S, A, and G represent respectively the number of processes allocated for the spatial, angular and energy decompositions. The KBA algorithm served as a basis for a large majority of parallel sweep algorithms on structured meshes.

Structured meshes

For structured meshes, Hoisie et al. [START_REF] Hoisie | Performance and Scalability Analysis of Teraflop-Scale Parallel Architectures using Multidimensional Wavefront Applications[END_REF] extended the KBA algorithm parallel performance model, by taking into account both communications and computations that fit the SWEEP3D1 MPI based application. Chaussumier in [START_REF] Silber-Chaussumier | Recouvrement des communications et des calculs, du matériel au logiciel[END_REF] studied the impact of software pipelining in the overlap of communications by computations. In [START_REF] James | Comparison via Parallel Performance Models of Angular and Spatial Domain Decompositions for Solving Neutral Particle Transport Problems[END_REF], Azmy et al. give a method for finding the best way to decompose a given fixed-size problem, between angular and spatial decompositions, through performance models.

Recently in [1], Adams et al. proposed a generalization of the KBA algorithm on structured meshes. As opposite to the KBA algorithm which parallelizes the sweep operation over planes, Adams et al. model defines a volumetric decomposition of the spatial domain according to the x, y and z axis. Indeed, they introduced a 3D spatial domain decomposition onto a 3D process grid. They have also presented scheduling algorithms, for the sweep operation, proved optimal. In the following, we will refer to this model as Adams et al. model. Its general description is given as following. Let G be the number of energy groups; A m and A g the numbers of angular directions and energy groups per task; A z the number of z-planes each process needs to compute before a communication step occurs. N k = L z /A z gives the number of communication steps per process. Hence, each task carries the computation of angular flux for A m directions, A g energy groups and L x × L y × A z cells. Note that the number of cells per task, according to the x and y dimensions, can be modified through aggregation factors A x and A y . According to Adams et al. model, the makespan of the the whole sweep is therefore:

T Adams = N stages (T task + T comm ), (3.1)
where T task is the cost of a computation step; T comm the time needed by a process to communicate its outgoing angular flux to its neighboring processes after each computation step; and N stages the total number of computation steps required to perform the whole sweep:

N stages = 2 (⌈ P 2 ⌉ -1 + ⌈ Q 2 ⌉ -1 + N k (⌈ R 2 ⌉ - 1
))

N fill + N p tasks 8M GN k /(A m A g ) . (3.2)
N fill is the minimum number of stages before a sweep front can reach the center-most processors, and N p tasks is the number of tasks per process. A key point to note here is that this Adams et al. model is able to predict the multigroup sweep computation time rather than just for a single-group.

Unstructured meshes

Flat performance models for the sweep operation, on unstructured meshes, have been largely studied. Thus, Kerbyson et al. in [START_REF] Kerbyson | Performance Modeling of Deterministic Transport Computations[END_REF] presented the first performance model for transport sweeps on unstructured meshes. A similar work has been carried out by Plimpton et al. as described in [START_REF] Steven J Plimpton | Parallel Sn Sweeps on Unstructured Grids: Algorithms for Prioritization, Grid Partitioning, and Cycle Detection[END_REF]. In [4], Kumar et al. explored scheduling strategies for a generalized sweep algorithm. Their algorithm is applicable to more general cases than radiation transport problems as it uses no geometric information about the mesh. In [START_REF] Colomer | Parallel Algorithms for SN Transport Sweeps on Unstructured Meshes[END_REF], the authors presented new algorithms for the sweep operation on unstructured meshes by designing algorithms that achieve overlap of communications by computations, as well as message buffering to reduce cost associated with the latency of parallel machines used.

Recent advances on parallel computer architectures motivated the research community to develop new strategies to leverage the full power of these modern architectures. This evolution trend, in the S N neutron transport community, has led to the development of Hybrid sweep algorithms.

Hybrid models

In [START_REF] Yan | Optimizing Parallel Sn Sweeps on Unstructured Grids for Multi-core Clusters[END_REF], Yan et al. proposed a performance model of the sweep operation for unstructured meshes, targeting multicore-based clusters. Their model favors the overlap of communications by computations by explicitly dedicating a single thread (master) per process to handle all communications performed locally by the parent process, while the remaining threads (workers) keep doing computations. To develop this model, they first represent the sweep algorithm as a Direct Acyclic Graph (DAG), before partitioning it into sets of vertices. Finally, these sets are evenly mapped on processes. According to this model, the sweep execution time is decided by the maximum running time of the master and workers. Their communication model is proportional to the number of processes, and does not integrates a dynamic overlap of communications by computations according to the problem size.

Except for the recent release of Partisn1 , as presented in [START_REF] Baker | PARTISN on Advanced/Heterogeneous Processing Systems[END_REF], which features a parallel implementation of the sweep using both MPI and OpenMP; and the Yan et al. sweep model, all the previously listed codes and theoretical models (see section 3.2.1) follow the classical Flat parallel programming model. Even if modern MPI implementations can directly take advantage of shared memory to synchronize two processes on the same node [START_REF] Richard L Graham | Open MPI: A high-performance, heterogeneous MPI[END_REF], it remains that some extra memory is used by MPI for managing communications [START_REF] Buntinas | Cache-Efficient, Intranode, Large-Message MPI Communication with MPICH2-Nemesis[END_REF] due to the double-buffering strategy. Moreover, as mentioned in section 1.4, each node of modern distributed multicore machines comprises an increasingly large number of cores. In addition, the main memory capacity per core is dropping by a factor of 30% every two years (see [START_REF] Mutlu | Memory Scaling: A Systems Architecture Perspective[END_REF]), and on-chip memory bandwith is not increasing as fast as the available computational power. Thus, it is essential to employ a programming model capable of taking into account these constraints. From this perspective, with a Hybrid model, that combines the use of MPI to manage inter-node communications, and threads within a multicore node to synchronize local data, we can benefit of lower memory latency and data movement on each node, as synchronizations can be done via the shared memory without extra copies. Bearing this idea in mind, we introduce in the next section a new performance model targeting distributed multicore machines. We will refer to this model as Hybrid model.

New performance model of a parallel sweep

Considering the need to shift the parallel programming models from Flat to Hybrid, we aim at designing and implementing a new hybrid sweep algorithm for modern architectures. The algorithm we propose is given as following: loop over all fronts, and perform the following actions on each front f . 1. On each node, a single multithreaded process executes all local tasks that belong to the front f ; each thread being bounded to a unique core. Thus, data-transfers between the tasks are achieved through shared-memory accesses.

2. On each node, a single communication thread is dedicated to process all communications involving that node. Each communication step, occurring on the front f , can be processed while that node executes the tasks of the front f + 1. Such a strategy enhances overlap of communications by computations as will be shown in section 3.3.2.

The Adams et al. model does not apply for such implementations. Thereby, we aim at extending this model to take into consideration the two levels of parallelism available on distributed multicore systems. The goal is twofold:

• validate the need for such an approach;

• provide us insight on the ideal process grid for a given problem and architecture.

To derive this model, we need to evaluate the new number of computation steps N stages , and a new global communication time T * comm . Indeed, in this model not all steps induce a communication step anymore thanks to the shared memory accesses. Furthermore, in this model, the constraint on the coupling between computation and communication steps is released, allowing to minimize explicit global synchronizations. Thus, the makespan of the sweep operation in this Hybrid model is given by:

T Hybrid = N stages T task + T * comm . (3.3) (

Computation steps

A computation step is defined as the execution of a set of tasks in parallel. This set may be reduced to a singleton. To derive the formula giving the total number of computation steps, let us consider the following notations:

• For a given process grid (or domain partitioning) defined by (P, Q, R), N (p, q, r) is the node to which belongs the process of rank r × P × Q + q × P + p. For the Hybrid model, there is a one-to-one mapping from the process list to the node list. Therefore, unless explicitly stated otherwise, we will use interchangeably process and node namings in the description of the Hybrid model

• N = ∪ (p,q,r)∈ 1,P × 1,Q × 1,R N (p, q, r),
is the set of all nodes

• T (p, q, r) is the set of all tasks belonging to the node N (p, q, r): • W f is the set of all nodes having at least one task of the front f , for all octants:

T (p, q, r) = { T o,a,b,c ∈ U | a o L x = p, b o L y = q, c o L z = r } ,
W f = { N (p, q, r) | ∃o, U o f ∩ T (p, q, r) ̸ = ∅ } ;
we will call each of those nodes a working node

• N w f is the number of tasks of the front f , for all octants, on the node w = N (p, q, r) ∈ W f :

N w f = |{U f ∩ T (p, q, r)}|
• N f is the number of computation steps to process all tasks of the front f

• N cores is the number of cores per node These notations are defined in a general 3D case, but in the following, for the sake of clarity, we will rather use 2D figures. Figure 3.1 shows the evolution of the sweep for all the four quadrants of a 2D spatial grid of 6 × 6 MacroCell. The computer used in this example is a cluster comprising nine nodes; each node having two computing cores dedicated to execute the tasks. For this example, there are eight nodes involved to process the front f = 2:

W 2 = {N (0, 0), N (0, 1), N (0, 2), N (1, 0), N (1, 2), N (2, 0), N (2, 1), N (2, 2)} .
The largest number of ready tasks on a single node, for the front f = 2, is max w∈W 2 (N w 2 ) = 2. Therefore, the number of computation steps required to process all tasks belonging to this front is N 2 = ⌈2/2⌉ = 1. In general, for a given working node w, if N w f is a multiple of N cores -1 1 , then the number of computation steps required to process all tasks of the front f on the node w is exactly equal to the quotient of N w f by N cores -1. Otherwise, we could at least give an upper bound on the number of computation steps according to the equation (3.4).

N f = ⌈ max w∈W f (N w f ) N cores -1 ⌉ .
(3.4)

According to the Proposition 1, we have

N f ≤ N y × N z N cores -1 .
Actually, considering the nearest superior integer when evaluating the number of computation steps is somewhat restrictive, as it overestimates this number, by implicitly requiring to finish all tasks of a given front before moving to the next. As an example, consider the sweep of the quadrant 0, over a 6 × 6 regular grid, using 2 parallel threads, as depicted in Figure 3.2. In this case, according to the equation (3.4), the complete sweep requires 21 steps. However, if we allow processing of tasks on U 0 f +1 , even if all tasks on U 0 f have not been yet finished, the same sweep takes 19 steps. This issue is solved by a simulator that will be presented in section 3.4. By summing up the equation (3.4) over all fronts, we obtain an upper bound of the number of computation steps required to process the whole sweep:

N stages = F -1 ∑ f =0 N f .
(3.5) Figure 3.3 shows a verification procedure of this formula, by counting the number of computation steps in a 2D case.

Evaluating the total number of the computation steps required to process the whole sweep algorithm is the first step toward a full modelization of the sweep execution time. Indeed, communications costs should be evaluated according to the process grid partitioning and the network characteristics of the target computer. The section 3.3.2 details strategies we developed to handle communications cost in the Hybrid model of the sweep.

Communication steps

We define a communication step as a data transfer, from a process sender A to a remote process receiver B. The amount of data sent is equal to the size of a MacroCell face: D u , u = x, y, z, times the number M of angular directions per octant, as presented by equation (3.6).

       D x = A y × A z × M D y = A x × A z × M D z = A x × A y × M (3.6)

Communication time on a single front

Let us first assume that every working node, w = N (p, q, r) ∈ W f , has finished processing its set of ready tasks, T (p, q, r)∩U f . These nodes can therefore start processing tasks in T (p, q, r)∩U f +1 , while local communication threads perform data transfer toward neighboring nodes.

Patterns of communications

We assume that all nodes communicate in parallel, as depicted in Figure 3.1b. Let us define:

• T * comp (w, f ) the time required to process all the N w f tasks;

• T * comm (w, f ) the time required by the node w for communications, when all its N w f tasks have been processed.

Then the global communication time for the front f is given by equation (3.7).

T * comm (f ) = max w∈W f T * comm (w, f ). (3.7)
Due to our hypothesis, requiring to finish processing all tasks of the front f before a communication step can occur, T * comm (f ) should be evaluated on a node having executed the highest number of tasks. In the following, this node is denoted by w max :

w max ∈ { w ∈ W f | T * comm (w, f ) = T * comm (f ) and T * comp (w, f ) = max w i ∈W f T * comp (w i , f ) } .
Proposition 2. For a given front, the set of nodes performing the highest number of communications and computations can be considered as a singleton.

Proof. Let w 1 , w 2 be two nodes belonging to this set. Recall that in this Hybrid model, it is required to finish all tasks in U f before moving to the next front. Also we assume that any two different nodes can communicate in parallel, with their respective neighbors; and that there is no network contention on the message delivery. According to these requirements, communication times for the nodes w 1 and w 2 shall not change. Consequently, no matter on which node, from this set, we evaluate the communication time on the front f . ■

Cost of a single data transfer For a given front f , let M w f (u)|u ∈ (x, y, z) be the number of MacroCell faces sent by the node w towards the u direction. Hence:

T * comm (w, f ) = ∑ u∈(x,y,z) M w f (u)τ u , ( 3.8)
where τ u u ∈ (x, y, z) is the time to required to complete the communication of a single face of a MacroCell, along the u direction. Estimating the communication time of a message over a distributed memory computer represents a tremendous research area. The linear communication model (see [START_REF] Hockney | The Communication Challenge for MPP: Intel Paragon and Meiko CS-2[END_REF][START_REF] Schmidt | Empirical Analysis of Overheads in Cluster Environments[END_REF]) is widely used. It is defined by equation (3.9).

τ u = α + D u β , ( 3.9)
where α and β are respectively the latency and bandwidth of the network interconnect, measured between any two nodes. This model is a good approximation of the communication time for two any computing nodes directly interconnected. However, on a typical cluster, the topology of the network interconnect may introduce a significant variation on this model. Indeed, in such an environment, the communication cost may be different depending on the considered pair of nodes as explained in [START_REF] Clematis | Modeling performance of heterogeneous parallel computing systems[END_REF], especially due to the message routing issues and the network contentions. A more elaborated communication model that takes into account concurrency over the network resources on nodes, useful to cope with large messages, is presented in [START_REF] Martinasso | A Contention-Aware Performance Model for HPC-Based Networks: A Case Study of the InfiniBand Network[END_REF][START_REF] Martinasso | Model of concurrent MPI communications over SMP clusters[END_REF].

For our study, we altered the linear communication model as following:

τ u =      o c • D u + T cr (D u ), D u < d 0 o c • D u + α + D u β , D u ≥ d 0 (3.10)
where T cr (d) returns the measured communication time between two nodes for a message of size d; and o c (expressed in seconds per byte) is a fixed communication overhead per transferred byte. We use the network benchmarking utility NetPIPE [START_REF] Quinn O Snell | NetPIPE: A Network Protocol Independent Performance evaluator[END_REF], for evaluating T cr (u). d 0 is experimentally set to the message size from which the messages are transferred with the bandwith of the communication network, so that the linear communication model can be used to predict the communication cost.

Global communication time

The global communication time for the sweep operation is finally given by equation (3.11).

T * comm = F -1 ∑ f =0 (1 -k f ) max w∈W f T * comm (w, f ) (3.11)
A part of them can be overlapped by computation threads. To take this effect into account, it is necessary to introduce an overlap rate, depending on the front being processed k f . Evaluating the overlap rate on a given front f , and on a given node w, requires to determine the amount of tasks on the front f + 1. This overlap rate is evaluated on the node w max as defined in section 3.3.2. If T * comp (w max , f + 1) > T * comm (w max , f ) then all communications on the node w max could be hidden. Otherwise, the amount of communications that could not be overlapped by computations is:

T * comm (w max , f ) -T * comp (w max , f + 1).
Consequently, the overlap rate is:

k f = max ( T * comm (w max , f ) -T * comp (w max , f + 1), 0 ) T * comm (w max , f ) , f = 0, • • • , F -2.
(3.12) Indeed, we do not evaluate k f on the last front because there is no communication at this step.

Asynchronous simulator of the sweep

The Hybrid theoretical performance model presented in the section 3.3, under-exploits the available parallelism in the sweep operation. Indeed, this model requires to execute all tasks on a given front before moving to the next one. Such a procedure imposes an explicit synchronization, which can therefore leads to an under-utilization of the computational resources. To overcome this issue, we propose a new theoretical model having the following properties:

• On each node, a multithreaded process executes local tasks asynchronously

• On each node, a dedicated thread handle all local communications, as for the Hybrid model, hence enhancing the overlap of communications by computations

• A task can be processed as soon as it has received all its dependencies Such a model allows to leverage more parallelism from an algorithm, and thus maximizes the occupation of computational resources. To study the efficiency and the performance of this parallelization strategy, we need to evaluate the number of computation steps, as well as the communication costs, by integrating the ability for a node to execute several tasks belonging to successive fronts. As one may observe, this model can not be easily described using a closed-form expression as we did for the Hybrid model. There are some existing tools aiming to simulate the execution of an application on distributed multicore systems. For instance, SimGrid, presented in [START_REF] Casanova | Versatile, Scalable, and Accurate Simulation of Distributed Applications and Platforms[END_REF], allows to have an abstract model of a large-scale distributed system, so that the execution of a parallel application can be simulated on that system. However, such tools are more generic and versatile, and thus do not explicitly take into account specific particularities of a given algorithm, which can afford to simplify the simulation. Thereby, to simulate the execution of Hybrid sweep on distributed multicore systems, we rather designed a simulator that mimics the behavior of an asynchronous scheduler, that keeps executing tasks at earliest time, and tailored for the sweep algorithm. Such a simulator is exploited in modern task-based runtime systems (see Chapter 4). This section presents the design of this simulator, referred as Hybrid-Async in the following. Without loss of generality, the simulation algorithm is presented in a 2D case.

General presentation of the simulation algorithm

The theoretical performance model presented here predicts the computation time of the sweep operation on a distributed memory computer. This is achieved by simulating the execution of an asynchronous implementation of the sweep operation on such a computer, according to the problem size and the characteristics of computational resources used (number of nodes; number of cores per node; latency and bandwidth of the network interconnect). Such a strategy has the potentiality of better estimating the computation time, compared to a formula-based model as presented in section 3.3.

Preliminary considerations

Let us consider:

• a discretization of the time variable into intervals of length equal to the computation time of a single task, T task :

     I l = [l • T task , (l + 1) • T task [ , l ≥ 0 [0, +∞[= ∪ l≥0 I l (3.13)
• a global time, clock, common to all nodes, and set to 0 at the beginning of the simulation;

• a local time for the communication thread, t (p,q) c , giving the date at which the communication thread on the node N (p, q) is available for performing next communications.

A task T q,a,b is represented by equation (3.14

) { T q,a,b := {start, end, dep, T task } T q,0,0 .start = 0.0 , ( 3.14)
where start and end represent the beginning and finishing execution dates of this task; dep is the number of dependencies1 for this task. One should note here that this modelization corresponds to a simultaneous sweep of all quadrants. In fact, the beginning execution dates for the the four corner tasks (T q,0,0 , q = 0, 1, 2, 3), are set to 0.0. Consequently these tasks could be executed at the beginning of the simulation.

Definition 6 (Task scheduling and execution). A task T q,a,b is scheduled if all its dependencies have been removed: T q,a,b .dep = 0.

If clock ≥ T q,a,b .start, then we say that the task T q,a,b is ready to be executed. The execution of a ready task T q,a,b allows to set its finishing execution date according to the equation (3.15):

T q,a,b .end = clock + T task .

(3.15)

The simulation algorithm

The simulation engine encapsulates the behavior of an asynchronous implementation of the sweep operation, over a multicore distributed machine: local tasks are executed by a set of working threads (we assume one thread per computing core), and communications are carried out sequentially by a separate thread per node. Proposition 3 defines the main idea behind the simulation algorithm. Proposition 3. On each time spacing I l , a given node could execute a maximum of N cores -1 tasks.

Proof. Let w be a node with a set of n tasks ready to be executed at time clock, and assume that n > N cores -1. On this node, there are N cores parallel threads, including a single one dedicated to communications. By assuming that the computation time of tasks is uniform, and equal to T task , it follows therefore that N cores -1 tasks, among the n ready tasks, are going to be executed by clock + T task / ∈ I l . The remaining tasks could not be executed since the global time is out of I l . ■

The main idea of the simulation algorithm (see Algorithm 7) is then to keep moving forward in time, and to perform the following action on each node: execute ready tasks on each node (Line 5) and update dependencies of adjacent ones; realize data transfers and schedule ready tasks for execution (Line 7 and Line 8). This process continues until all tasks are executed. The value of clock at the end of the simulation gives the sweep time.

Algorithm 7: General simulation algorithm of an asynchronous sweep in 2D In :

• Latency (α)

• Bandwidth (β)

• Scheduling overhead (δ)

• Communication overhead (o c )
Out: Sweep execution time (clock)

1 clock = 0.0; ▷ Global clock 2 RemainingTasks = 4 × S x × S y ; 3 while RemainingTasks > 0 do 4 for w = N (p, q) ∈ N do 5 ▷ Execute n ready tasks on N (p, q): T q,a,b 1 , • • • , T q,a,b n (n ≤ N cores -1) 6
RemainingTasks -= n; Scheduling and execution of tasks We recall that, as the sweep goes through the grid from a given corner to the far opposite corner, each processed task removes one dependency to one of its adjacent tasks. When all dependencies of a given task have been removed, then this task is scheduled for execution. It is worth to note here that for a sequential execution, only a single task can be released for execution per time step. In this case, the makespan of the problem is simply the total number of tasks times the computational cost of a single task. Opposingly, for a parallel execution, we can have more ready tasks to be executed than available computational resources. It therefore raises the following problem:

Problem 2 (Scheduling). Let an algorithm described by a task graph and a set of computational resources. Given a set of ready tasks at some execution time, find the most efficient way to schedule those tasks in order to minimize the makespan.

This scheduling problem is a more general concept and it is subject to numerous studies in the parallel computing area. For the particular case of the Hybrid-Async model, the available parallelism during the sweep operation is dynamically discovered as the sweep moves throughout the spatial grid. We state the following requirements to achieve an optimal scheduling: 1. The execution of a set of ready tasks should contribute to release the highest number of tasks that can be released at this step;

2. Locally executed tasks should contribute to minimize idle time of neighboring nodes.

In the Hybrid-Async model, tasks are sorted in respect to these requirements. However, for a real implementation implementation on top of a generic runtime system (see section 4.3) for instance, such scheduling policies may introduce an overhead per task, depending on the heuristics used internally in the runtime in order to decide which tasks should be scheduled for next execution.

We integrate this property in the simulator as following: once a given task T q,a,b has received all its incoming dependencies from tasks T q,a,b-1 and T q,a-1,b (see Figure 3.4), at a given time t r , then it can be scheduled at the time t s = t r + δ, where δ is the average scheduling overhead per task. This scheduling overhead can be interpreted as the delay needed by the runtime system to select the task for execution according to a given scheduling policy. Hence, the beginning execution date of the task T q,a,b is updated using the equation (3.16):

T q,a,b .start = max ( T q,a-1,b .end, T q,a,b-1 .end ) + δ. ( 3.16)
It is important to note here that this scheduling overhead per task induces some requirements on the task granularity: a task should be coarse enough such that the value of the scheduling overhead can be negligible compared to the task computation time. Otherwise, the performance of the implementation can be very poor. By the same time, the granularity of the task should allow to extract a large number of tasks, from the considered problem, to maximize occupation of computational resources. We will discuss the impact of this scheduling overhead on the performance in section 3.5.

Communication costs

When a given node N (p, q) has processed some of its ready tasks, the global time will be clock + T task . At this time, a set of communications to be processed on this node is notified to the local communication thread. Because there may be some pending communications when this notification is sent, and due to the fact that local communications are handled sequentially, the communication thread can only start processing those communications at time:

t (p,q) c = max ( clock + T task , t (p,q) c ) .
To show how each communication is handled inside the Hybrid-Async model, let us consider a task T 0,a,b ∈ T (p, q) which releases one dependency to its neighboring task T 0,a+1,b , along the x dimension. A communication must occurs therefore from the node N (p, q) to the node N (p + 1, q). As said previously, a dedicated thread is going to handle sequentially all communications initiated on each node. For this reason, the considered communication on the node N (p, q) will be issued at time: max

( t (p,q) c , T 0,a,b .end ) ,
and t

(p,q) c is updated to this time, to handle the sequentiality of communications. Finally, the cost of a message transmission between any two nodes follows the same communication model as we used for the Hybrid model (see section 3.3.2). Again, the communication thread is updated to take into account transmission of the message:

t (lp,lq) c + = τ x .
Algorithm 8 presents the full Hybrid-Async model.

In this section, we presented two different but related algorithms of parallelizing a sweep, on distributed multicore systems: the Hybrid model and the Hybrid-Async simulator. In the next section, we present some performance studies to highlight the advantage of these two algorithms over classical Flat approaches.

Comparative study of the performance models

The goal of this section is to compare all performance models described in previous sections. We first explain how parameters used in the performance models are obtained, then we discuss the efficiency and performance of each model. These studies are realized using two different 3D test cases, and model parameters are set for the ivanoe platform (see Table A.1).

Parameters of the models

This section specifies measurements of the parameters used to benchmark the theoretical performance models.

Computer network performances For a given machine, we use the NetPIPE [START_REF] Quinn O Snell | NetPIPE: A Network Protocol Independent Performance evaluator[END_REF] utility to evaluate communication time for message of increasing sizes. This allowed us to evaluate the network performance of the ivanoe platform, and the result is presented in Figure 3.5.

Task computing time

In the neutron transport community, the computation time per spatial cell, per angular direction and per energy group is usually called grind time, denoted as T grind in the following. Thus, the computation time of a single task can therefore be evaluated according to the number of spatial cells, angle and energy groups aggregated into it. This linear model of computation time per task ignores cache effects and bus memory contentions due to concurrent accesses; but it gives at least a lower limit on the task computation time.

In this study, the grind time is evaluated through experimental measurements. For each test case, we set this value to the average computation time per spatial cell and per angular direction, using a sequential implementation of the sweep operation. It is worth to note here that in practice, the grind time is not a linear function of the number of cells. To assess this hypothesis, we conduct an experiment consisting to evaluate the average performance per MacroCell, from the performance of a sequential implementation of the sweep operation, on the bigmem computing node (see Table A.1). To achieve this, we considered a test case featuring 5 × 5 × 5 MacroCells of increasing sizes. This study is presented on Figure 3.6. We observe that the performance per MacroCell (or indirectly, per spatial cell) increases, non-linearly, with the number of cells and increases slightly starting from MacroCells of size 60 × 60 × 60 cells, and corresponding to an average performance per MacroCell of 10.6 GFlop/s. As the MacroCell sizes (A x , A y and A z ) define the task granularity in all the performance models, one should be aware on the choice of the grind time. In the following, we use a single grind time per test Algorithm 8: Asynchronous sweep simulator in 2D

In :

• Latency (α)

• Bandwidth (β)

• Scheduling overhead (δ) ▷ Perform communications and schedule ready tasks

•
for T q,a,b ∈ { T q,a,b 1 , • • • , T q,a,b n } do ▷ Communications if T q,a,b .comm x then t (lp,lq) c + = τ x ; T q,a+1,b .start = max ( T q,a+1,b .start, t (lp,lq) c ) ; if T q,a,b .comm y then t (lp,lq) c + = τ y ;
T q,a,b+1 .start = max ( T q,a,b+1 .start, t

(lp,lq) c ) ; ▷ Scheduling if T q,a+1,b .ready then T q,a+1,b .start + = δ if T q,a,b+1
.ready then T q,a,b+1 .start + = δ clock+ = T task ; return clock; case, obtained by averaging the computation time for a sequential run of sweep operation (see Table 3.2).

Communication and scheduling overheads

For a given computer, we choose the best values of these parameters that fit the sweep models to real execution times. It is worth to note that the scheduling overhead (δ) is strongly dependent on the runtime being used. We will use different values of δ to illustrate its impact on the performance when the task granularity is not quite large enough to amortize it.

We will first evaluate the performance models by considering that scheduling overhead per task, and communication overhead, are null.

Evaluation of the performance models

We consider two different test cases: small and big as described in Table 3.2. These test cases differ on their spatial mesh sizes; the total number of angular directions remain the same, and corresponds to a S 16 Level Symmetric quadrature (288 angular directions). The small test case is intended to stress the behavior of theoretical performance models when there is not enough parallelism; while the big test case allows to perform a strong scalability study on large number of cores.

Optimizing process grid and task-granularity

We recall that the main goal of the developed performance models is to optimize the process grid partitioning for a set of computational resources and for a given test case. Thereby, for all the theoretical models, for a given test case and a given number of nodes (or processes), we select the optimal process grid partitioning that minimizes the global sweep time by exploring all the possible configurations. Furthermore, for the Adams et al. model, we select the best value of A z that minimizes the total number of stages (equation (3.2)), by solving the Problem 3 (notations are defined in Table 3.1).

Problem 3 (Optimal sweep). The minimal computational time for the sweep operation for a given test case and a set of computational resources is reached for a given domain partitioning

(P opt , Q opt , R opt ),

and a given aggregation factors (A opt

x , A opt y , A opt z ) such that for all domain partitioning (P, Q, R) and for all MacroCell size (A x , A y , A z ):

T sweep { (P opt , Q opt , R opt , A opt x , A opt y , A opt z ) } ≤ T sweep {(P, Q, R, A x , A y , A z )} .

Adams et al. model

It is worth to note here that, for the Adams et al. model, for a given process grid partitioning, the MacroCell sizes along the x and y dimensions are set respectively to N x /P and N y /Q as the authors did. Indeed, optimizing A x and A y makes more sense when the local subdomain has to be processed by a multithreaded process, as for Hybrid and Hybrid-Async models, in order to adjust the amount of local tasks, for efficiency purposes.

Hybrid model and Hybrid-Async simulator

To determine the optimal MacroCell size for the Hybrid model and the Hybrid-Async simulator, we must have a measure of the scheduling overhead per task. Indeed, assuming that the task scheduling is negligible, it is clear that the best decomposition is therefore the one that will extract a large amount of parallelism from the considered problem. Therefore, the optimal MacroCell will be a single cell, which does not correspond to our experimental studies as we are going to see in Chapter 4. This suggests that there is an optimal lower bound for the MacroCells, which can be determined for instance by an estimatation of the scheduling overhead per task. Without having to define this overhead, we can meanwhile compare the performance predictions from our models (Hybrid model and Hybrid-Async simulator) with that of Adams et al. model, by searching for the optimal MacroCell size to be not less than the optimal task size returned by the Adams et al. model, and defined by optimal aggregation factors (A opt x , A opt y , A opt z ).

Performance comparisons of the models

In this section we propose a comparative study of theoretical models presented in previous sections, regarding performance and scalability using both small and big test cases (see Table 3 Hybrid model Figure 3.7 shows the performance of the sweep operation as predicted by the Hybrid model, described in section 3.3, using the small test case. This experiment was carried out using parameters of the ivanoe platform, and N cores -1 computing cores per node; 1 core is set free for the communication thread. This figure presents the performances of the sweep operation for three different settings, regarding the overlap rate of communications by computations. The first case corresponds to an implementation where there is no overlap of communications by computations, which underestimates the performances. In this case, the predicted performance when using 768 cores is 1.96 Tflop/s, corresponding to a parallel efficiency of 23.8%, relative to the performance on a single node (12 cores). In the second case, where we assume a full overlap of communications by computations, the predicted performances are overestimated. More precisely, the performance and the efficiency become respectively 3.5 Tflop/s and 42.3% representing an improvement of a factor of 1.7. However, it is not always possible in practice to overlap all communications by computations. Indeed, as explained in section 3.3.2, this is only possible when there are enough computational tasks per node to hide the progress of data transfers through the network interconnect. To illustrate this, let us consider the Hybrid model prediction with an adaptive overlap rate for 64 nodes (768 cores), which slightly underestimates the performances because the computation steps is overestimated (see 3.3.1). For that number of nodes, the Hybrid model found that [START_REF] Buntinas | Cache-Efficient, Intranode, Large-Message MPI Communication with MPICH2-Nemesis[END_REF][START_REF] Barrault | Efficient Parallel Resolution of the Simplified Transport Equations in Mixed-Dual Formulation[END_REF]5) is the optimal size of the MacroCell (see Table 3.4), which gives a total of 18432 tasks, and thus only 288 tasks per node to be shared by 11 computation cores (≈ 26 tasks per core). Bearing in mind that these tasks are not ready at the same time, it is therefore not possible to achieve a high overlap rate when the number of cores is large. This is what is shown by the predicted performances from Hybrid model with an adaptive overlap rate. Up to 96 cores (960 tasks per core), all communications are hidden; and beyond that point, the overlap rate decreases until 37.5% at 768 cores. Indeed, for this point, there are only 18 tasks per core.

Note 1.

Considering the results of Hybrid using the small test case on 768 cores, we observe that when using MacroCells of size (5,5,5) there are 216 tasks per process, that is ≈ 20 tasks per core, which is largely enough to keep busy all the computing cores. However, even without taking into account scheduling overhead, the Hybrid model returned [START_REF] Buntinas | Cache-Efficient, Intranode, Large-Message MPI Communication with MPICH2-Nemesis[END_REF][START_REF] Barrault | Efficient Parallel Resolution of the Simplified Transport Equations in Mixed-Dual Formulation[END_REF]5) as the best MacroCell size, which corresponds to 36 MacroCells per process. The reason is that communications are much more dominant at this point, and therefore by decreasing task granularity, the amount of communications is larger.

In the following, all results with the Hybrid model are evaluated using an adaptive overlap rate. 3.8 shows a comparison of the performance as predicted by the Adams et al., Hybrid and Hybrid-Async models on the small test case. On a single node, the performances of all the three models are similar. Up to 96 cores, Hybrid and Hybrid-Async achieve nearly the same performances. Indeed, as explained previously, when the number of tasks per core is relatively high, the Hybrid model can perfectly overlap almost all communications by computations. For the same number of cores, we observe that the performance of Hybrid and Hybrid-Async is 1.4 times faster than the performance of Adams et al. model. This is due to the fact that the Adams et al. model involves much more communications than our models. Beyond 384 cores, the Adams et al. model is faster than the Hybrid because the latter overestimates the number of computation steps, as explained in section 3.3.1. However, in the Hybrid-Async model, those constraints are removed and no explicit synchronization exists anymore. For this reason, the predicted performances from the Hybrid-Async simulator is higher that those of Hybrid and Adams et al. models. At 768 cores, the Partitioning Aggregation N tasks Message size (KB) It is worth to note that the Hybrid model is simpler than Hybrid-Async simulator, mainly because it does not take into account asynchronous execution of tasks and evaluates the computation steps using a closed-form expression. Consequently, it is able to quickly predict the sweep computation according to parameters given as input. Opposingly, the Hybrid-Async simulator is more complex and gives better performances, at a cost of the run time of the simulation which is higher than that of Hybrid model (by a factor of ten in average). We performed the same study on the big test case, which is based on a larger spatial mesh. We observe that the Hybrid model is able to hide almost all communications (Figure 3.9), and the predicted performances is slightly higher than that of the Adams et al. model (Figure 3.10), by a factor of 1.02 at 768 cores, even with the overestimation of the computation steps. The Hybrid-Async exploits better the available parallelism, and is able to achieve a very good scalability. The predicted performances from this model is higher than Hybrid and Adams et al. models, respectively by ×1.13 and ×1.16 at 768 cores. This study is carried out with the small test case, using 64 nodes of the ivanoe platform. We fixed the task granularity and the process grid partitioning to the optimal values obtained from a run without any scheduling overhead (see Table 3.4). We made the following observations: the performance of the sweep operation decreases as the scheduling overhead increases, with a large drop of performance (from 2.1 Tflop/s to 1.6 Tflop/s) when using δ = 57.43 µs. This value coincides with the computation time of a single task for the considered task granularity. This behavior is related to the fact that our simulation algorithm executes tasks by time spacing equal to T task . Process grid partitioning and MacroCell sizes are fixed to the optimal values obtained from the run without any overhead (see Table 3.4). The vertical line corresponds to the computation time of a single task.

Adams et al., Hybrid and Hybrid-Async models Figure

N cores P Q R A x A y A z s x s y s z 12 3 2

Impact of the scheduling overhead in the Hybrid-Async model

▶ In this chapter, we have theoretically studied the performance of the sweep operation on multicore distributed systems, predicted by different classes of models. The classical Flat model, which assumes a parallel computer as a collection of independent computing cores, is widely used for designing parallel algorithms. We reviewed one such model, namely Adams et al. model, providing optimal results for parallel sweeps. Due to the shift on parallel computer architectures, moving from uniprocessor-based clusters to multi/many-core based clusters, traditional approaches are not anymore suited to leverage the full power of today computers. Given this situation, we developed two new performance models that take better into account architectures of hybrid machines: these are Hybrid and Hybrid-Async models. We found the predicted performances from the Hybrid model, acting as an implementation using explicitly a combination of message-passing and threading, is higher than that of the Adams et al. model when the considered test case provides enough parallelism to overlap communications by computations. The Hybrid-Async simulator, which simulates the execution of an asynchronous implementation of the sweep algorithm, removes all explicit synchronizations and thus is able to better exploit the available parallelism. This latter strategy increases the occupation of computational resources, and we found that the predicted performances from this simulator remains higher than predicted performances from Adams et al. and Hybrid models.

These results justify our strategy to design a Hybrid implementation of the sweep operation, targeting distributed multicore-based machines. In Chapter 4, we will present an implementation of the sweep operation, following the Hybrid-Async simulator, and using a generic task-based runtime system. We will see that the predicted performances with Hybrid-Async simulator are closer to the measured performances as compared to the Hybrid model. However, the latter model is nevertheless capable to quickly predict the optimal process grid partitioning, and thus it can be used to parametrize a sweep run in production contexts.

Chapter 4 A Massively Parallel Implementation of the Cartesian Transport Sweep Contents

To design an efficient massively parallel neutron transport solver, we considered the discrete ordinates method for solving the neutron transport equation (see Chapter 1). The vast majority of computations performed with this method are spent in the sweep operation. In Chapter 3, we have theoretically studied the performances of that operation on distributed multicore architectures, through the Hybrid model and the Hybrid-Async simulator. We established that the predicted performances of the Hybrid-Async simulator are higher than that by the classical Flat approaches such as the Adams et al. model. In this chapter, we will present a practical implementation of the Hybrid-Async model, using three different task-based runtime systems: PaRSEC [START_REF] Bosilca | DAGuE: A Generic Distributed DAG Engine for High Performance Computing[END_REF], StarPU [5] and Intel TBB [START_REF] Reinders | Intel Threading Building Blocks[END_REF] frameworks.

The emergence of generic task-based runtime systems

As presented in section 1.4, modern parallel computers are built as an interconnection of several heterogeneous computing nodes, each of which comprises a number of computing devices such as multicore processors, accelerators such as Graphics Processing Units (GPUs) or manycore devices. Each of these devices must be addressed with a specific programming paradigm, or using extensions of computer programming languages. Consequently, such an approach implies mixing several programming paradigms together in the same application, through the well known MPI+X programming model.

The traditional MPI+X model

MPI+X is the most popular programming model for modern parallel computers. It consists in using the Message Passing Interface (MPI)1 standard to handle inter-nodes communications; while computing devices on each node are explicitly addressed via programming models specifically tuned for considered node architectures. For instance, to write a computational kernel targeting a computing node equipped with traditional multicore processors, X can be either OpenMP [START_REF] Dagum | OpenMP: An Industry-Standard API for Shared-Memory Programming[END_REF], Intel TBB [START_REF] Reinders | Intel Threading Building Blocks[END_REF] or PThreads. On the other hand, if the computing node comprises GPUs, X can be NVIDIA CUDA [START_REF] Nickolls | Scalable Parallel Programming with CUDA[END_REF] or OpenCL [START_REF] Stone | OpenCL: A Parallel Programming Standard for Heterogeneous Computing Systems[END_REF].

This MPI+X programming model generally follows a fork-join model, which is no longer relevant for today's parallel computers. Indeed, with this model, computation and communication phases are usually serialized, and thus global synchronizations are implicitly introduced in the execution flow. Thereby, the available parallelism inside an algorithm is underutilized, and computational resources are not efficiently exploited. Indeed, some workers go to idle state while they could be making progress on some other work (e.g. for MPI+OpenMP: master thread performs communications while other threads are waiting). These issues motivated the emergence of task-based models on top of generic runtime systems.

Task-based models

Generic concepts of task-based runtime systems

Regarding modern heterogeneous clusters, many initiatives have emerged in previous years to develop efficient runtime systems. Most of these runtime systems use a task-based paradigm to express concurrency and dependencies by employing a task dependency graph to fully represent the application to be executed. This graph is a directed acyclic graph (DAG) where nodes are computational tasks and edges represent data flows and dependencies. The major interest of this approach is the separation of major concerns arising when designing parallel programs:

• Description of the algorithm • Writing optimized computational kernels • Efficient scheduling of the tasks over the underlying hardware These task-based models offer an elegant and efficient way to express parallelism inside an algorithm, while avoiding the cumbersome hand-coding of all the communication primitives needed to realize data transfers. Such models allow the removal of all unnecessary and global synchronizations between tasks; hence exposing a larger number of parallel tasks.

State-of-the-art on task-based models

The past decade has witnessed the development of several approaches of task-based models on top of generic runtime systems. The main differences between these approaches are related to their representation of the graph of tasks, whether they manage data movements between computational resources, the extent to which they focus on task scheduling, and their capabilities of handling distributed and heterogeneous systems. Runtime systems such as Quark [START_REF] Yarkhan | Dynamic Task Execution on Shared and Distributed Memory Architectures[END_REF], StarPU [5], or StarSS [START_REF] Badia | Parallelizing Dense and Banded Linear Algebra Libraries using SMPSs[END_REF] propose an insert task paradigm where a sequential code submits all computational tasks. In this case, the dependency graph is dynamically discovered at runtime according to how the data is used, which is indicated through keywords such as INPUT, OUTPUT or INOUT. Intel TBB is another framework which also allows the use of a taskbased approach. Within this framework, the task dependency graph is constructed according to high level programming paradigms, or algorithms such as parallel_for, parallel_do or parallel_reduce. Another framework is Charm++ [START_REF] Kalé | CHARM++: A Portable Concurrent Object Oriented System Based on C++[END_REF] which is a parallel variant of the C++ language, and allows the writing of programs where a flow of tasks is applied to each piece of data. In addition, it provides sophisticated load balancing and a large number of communication optimization mechanisms. Intel CnC [START_REF] Burke | The Concurrent Collections Programming Model[END_REF] and PaRSEC [START_REF] Bosilca | DAGuE: A Generic Distributed DAG Engine for High Performance Computing[END_REF] construct an abridged representation of the DAG (with its tasks and their dependencies) with a structure agnostic to algorithmic subtleties, where all intrinsic knowledge about the complexity of the underlying algorithm is extricated, and the only constraints remaining are annotated dependencies between tasks [START_REF] Cosnard | Compact DAG Representation and its Symbolic Scheduling[END_REF]. However, as noted in [START_REF] Tang | Measuring the overhead of Intel C++ Concurrent Collections over Threading Building Blocks for Gauss-Jordan elimination[END_REF], the Intel CnC framework which is built on top of Intel TBB incurs some overhead compared to Intel TBB.

For designing and implementing our sweep operation on top of generic task-based runtime systems, we compare the following approaches: hand-written (Intel TBB), parametrized DAG (PaRSEC) and insert task (StarPU). In the following, we will present the implementation of the sweep on top these frameworks. For the sake of clarity, the following presentation will be given in a 2D case, and we consider the case of vacuum boundary conditions which enables a concurrent sweep of all the quadrants.

Implementation of the sweep algorithm with Intel TBB

To implement the sweep operation with Intel TBB, we rely on the parallel_do primitive [START_REF] Robison | Optimization via reflection on work stealing in TBB[END_REF], that enables a dynamic scheduling of the parallel tasks. This parallel function allows a pool of threads to execute the tasks from a list, which is dynamically updated at runtime. At the beginning of the sweep, the task list is composed of four tasks, one for each quadrant: {T q,0,0 , q = 0, 1, 2, 3}. Some of the running threads process these tasks and update the task list to: {T q,1,0 , q = 0, 1, 2, 3} ∪ {T q,0,1 , q = 0, 1, 2, 3} and so on. The Intel TBB parallel_do primitive we used, dynamically explores, executes and updates the list of ready tasks. The computations performed inside a task consist in updating the outgoing angular fluxes in one quadrant (or octant in 3D), for a single MacroCell. These computations are vectorized over angular directions thanks to the generic C++ library Eigen [48] and was presented in Chapter 2.

Note 2. It is worth noting that each

MacroCell object encapsulates a scalar flux field, common to all quadrants. Therefore, for a concurrent sweep of all quadrants, each MacroCell must possess a mutex, to prevent two threads from processing the same data simultaneously. When a thread starts processing a given task T q,a,b , it acquires the mutex corresponding to the MacroCell(a, b). When that processing is finished, the mutex is released, so that another thread can start working on that MacroCell.

Implementation of the sweep algorithm with PaRSEC

PaRSEC is a framework intended to develop parallel applications on distributed heterogeneous architectures. It features a generic data-flow runtime system, supporting a task-based implementation and targeting distributed hybrid systems. This framework relies on the dynamic scheduling of a directed acyclic graph of the considered algorithm: the nodes represent the computational kernels (tasks), and edges represent data transfers between tasks. Thanks to an algebraic description of the task dependencies, the scheduling is completely asynchronous and fully distributed. Moreover, it takes into account user defined priorities and overlaps communications by computations.

As mentioned above, the programming model exploited in PaRSEC enables the separation of the major concerns in distributed computing: the kernels, the algorithm, and the data distribution. Here, the computational kernel is the same as for the Intel TBB implementation: update of angular fluxes in one quadrant (or octant in 3D) for a single MacroCell.

Task-graph of the sweep operation

To use the PaRSEC framework, the algorithm must be described as a DAG using the symbolic representation specific to PaRSEC in a Job Data Flow (JDF) file. In this file, all tasks of the algorithm are defined by their execution space; their data placement or affinity; their input, output or in-out data-flows; and body. Each data-flow has incoming and outgoing edges connecting them to other tasks of DAG, or directly to memory accesses. The body specifies the computations carried out by the task. This file is later compiled into a C code with a set of functions that will submit the tasks to runtime system, check the dependencies, release the data to the following tasks, and execute the body. The Cartesian sweep operation, either 2D or 3D, by its geometric structure, is a natural and simple candidate for this formalism. A simplified version of the 2D sweep operation, without boundary cases and one single quadrant, is given in the Listing 4.1. Only one task ComputePhi, described by its position on the grid (a, b), composes the algorithm. The execution space specifies that there are as many tasks as cells on the grid of size ncx × ncy. The parallel partitioning argument instructs the runtime to run the task ComputePhi(a,b) on the node where the data mcg(a, b) is located (see section 4.3.2). mcg is the structure describing the data distribution, and it is called data descriptor in the PaRSEC taxonomy. Thus, mcg(a, b) represents here the data structure of the MacroCell of coordinates (a, b). Each task ComputePhi has three in-out data dependencies:

• PSIX and PSIY are aliases that correspond to the neutron current associated to one MacroCell, respectively along the x and y dimensions. These two variables are labeled as read/write (RW) to indicate that the neutron current on the incoming faces are overwritten by those on the outgoing faces, that is the neutron currents are not stored. PSIX is the current along the x dimension. It comes from the previous task on this dimension ComputeP hi(a -1, b), and is forwarded to the next task ComputeP hi(a + 1, b). PSIY is the current on the second dimension, and its flow is similar to that of PSIX. On the initial border of the domain, they are directly read from the main memory, and on the final border of the domain, they are written to the initial storage space. In the case of non reflecting nor periodic boundaries, they are directly initialized through a set of initial tasks at the beginning, and destroyed at the end.

• MCG is the alias on the address of the MacroCell object associated to the considered task.

It is directly read from the main memory using the data descriptor mcg. This alias is used to discover some physical data encapsulated inside the MacroCell structure and needed to execute the computational kernel.

The priority line allows the developer to provide a hint to the scheduler helping it to prioritize the most important tasks. We use this feature to optimize the scheduling of the sweep operation (see section 4.3.3). Finally, The BODY section contains the computational task itself exploiting the parameters and flow aliases of the task to access the data.

Data distribution

Once the algorithm has been described in the PaRSEC language, the runtime needs to know how the data is distributed. A simple API must be implemented to provide this information to the scheduler. This is shown in Listing 4.2 for the case of a 2D block distribution of ncx × ncy spatial grid over a P × Q grid of processes, as shown in Figure 3.1a. The rank_of() function is used by the scheduler to know to which process data the belongs to, but also, in our case, for task mapping over the node as previously shown. Two tasks on different nodes can then be detected by the locality of the data they use, and communications are automatically generated: through direct accesses in the case of shared memory, or MPI asynchronous communications on distributed memory. This separation between algorithm and data, specific to PaRSEC, allows a quick evaluation of several data distributions for the same algorithm. Note that it is also possible to have data distribution depending upon subgroup of threads, useful for minimizing the data traffic on NUMA architectures. But for the experiments that will follow, we did not use this feature1 . The data_of() function returns the pointer to the MacroCell object when this one is local. Addresses of transferred objects are internally handled by the runtime.

Optimization of the scheduling through priorities

The PaRSEC runtime implements several schedulers, which may impact the order of task processing and the performance. The default scheduler, Local Flat Queue (LFQ), favors the memory affinity by maximizing memory reuse, by following the dependencies of tasks as defined in the JDF file (Listing 4.1). Tasks released by a dependency are added to the local queue of the working thread. In the case of the 3D sweep, it leads to a prioritization of the sweep per columns of cells (Figure 4.1a). Such a scheduling is not an efficient policy for the sweep since it does not try to maximize the wave front, and thus the number of parallel tasks available. To tackle this problem, we considered the Priority Based Queue (PBQ) scheduler. This scheduler, similar to LFQ, adds ordering of the tasks in the local queue based on optional user defined priorities. We have studied two different priorities: PlaneZ and Front. The first favors tasks belonging to the same z-plane (Figure 4.1b), while the latter favors a progression by front (Figure 4.1c). The first one gives good results when the number of process over z, R, is equal to 2, because it reduces the idle time of processors in the same (x, y) plane; while the second is more generic as it frees tasks quickly along the 3 dimensions.

Implementation of the sweep algorithm with StarPU

Basic concepts

There are two different ways to implement an algorithm on top of StarPU: either using C extensions (pragmas) to annotate a sequential code, or directly using the StarPU's API. In the following, we will focus on the latter. Thus, the implementation of an algorithm with StarPU is relatively straightforward and consists of the following two steps:

1. Create a codelet for each computational kernel in the considered algorithm. A codelet is a structure that holds various implementations of the kernel that it represents, and buffers manipulated by that kernel. 2. Create and submit sequentially all tasks to the runtime system; the runtime can then infer the whole task-graph corresponding to the algorithm. Each task is an association of a codelet and a set of parameters manipulated by the kernel corresponding to that task.

For implementing the sweep operation with StarPU, we consider a single codelet, which corresponds to the update of the angular fluxes on a single MacroCell, for one octant. The MacroCell object associated with the codelet, and the process to which it belongs, is discovered using the same data descriptor we used in the PaRSEC implementation: mcg. In 2D, this codelet defines two read/write buffers, each of which being attached to a single MacroCell face. These are similar to the PSIX and PSIY aliases we used in the PaRSEC implementation (see Listing 4.1). Each buffer is registered to the StarPU runtime with a data handle, which is an opaque pointer that designates an array. As in the PaRSEC implementation, we declare that each of these buffers is allocated on-the-fly at runtime. This codelet defines one computational kernel, which is the same as in the PaRSEC implementation.

Submission of tasks

Let us consider the basic implementation of the sweep operation with StarPU, in 3D, as described in Listing 4.3. This implementation submits all tasks of the first octant, before submitting those of the next one, and so on. Moreover, the submission of tasks for a single octant favors the tasks in the same column (z dimension). However, this implementation has two issues. Firstly, if the number of tasks per octant is large and exceeds the size of the window of visible tasks defined by StarPU, then the octants are going to be processed sequentially, because the StarPU runtime will not be able to visualize all the available parallelism. Thereby, the discovery of the available parallelism is limited, and the computational resources are underexploited. Secondly, due to the fact that StarPU considers the tasks according to their submission order, the considered order of tasks submission (column-wise) does not guarantee that the execution will maximize the front exploration which is the optimal evolution of the sweep as mentioned in section 4.3.3. Thus to enforce the discovery of the available parallelism and to enforce the front exploration, it is essential to submit (or prioritize) the execution of tasks by following the sweep front. To this end, we loop over all fronts and we successively submit tasks belonging to each front for all octants. The Listing 4.4 illustrates this strategy.

Experiments

This section presents the performances of the PaRSEC, StarPU and Intel TBB implementations of the sweep operation, on shared and distributed multicore systems. Performance measurements were carried out in single precision, on the ivanoe and athos clusters (see Table A.1), using the 3D test cases small and big (see Table 3.2). All the experiments were performed with a concurrent sweep over all octants.

To evaluate the efficiency of the task-based implementations of the sweep operation, we will rely on the theoretical performance models presented in Chapter 3. From this perspective, it is necessary to first calibrate those models for the target architectures.

( void *) , STARPU_VALUE , & fx , sizeof ( int ) , STARPU_VALUE , & fy , sizeof ( int ) , STARPU_VALUE , & fz , sizeof ( int ) , STARPU_VALUE , & lnx , sizeof ( int ) , STARPU_VALUE , & lny , sizeof ( int ) , STARPU_VALUE , & lnz , sizeof ( int )
, STARPU_RW , mcg -> getPsi ( mcg , 0 , fx , fy , fz , 0 , y , z ) , STARPU_RW , mcg -> getPsi ( mcg , 1 , fx , fy , fz , x , 0 , z ) , STARPU_RW , mcg -> getPsi ( mcg , 2 , fx , fy , fz , x , y , 0 ) , STARPU_EXECUTE_ON_NODE , mcg -> rank_of ( mcg , x , y , z ) , 0 );

} } } } } }

Task-granularity selection and parameters of the performance models

Task-granularity

To determine the best task-granularity using the Hybrid model or the Hybrid-Async simulator, it is necessary to have a good estimation of the scheduling overhead associated with the runtime in use. Moreover, it is also required to model the cache effects of the target processors, because the performance of the sweep kernel depends on whether or not the MacroCells fit into the processor cache memory. Such a strategy will introduce some low level details into the performance models. Recalling that our goal is to design a simple performance model capable of returning optimal domain partitioning, as a function of the problem size and target machine parameters, we choose to rather adopt an experimental process for evaluating the optimal task-granularity.

Figure 4.2 presents the single-core performances per MacroCell, of the PaRSEC, StarPU and Intel TBB implementations of the sweep operation. There are 5 MacroCells along each dimension, whose size varies from 5 × 5 × 5 to 100 × 100 × 100 (we consider a cubic MacroCell). This experiment was conducted on the athos platform. We observe that for a MacroCell of size 5 × 5 × 5, the performance of the Intel TBB implementation is better than that of PaRSEC and StarPU, and reaches a plateau of 16.8 Gflop/s, starting from MacroCell of size 10×10×10. However, the performance of the PaRSEC implementation increases quickly and stabilizes at 18.1 Gflop/s starting from a MacroCell of size 20 × 20 × 20. Surprisingly, the performance of the StarPU implementation is lower than that of PaRSEC, up to MacroCell of size 50 × 50 × 50. As the computational kernel used in both implementations is the same, this result shows that the scheduling overhead per task for StarPU and Intel TBB is larger than that of PaRSEC, and thus StarPU is suited to schedule coarser tasks. In the following of this study, for all the three implementations, we will use MacroCells of size 20 × 20 × 20 for the big test case, and 10 × 10 × 10 for the small test case.

Parameters of the performance models

On a given machine, we recall that the input parameters to be used for the performance models are the following:

• The average computation time per spatial cell, per angular direction and per energy group, obtained from a sequential run of the sweep on the target machine: T grind

• The scheduling overhead per task, associated with the considered runtime: δ

• The network performance of the target machine and the communication overhead (o c )

The evaluation of T grind has been presented in section 3.5.1. Here, we recall its value for the ivanoe machine and we give its value for the athos machine (see Table 4

Shared memory performances

We have conducted a comparative study of the performances obtained with the three implementations of the sweep, using respectively Intel TBB, PaRSEC and StarPU. This study was realized using a single computing node of the athos machine, and the result is presented on Figure 4.4. The following comparison will be on 23 cores. Indeed, StarPU allocates one thread per process for its internal runtime management which causes a performance drop when using all the 24 cores of the computing node. The three implementations of the sweep were run with the same task-granularity: 20×20×20 cells per MacroCell. We evaluated the Hybrid-Async simulator, without any scheduling overhead, and using T grind obtained from a sequential run of the PaRSEC implementation. The maximum performance predicted by the Hybrid-Async simulator is 404.2 Gflop/s at 23 cores. This performance gives an upper-bound on the performance that can be obtained from an implementation of the sweep operation, on the considered computing node. All three implementations of the sweep achieve nearly the same performances for a run on a single core. This is unsurprising because only a single task is processed; and consequently neither the scheduling policy nor the data locality can impact the performances. However, we found that the performance of the PaRSEC implementation of the sweep, for a run using 23 cores, reaches 367.6 Gflop/s. This corresponds to 90.8% of the predicted performance by the Hybrid-Async simulator, and to 35.45% of the theoretical peak performance of the considered computing node. Moreover, the parallel efficiency of the PaRSEC implementation is 90.9% using 23 cores, and it is respectively 6.76% and 13.06% faster than the implementations with Intel TBB and StarPU. The parallel pattern being constant, we interpret this improved speed-up as a sign for a reduced scheduling overhead for the PaRSEC framework. Indeed, with the PaRSEC framework, the graph of task is not unfolded in memory, thanks to the parametric representation of the DAG: only ready tasks to be executed exist in the system (see [START_REF] Bosilca | DAGuE: A Generic Distributed DAG Engine for High Performance Computing[END_REF]). In addition, the PaRSEC implementation does not have any additional cost associated with insertion of tasks as in StarPU.

The results presented in this section show that the PaRSEC runtime is able to give good performances per core, even when the task granularity is small, as compared to StarPU. A small MacroCell size allows to extract a larger number of tasks from the sweep operation, which is a necessary condition for maintaining a good strong-scalability on large number of cores. For this reason, we are going to focus on the PaRSEC framework for distributed memory performance studies.

Distributed memory performances

On a distributed memory machine, the performance of the sweep operation depends strongly on the domain partitioning. Using Hybrid and Hybrid-Async performance models, we want to determine the optimal process grid partitioning, according to the problem size and the characteristics of the target computer. In this section, we present the performances of the PaRSEC implementation of the sweep operation on distributed multicore architectures, using optimal 3.2), on a single 24-cores computing node of the athos supercomputer. All measurements were obtained using the same task-granularity of 20 × 20 × 20 cells per task. process grid partitioning.

Optimal partitioning of the spatial domain

For each number of cores, we use the Hybrid model and Hybrid-Async simulator to give us the optimal data distribution (P, Q, R) that minimizes the sweep execution time. Figure 4.5 compares the experimental results against the predicted results from Hybrid model and Hybrid-Async simulator, running the big test case with MacroCells of size 20 × 20 × 20, on 48 nodes of the ivanoe platform. We first set the communication overhead to o c = 0 and the scheduling overhead to δ = 0. We found that the predicted performances by the Hybrid-Async follows the same trend as actual measurements from PaRSEC, except for the partitioning (24, 2, 1). Moreover, the PaRSEC implementation and Hybrid-Async simulator gave the same optimal partitioning: (12, 2, 2). Then, we fitted the Hybrid-Async data to the experimental measurements and we found that o c = 1.8 • 10 -10 s/Byte is the value of the communication overhead that minimizes the discrepancy between the Hybrid-Async data and actual measurements. Using this value of o c , the predicted performances comply with actual measurements for each of the partitioning. However, the predicted performances by the Hybrid model, while following the same trend as actual measurements, remain lower than actual measurements. This is unsurprising because, this Hybrid model overestimates the computation step number as explained in section 3.3.1.

The experiment we conducted here has shown that both Hybrid model and Hybrid-Async simulator are able to give optimal data distribution. Thus, all the performance measurements that we will present in the following, are obtained using optimal domain partitioning.

Comparison of Hybrid and Flat approaches

In order to compare Hybrid and Flat approaches, we would have needed a hand-written MPI implementation of the sweep. Nonetheless, we can mimic its behavior with the PaRSEC implementation. To achieve this, we used as many processes as available cores to perform the experiment; each process being bound to one core. In reality, the PaRSEC framework runs an extra thread per process to manage the communications. Thus, for preventing the communication thread from disrupting the computation progress, we dedicated two cores per process: one for the computation, and one for the communication. Therefore, to run the Flat experiments with PaRSEC, the number of processes started on a node is equal to half the number of cores on the node. To ensure a fair comparison of Hybrid and Flat approaches, Hybrid performance measurements were performed by launching one process per node and a number of computational threads equal to half the number of cores. Hence, with both models, the same amount of data is going through the interconnection network.

We performed this experiment on the ivanoe machine, and the result is presented on Figure 4.6. In this particular case, 6 computation cores are being used on each node. Performance measurements using two nodes, that is 12 cores, show that the Hybrid implementation is 30% times faster than the Flat one. To justify this discrepancy, first it must be noted that although intra-node MPI communications can be done via shared memory, it remains that they conduct to a poor usage of caches and saturation of the memory buses. In contrast, these effects are much reduced when using threads in shared memory, and this may justify the observed difference. At 384 cores, the Hybrid implementation is 1.9 times faster than the Flat one. These results confirm that the Hybrid approach is more efficient than the Flat one. This experiment was conducted by running the small test case on the ivanoe platform. By setting δ = 0 and o c = 0, we observe that up to 96 cores (8 nodes partitioned into (2, 2, 2)), the predicted performance from Hybrid-Async simulator complies well with the actual measurements from the PaRSEC implementation. For this number of cores, the Hybrid-Async simulator predicts a performance of 1.01 Tflop/s, while the PaRSEC implementation achieves 0.84 Tflop/s and corresponds to a parallel efficiency of 81.91%. At 768 cores the PaRSEC performances is 1.56 Tflop/s and the corresponding parallel efficiency drops to 18.85%. Moreover, the performance predicted by the Hybrid-Async simulator, for 768 cores, is 2.54 Tflop/s corresponding to a theoretical parallel efficiency of 30.68%, which is 1.6 times more efficient than actual measurements. However, one should be aware that when using 768 cores with the small test case and MacroCells of size 10 × 10 × 10, the average number of tasks per core is 18; and thereby the runtime scheduling overhead is not negligible. To assess this hypothesis, we performed a fit of the experimental data to the simulator data. We found that using δ = 8.1 • 10 -6 s and o c = 4.2 • 10 -10 s/Byte, the predicted performances agrees well with actual measurements. Note that this value of the communication overhead is different from that which was obtained in the big test case (o c = 1.8 • 10 -10 s/Byte). We used this value of o c to run the Hybrid model. We observe that the optimal partitioning is, each time, the same as that of Hybrid-Async, but the predicted performance is lower. Moreover, starting from 384 cores, predicted performance from the Hybrid model is less than that which is obtained from actual measurements. As already discussed previously, this behaviour is related to overestimation of the computation steps in Hybrid model. .8 presents the same study as presented in the previous paragraph, but using the big test case. The computer used is the same ivanoe platform. This test case comprises 110592 tasks, that is 8 times larger than the number of tasks in the small test case, and the average number of tasks per core when using 768 cores is 144. In addition, the larger MacroCell size used for this test case, 20 × 20 × 20, improves the performance per core of the sweep kernel. Thus, on a single computing node, the performance of the PaRSEC implementation reaches 143.9 Gflop/s when using 12 cores, and corresponds to 51.1% of the theoretical peak performance of the node. This high performance per node is explained by the usage of SIMD units and a good data locality exposed by our sweep kernel, thus improving the arithmetic intensity. When using 768 cores, the performance of the sweep implementation with PaRSEC is 6.1 Tflop/s, which corresponds to 33.9% of the theoretical peak performance of the 64 nodes of ivanoe, and the parallel efficiency is 68%.

To assess the PaRSEC performances on this test case, we ran both the Hybrid model and Hybrid-Async simulator using the same value of the communication overhead as used when we studied the optimal partitioning: o c = 1.8 • 10 -10 s/Byte. Unsurprisingly, the predicted performances from the Hybrid-Async model are close to actual measurements and at 768 cores the PaRSEC performances is 88.3% of the predicted value.

▶ In this chapter, we have presented a massively parallel implementation of the sweep operation, using task-based models on top of generic runtime systems: Intel TBB, PaRSEC and StarPU. We compared actual performances obtained from these implementations to the performances predicted by the theoretical models. First, we found that the single-core performance of the StarPU implementation is lower than that of PaRSEC, suggesting that the internal cost of the StarPU runtime is higher than that of PaRSEC. Moreover, on a 24-cores NUMA node based on Intel Ivy-Bridge processors, the performance of the PaRSEC implementation reaches 367.6 Gflop/s, corresponding to 35.4% of the theoretical peak of the node. This performance is higher than those of StarPU and Intel TBB implementations (358.7 Gflop/s and 287.57 Gflop/s respectively).

On distributed multicore machines, we assessed that our Hybrid-Async simulator and Hybrid model are capable to predict the optimal partitioning as a function of the machine parameters and the test case in use. However, the time required to run the Hybrid-Async simulator is higher than that is required to run the Hybrid model because the simulator is more precise, and thus better predict the sweep computation time. Therefore, we will use Hybrid model for evaluating the optimal partitioning and the Hybrid-Async simulator for evaluating the maximum sustainable performances. Using 64 nodes of ivanoe platform, based on Intel Nehalem processors, the PaRSEC implementation with optimal partitioning, achieves 6.1 Tflop/s corresponding to 33.9% of the theoretical peak of the considered nodes.

Having assessed the performance of the sweep operation on modern massively parallel architectures, we can therefore build a massively parallel neutron transport solver on these architectures. This solver presented in Chapter 5.

87

In this thesis, we consider the discrete ordinates method for solving the neutron transport equation as presented in Chapter 1. We have seen that this method involves iterating over the scattering source for solving the monokinetic neutron transport equations. Each scattering iteration (also referred to as source iteration) involves a sweep over the spatial domain for a set of angular directions. As mentioned, this sweep operation gathers the vast majority of computations within the discrete ordinates method. In Chapter 4, we presented a massively parallel implementation of this sweep operation on distributed multicore-based supercomputers, using the PaRSEC task-based runtime system. We have shown that the efficiency of this implementation compares well with the predicted performances from theoretical performance models. In this chapter, we consider the whole S N algorithm as implemented in the Domino solver which integrates our task-based implementation of the sweep.

In section 5.1, we first recall the source iteration (SI) scheme and its convergence acceleration using the classical DSA in order to deal with strongly diffusive problems. Then, after recalling some limitations of the DSA in a parallel context, we present an implementation of a new acceleration scheme: PDSA, which extends the DSA scheme. In section 5.2, we present a validation study to assess the accuracy of the S N method and its performance on a class of benchmarks including 3D PWR full-core models.

Acceleration of source iterations (SI)

As mentioned in section 1.3.1, the convergence of the scattering iterations (Algorithm 3) is very slow in highly diffusive media (Σ s ≈ Σ t ), and an acceleration scheme must be used to remedy this issue. One of the widely used acceleration scheme in this case, is the Diffusion Synthetic Acceleration (DSA) [START_REF] Edward | Advances in Discrete-Ordinates Methodology[END_REF].

Diffusion Synthetic Acceleration (DSA) method

General presentation of the DSA

Here we just recall the basics of this method, and the reader can refer to the paper [START_REF] Edward | Advances in Discrete-Ordinates Methodology[END_REF] for more details regarding its effectiveness and the Fourier analysis characterizing its convergence properties. Let us define ϵ k+ 1 2 = ψψ k+ 1 2 as the error on the solution obtained after the k + 1 2 th iteration of the source iterations (SI) scheme, relative to the exact solution ψ, as defined by equation (1.15). By subtracting the equation on Line 3 of the Algorithm 3 in equation (1.15), the error ϵ on the angular flux satisfies the following transport equation:

Lϵ k+ 1 2 ( #" r , #" Ω) = ∫ S 2 d #" Ω ′ Σ s ( #" r , #" Ω ′ • #" Ω)ϵ k+ 1 2 ( #" r , #" Ω ′ ) + Σ s0 ( ϕ k+ 1 2 -ϕ k ) , ( 5.1)
which is as difficult to solve as the original fixed-source transport problem (1.15). However, if an approximate solution εk+ 1 2 of this equation was available, the scalar flux could be updated to:

ϕ k+1 = ϕ k+ 1 2 + εk+ 1 2 .
The idea of the DSA method is to use a diffusion approximation instead of solving transport equation (5.1). In Domino, the diffusion approximation relies on the simplified P N (SP N ) method. This method was firstly introduced in [START_REF] Gelbard | Simplified Spherical Harmonics Equations and Their Use in Shielding Problems[END_REF], but in our case we focus on the mixed-dual formulation as presented in [START_REF] Jj Lautard | Mixed dual methods for neutronic reactor core calculations in the CRONOS system[END_REF]. The choice of the diffusion approximation is firstly motivated by the simplicity of the diffusion operator, hence allowing to solve more easily the transport equation on the error (5.1). The second advantage of using this approximation is related to the error attenuation between successive iterations of the SI scheme. Indeed, the Fourier analysis performed in [START_REF] Edward | Advances in Discrete-Ordinates Methodology[END_REF] shows that DSA attenuates low-frequency errors left by the transport sweep operation. The latter operation attenuates almost only high-frequency errors in highly diffusive problems.

Solution of the diffusion problem

Equation (5.1) is approximated by a monogroup SP1 problem with an isotropic cross-section (Σ s0 ), as presented by the following problem in mixed dual formulation:

Problem 4. Find (ϵ, #" J ) ∈ L 2 (D) × H(D, div) such that:          div #" J ( #" r ) + Σ a ϵ( #" r ) = S( #" r ) in D #" ∇ ϵ( #" r ) + 1 D #" J ( #" r ) = #" 0 in D ϵ = 0 on ∂D (5.2)
where:

S( #" r ) = Σ s0 ( ϕ k+ 1 2 ( #" r ) -ϕ k ( #" r ) )
is the source term; Σ a = Σ t -Σ s0 the absorption cross-section, and D = 1 3Σa the diffusion coefficient.

The complete description of the SP N method is beyond the scope of this dissertation, but rather we just give its main characteristics as implemented in our SP N solver Diabolo [START_REF] Plagne | Generic Programming for Deterministic Neutron Transport Codes[END_REF]. Let us consider the general case of Problem 4 where the flux at the boundary is defined as:

ϵ = ϵ b on ∂D.
Thus, multiplying the first line of equation (5.2) by v ∈ L 2 (D), the second by #"

w ∈ H(D, div), and finally applying the Green formula, we obtain the following mixed-dual variational problem [START_REF] Coulomb | Mixed and Mixed-Hybrid Elements for the Diffusion Equation[END_REF]:

Problem 5. Find (ϵ, #" J ) ∈ L 2 (D) × H(D, div) such that:        ∫ D div( #" J ( #" r )) v( #" r )d #" r + ∫ D Σ a ( #" r )ϵ( #" r )v( #" r )dD = ∫ D S( #" r )v( #" r )d #" r , ∀v ∈ L 2 (D) ∫ D 1 D( #" r ) #" J ( #" r ) • #" w( #" r )d #" r - ∫ D ϵ( #" r )div( #" w( #" r ))d #" r = - ∫ ∂D ϵ b ( #" r ) #" w( #" r ) • #" n dΓ, ∀ #" w ∈ H(D, div),
(5.3) where #" n is the unit normal vector to the boundary ∂D.

Then, the Problem 5 is discretized spatially using the RTk mixed-dual finite element described in [START_REF] Nédélec | A New Family of Mixed Finite Elements in R3[END_REF][START_REF] Raviart | A Mixed Finite Element Method for 2nd Order Elliptic Problems[END_REF]. To apply the DSA, computations are made with the RT0 element, which has 7 (resp. 4) degrees of freedom (DoFs) per cell in 3D (resp. 2D): 1 DoF for the scalar unknown and 6 (resp. 3) DoFs for the vector unknown (see Figure 5.1).

The DSA acceleration scheme [2] was proven effective, provided that the spatial discretization of the transport equation is consistent with the spatial discretization of the diffusion solver [START_REF] Larsen | Unconditionally Stable Diffusion Synthetic Acceleration Methods for the Slab Geometry Discrete Ordinates Equations[END_REF]. In [START_REF] Hébert | High order diamond differencing schemes[END_REF], the author proved that this consistency requirement is lifted when using a Diamond Differencing scheme of order k for the transport equation, and RTk finite elements as a discretization scheme for the diffusion solver. As our SP N solver uses also RTk finite elements, we are assured of the stability of the DSA scheme. However, when integrated into a parallelized transport solver, DSA may become a bottleneck for the scalability of the transport solver if, for instance, a serial implementation of the diffusion solver is used. On the other hand, if the diffusion solver is parallelized, as presented in [START_REF] Barrault | Efficient Parallel Resolution of the Simplified Transport Equations in Mixed-Dual Formulation[END_REF][START_REF] Jamelot | Fast non-overlapping Schwarz domain decomposition methods for solving the neutron diffusion equation[END_REF] using a domain decomposition method, the iteration count to the solution increases with the number of subdomains, and can lead to a poor global scalability [START_REF] Yavuz | Iterative Methods for Solving x-y Geometry SN Problems on Parallel Architecture Computers[END_REF]. To remedy this issue, a variant of the DSA has been recently proposed by F. Févotte in [START_REF] Févotte | PDSA: a Piecewise Diffusion Synthetic Acceleration Scheme[END_REF].

Piecewise DSA method (PDSA)

The general presentation and the convergence proof of the PDSA method are given in [START_REF] Févotte | PDSA: a Piecewise Diffusion Synthetic Acceleration Scheme[END_REF].

General algorithm

We assume that the spatial domain D is split, along the 3 dimensions of the space, into N = P × Q × R non-overlapping subdomains D I such that: D = ∪ I∈I D I , where

I = 1, P × 1, Q × 1, R .
We set: Γ IJ = ∂D I ∩ ∂D J the non-empty interfaces between subdomains of index I and J; Γ I = ∂D ∩ ∂D I ; and #" n I the unit normal vector to ∂D I .

The first step of the PDSA method is, as in the case of the DSA method, an S N transport sweep operation on the whole spatial domain. The second step, consisting of setting an approximation of the error on the scalar flux, is split in two sub-steps: we successively solve two diffusion problems, on each subdomain, respectively with homogeneous Neumann boundary conditions (equation (5.4)), and non-homogeneous Dirichlet boundary conditions (equation (5.5)).

                 div #" J I N ( #" r ) + Σ a ϵ I N ( #" r ) = S I ( #" r ) in D I #" ∇ ϵ I N ( #" r ) + 1 D #" J I N ( #" r ) = #" 0 in D I ϵ I N = 0 on ∂Γ I #" ∇ϵ I N • #" n = 0 on Γ IJ (5.4)                    div #" J I D ( #" r ) + Σ a ϵ I D ( #" r ) = S I ( #" r ) in D I #" ∇ ϵ I D ( #" r ) + 1 D #" J I D ( #" r ) = #" 0 in D I ϵ I D = 0 on Γ I ϵ I D = ϵ I N + ϵ J N 2 on Γ IJ (5.5)
This is a major shift from the classical DSA method, as we are no longer required to get the solution of the diffusion problem on the whole spatial domain. The first advantage of this method is that the explicit global synchronizations between the resolutions of the piecewise diffusion problems are largely reduced, hence allowing to fully parallelize the DSA method without efficiency loss. In addition, as we are going to see in section 5.2.3, the effectiveness of the PDSA method is comparable to that of the classical DSA method on a class of benchmarks.

The accelerated S N flux is finally given by:

ϕ accel = ϕ + ϵ D .

Treatment of boundary conditions

The Neumann boundary conditions in equation (5.4) correspond to a reflective boundary condition #" ∇ϵ • #" n = 0, and its implementation poses no difficulty. However, as noted in [START_REF] Babuška | The Treatment of Nonhomogeneous Dirichlet Boundary Conditions by the p-Version of the Finite Element Method[END_REF], imposing non-homogeneous Dirichlet boundary condition such as in equation (5.5), represents several numerical challenges. Meanwhile, thanks to the mixed-dual variational formulation (5.3) we are able to bypass this issue. Indeed, considering the particular case of the second line of this variational formulation with the Dirichlet boundary conditions of equation (5.5), we obtain:

∫ D I 1 D( #" r ) #" J ( #" r )• #" w( #" r )d #" r - ∫ D I ϵ D ( #" r )div( #" w( #" r ))d #" r = - ∫ Γ IJ ϵ b ( #" r ) #" w( #" r )• #" n I dΓ, ∀ #" w ∈ H(D, div), (5.6)
where

ϵ b ( #" r ) = ϵ I N ( #" r ) + ϵ J N ( #" r ) 2 .
The functions ϵ I N , which come from the first SP N resolution, are naturally expanded on the RT0 finite element basis of as follows:

ϵ I N = ∑ i E i I v i ,
where E i I is the unknown vector of the flux over the subdomain, and v i are the basis functions for the flux (see Figure 5.2). On the other hand, equation (5.6) is also discretized using RT0 finite elements. The test functions for which it is evaluated are thus basis functions for the current unknown, hereafter denoted by ( #" w j ) j . Without loss of generality, we suppose that the interface Γ IJ is vertical (as seen for example on Figure 5.2). Therefore, #" w j • #" n = w x j on Γ IJ , and the integral in the right-hand side of equation (5.6) becomes:

B j = ∫ Γ IJ ϵ b ( #" r )w x j ( #" r )dΓ
Due to the expression of RT0 basis functions, we can notice that B j = E i I , where indices i and j are related as shown on Figure 5.2: the flux DoF indexed by i is associated to the cell whose boundary supports the current DoF indexed by j.

An illustration of the processing of the boundary conditions in the case of two subdomains is on Figure 5.3.

Parallelization of the PDSA Method

Figure 5.4 illustrates a parallel implementation of the PDSA method in 2D, when the global domain is partitioned in two subdomains. The partitioning of the global domain uses the same block data distribution as for the sweep operation. As we mentioned previously, the diffusion problem on each subdomain is solved using our SP N solver Diabolo which is parallelized on shared memory system using Intel TBB framework.

Hence, by mapping each subdomain to a single process, the resolution of the diffusion problems on D 1 and D 2 , when applying the PDSA method, is naturally done in parallel. Moreover, for the first step, the use of Neumann boundary conditions requires no communications with the neighboring processes. However, in the second step, each process needs to have the average value of the scalar flux at the interfaces between its neighbors. Therefore, each process must perform send and receive operations to exchange data with its neighbors. These data exchanges are point-to-point communications as only two processes are involved for each data exchange. To achieve this, we allocate two extra buffers per process: the first one is dedicated to store the extracted scalar flux at the interface between the subdomains which is then communicated to the neighboring process; while the second one is used as a reception buffer. We use asynchronous MPI communication primitives to exchange the flux at the interfaces.

Validation and performances of Domino

We assessed both the accuracy of the discretizations used in Domino and the efficiency of the PDSA method. This section presents experiments carried out to this end. We first start by presenting the different benchmarks used to perform our experiments.

Benchmarks

We used three different benchmarks to assess the accuracy and parallel performance of the Domino solver.

Kobayashi benchmarks [67]

These one-group benchmarks are used to assess the accuracy of the flux distribution, on geometries having void regions, in a highly absorbing medium. A full description of these benchmarks is available in [67]. Among these benchmarks, we considered the Problem 1 as depicted on Figure 5.5, and the case ii characterized by 50% of scattering. Table 5.1 recalls the values of the cross-sections and neutron source used.

Takeda benchmarks [116]

These benchmarks, consisting of four core models, are widely used for checking the validity of 3D neutron transport solvers. We considered the Model 1, which corresponds to a small Light Water Reactor (LWR) core. This model can be used with control rods inserted or not. In this study, we consider the former case, as depicted on Figure 5.6 (case 2). The two-group cross sections used in this model are characterized by a high scattering ratio (c ≃ 0.98765). Therefore, this model is a good candidate to show the effectiveness of the scattering acceleration methods. The boundary conditions associated with kobayashi and takeda benchmarks are reflective. However, since we consider only vacuum boundary condition in this study, we have replicated the spatial domain for both benchmarks in order to use vacuum boundary conditions. PWR 900 MW benchmarks These benchmarks correspond to a PWR 900 MW core, and enable to perform 2, 8 and 26 energy groups calculations. A full description of these benchmarks is available in [START_REF] Courau | Specifications of a 3D PWR Core Benchmark for Neutron Transport[END_REF]. It corresponds to a simplified 3D PWR first core loaded with 3 different types of fuel assemblies characterized by a specific Uranium-235 enrichment (low, medium and highly enriched uranium). There are no control rods inserted in this core model. Along the z-axis, the 360 cm assembly is axially reflected with 30 cm of water which results in a total core height of 420 cm. The 3 types of fuel assemblies appear on Figure 5.7 where the central assembly corresponds to the lowest enrichment, while the last row of fuel assemblies has the highest enrichment to flatten the neutron flux. Each fuel assembly is a 17 × 17 array of fuel pins, with a lattice pitch of 1.26 cm that contains 264 fuel pins and 25 water holes. The boundary condition associated with this benchmark problem is a pure leakage without any incoming angular flux. The associated nuclear data, 2-group, 8-group and 26-group libraries, derive from a fuel assembly heterogeneous transport calculation performed with the cell code dragon [START_REF] Marleau | A User's Guide for DRAGON 3[END_REF]. Table 5.2 summarizes the discretization parameters for the considered benchmarks.

Region

S Σ t Σ s (n cm -3 s -1 ) (cm -1 ) (cm -1 ) Problem i Problem ii 1 1 0.1 0 0.05 2 0 10 -4 0 0.5 × 10 -4 3 0 0.1 0 0.05
• For the kobayashi benchmark, we considered a uniform 2-cm mesh step, leading to a spatial mesh of 100 × 100 × 100 cells. We set the MacroCell sizes to 20 × 20 × 20, and we considered a S 16 angular quadrature.

• The spatial mesh size for the takeda benchmark is the same as that of the small test case we used in Chapter 4, and we use the same MacroCell size. The angular quadrature used for the takeda benchmark is a S 8 (80 angular directions).

• The spatial mesh used for the PWR benchmarks is based on a pin-cell mesh in the xy plane. Each pin-cell is then subdivided into 70 (resp. case enables to study the strong scalability of our implementation at high core count. For all the three benchmarks, we use MacroCells of size 20 × 20 × 20.

The calculation of DoF numbers consider 3 DoFs per cell, per energy group and per angular direction. ϵ k eff and ϵ ψ define the thresholds used to check the stopping criteria at iteration n + 1 of the power algorithm, respectively on the eigenvalue and on the fission source as follows:

|k n+1 eff -k n eff | k n eff < ϵ k eff , ||Fψ n+1 -Fψ n || ||Fψ n || < ϵ ψ .
(5.7)

The following experiments were conducted by launching one MPI process per computing node and as many threads as available cores; keeping one core per node for the communication thread. All experiments were conducted in single precision. Computation times do not include setup (reading of cross-section files from the hard disk), but include all communications and stopping criterion checks. For all the experiments presented in the following sections, the setup time is less than a minute.

Validation and performances of the source iterations scheme

In this section, we present the performances of the source iterations scheme (S N -only), without using the acceleration.

Kobayashi benchmark

To check the accuracy of the source iterations as implemented in Domino, we first consider the Problem 1ii of the kobayashi benchmarks. The computation settings are defined in Table 5.2. We ran this benchmark on a single 24-cores computing node of the athos platform, and the convergence was reached in 20 iterations. Figure 5.8 presents solutions obtained from Domino and the extracted reference results from [67]. The reference fluxes are defined for three set of points: A (Figure 5.8a and Figure 5.8b), B (Figure 5.8c and Figure 5.8d) and C (Figure 5.8e and Figure 5.8f). We observe that the flux obtained from Domino follow the same trend as reference values. However, the relative errors on the flux are higher: 0.9 along the line x = y = z (Figure 5.8d), at the mesh point defined by x = 95 cm. This observation is explained by the ray-effects as shown in [67]. It is possible to mitigate the ray effects using for example firstcollision source approximation methods or the Gauss-Legendre (GL) quadrature formula with large number of directions. Such a strategy is implemented in the radiation transport code Denovo [START_REF] Thomas M Evans | Denovo: A New Three-Dimensional Parallel Discrete Ordinates Code in SCALE[END_REF], enabling to get accurate solutions. Meanwhile, for our target applications (PWR core simulations), the ray-effects are negligible. Indeed, in a reactor core the neutron sources are uniformly distributed inside the reactor core, enabling to dramatically mitigate the ray effects.

Table 5.3 presents the computation time of the kobayashi problems. The computation time required to solve problem 1ii with Domino is 0.67 s, using a single 24-cores computing node of the athos supercomputer. For the same problem, Denovo code requires 3.1 s to get the solution using 16 processors [START_REF] Thomas M Evans | Denovo: A New Three-Dimensional Parallel Discrete Ordinates Code in SCALE[END_REF] of the Jaguar XT5 supercomputer. Meanwhile, one should keep in mind that because we do not take into account symmetry boundary conditions, the spatial mesh that we used is eight times larger. In addition, the processor architectures used in both cases are different. Therefore, the aim of this comparison is to give an order of magnitude of how Domino compares to another existing neutron transport code. 5.3: S N -only solution times for the Kobayashi problem 1ii. Domino was run on a single 24cores computing node (dual Intel Xeon E5-2697v2 processors) of the athos platform. Denovo was run on 16 processors of the Jaguar XT5 supercomputer, and correspond to the case where weighted diamond difference (WDD) spatial discretization scheme is used [START_REF] Thomas M Evans | Denovo: A New Three-Dimensional Parallel Discrete Ordinates Code in SCALE[END_REF].

Takeda benchmark

The second problem we considered for studying the performance of the source iterations is the Model 1 of the takeda benchmark. We used a spatial mesh resolution of 0.416 cm, and a S 8 Level Symmetric quadrature. The resulting spatial mesh contains 120 × 120 × 120 cells. For the external iterations, the stopping criterion on the eigenvalue and on the fission source term are respectively set to: ϵ k eff = 10 -6 and ϵ ψ = 10 -5 . As there is no up-scattering in the takeda benchmark, we set the number of Gauss-Seidel iterations (I g ) to one. Also, the number of source iterations (N src ) is fixed at one. Thereby, the convergence of the source iterations is determined by the convergence of the power algorithm on external iterations. We carried out this experiment on a single computing node of the athos platform, and the result is presented in Figure 5.9. The convergence is reached in 282 iterations and the obtained eigenvalue is 0.962386. This value is well within the Monte Carlo reference error bar (0.9624 ± 0.0005, see [START_REF] Takeda | 3-D Neutron Transport Benchmarks[END_REF]), therefore justifying the convergence of the source iterations scheme. The drop point in the convergence curve of the k eff , as depicted in Figure 5.9a, highlights a non-monotonic convergence on the eigenvalue. Indeed, the relative error on the k eff decreases until the iteration 7 and then changes the sign. In Figure 5.9b, the drop points correspond to the iterations where the Chebyshev acceleration method is applied (more details on the implementation of this method is available in [START_REF] Ponçot | Note de principe du solveur SN DOMINO[END_REF]). We evaluated the region-averaged flux obtained from Domino on the Model 1 of the takeda benchmarks. The relative discrepancies between the Domino solution and the reference solution is presented in Table 5.4. The maximum value of the discrepancies between the Domino solution and the reference solution is 0.44% (resp. 0.16%) on the fast (resp. thermal) group (see Table 5.4). These lower discrepancies illustrate that the convergence on the scalar flux is reached. [START_REF] Takeda | 3-D Neutron Transport Benchmarks[END_REF], for the Model 1 of the takeda benchmark.

|δϕ 1 | (%)
The study we presented in this section has validated the accuracy of the discretization schemes and the convergence of the source iterations scheme. In the next section, we are going to study the capabilities of the PDSA method for the reduction of the number of source iterations number.

Efficiency of the PDSA scheme

As mentioned in section 5.1.2, the application of the PDSA method after each source iteration, requires solving two successive SP N problems. The solutions to the SP N problems need to be converged to ensure the convergence of the PDSA scheme. As explained in [START_REF] Févotte | PDSA: a Piecewise Diffusion Synthetic Acceleration Scheme[END_REF], the convergence of the PDSA scheme can be determined using a parameter ρ max PDSA which depends on the scattering ratio c of the problem, and on the domain partitioning. This parameter, which represents the amplification factor of the PDSA scheme, is defined by equation (5.8):

ρ max PDSA = ρ max DSA + ρmax d R PDSA (θ), (5.8)
where: ρ max DSA is the spectral radius of the DSA scheme; ρmax d is an expression depending on the scattering ratio of the problem; and R PDSA (θ) is a closed-form formula which depends on a parameter θ, characterizing the optical thickness of subdomains and defined by:

θ = √ 3(1 -c)τ, ( 5.9)
where τ is the optical thickness of a subdomain1 . According to the study performed in [START_REF] Févotte | PDSA: a Piecewise Diffusion Synthetic Acceleration Scheme[END_REF], R PDSA converges towards 0 for optically thick subdomains. In this case, ρ max PDSA ≈ ρ max DSA , so that the efficiency of the PDSA scheme is comparable to that of the classical DSA scheme. It should be noted that the indicator in equation (5.8) is defined for 1D homogeneous cases. For a 3D, heterogeneous problem, R PDSA is evaluated by considering the lowest optical thickness over all subdomains, and over the three dimensions. The uniform optical thickness obtained in this way is somewhat penalizing; a more precise approach would consist in evaluating the optical thickness using a homogenized set of cross-sections per subdomain. In the following, we will use the ρ max PDSA indicator to highlight the efficiency of the PDSA scheme for multidimensional partitionings. Let us first study the case of a single subdomain.

Case of a single subdomain

When there is only a single subdomain, then the PDSA scheme is equivalent to the classical DSA scheme, except that the former requires one extra SP N resolution. Figure 5.10 shows the variation of outer iterations (N outer ) as a function of the number of SP N iterations (N SPN ), using the Model 1 of the takeda benchmarks (see section 5.2.1). We observe that by imposing one SP N iteration, the convergence is reached in 38 outer iterations. Furthermore, N outer decreases when N SPN increases, and reaches a minimal value of N outer = 14 from N SPN = 4, suggesting that the convergence of the SP N solver is reached with this value of N SPN . The corresponding eigenvalue is 0.962437 (see Table 5.5), which is slightly different (0.962386) than that was obtained without using the acceleration method.

Figure 5.10 also shows the total computation time as a function of N SPN . As expected, the computation time follows the same trend as N outer for the first four iterations, and reaches a minimal value of 2.43 s. However, after N SPN = 4 the total computation time increases because the number of outer iterations is then constant. This study illustrates that for a given problem, there exists an optimal value of N SPN which minimizes the total computation time.

k

Multidimensional case

According to a preliminary experimental study, we found that in a multidimensional case, the optimal value of N SPN required to minimize N outer may be greater than that obtained in the case of a single domain. Therefore, for the following study, we used N SPN = 15 to ensure the convergence of the SP N solver. The goal of this study is to assess the convergence of the PDSA scheme for multidimensional partitionings. We considered the takeda benchmark presented in section 5.2.1 and we evaluated the convergence of the solver for several partitionings (between (1, 1, 1) to (10, 10, 10)) of the global domain. For each partitioning, we evaluated ρ max PDSA , and N outer . The result is depicted in Figure 5.11. Each point corresponds to a given partitioning and thus is associated with a specific value of ρ max PDSA . The gradient colors are used to indicate the value of P + Q + R, for a partitioning (P, Q, R). For all the considered partitionings, the PDSA scheme converges and enables to reduce N outer as compared to the source iteration (SI). Furthermore, each value of ρ max PDSA is associated with several partitionings. Partitionings with larger values of P + Q + R are those giving larger number of iterations. This latter observation suggests that the ρ max PDSA indicator does not capture 3D effects. However, it seems to conservatively indicate partitionings where PDSA will converge.

Important remark. For a fixed number of cores, the partitioning of a spatial domain according to a Flat implementation will feature a larger number of subdomains as compared to that of a Hybrid approach. Therefore, the sum P + Q + R will be larger, and consequently the convergence of a Flat implementation will be slower (Figure 5.11) than that of a Hybrid implementation.

Figure 5.12 presents the convergence of the PDSA scheme for a multidimensional partitioning of (4,4,2). For this partitioning, we found that N SPN = 3 iterations for the SP N solver minimizes the total computation time which is 2.38 s, and the corresponding number of outer iterations is N outer = 62. According to this observation, the following results are obtained by selecting optimal values of N SPN .

Full-core 3D PWR calculations

In this section, we present full-core k eff computations using the 3D PWR core models described in section 5.2.1. The calculation parameters are defined in Table 5.2. From a preliminary study with a single subdomain, we found that the optimal number of SP N iterations is one, for each of the three benchmarks. Therefore, all the following results are obtained using this value.

2-group PWR core model

Table 5.6 presents a strong scalability study on the athos platform, using the 2-group PWR core model. As expected, for this benchmark characterized by a high scattering ratio (c ≈ 0.97), the convergence of the source iterations is very slow. For instance, without using the PDSA scheme, 2116 external iterations are required to reach the convergence with a partitioning of (2, 2, 1). The total computation time associated with this partitioning is 7724.4 s, of which 81.1% represents the time spent in the sweep operation, highlighting that the sweep operation is still dominant for the considered benchmark. Using the PDSA scheme, the number of external iterations is largely reduced to 92 corresponding to an acceleration of 23. Indeed, for this benchmark ρ max DSA ≈ 0.21 which is similar to the value of ρ max PDSA associated with the partitioning (2, 2, 1). Thereby, the convergence of the PDSA scheme is optimal. Furthermore, the total computation time drops to 451.6 s representing a speed-up (relative to the case without PDSA) of 17.1. It should be noted that this speed-up is less than the optimal value of 23, even if this is the case for the sweep operation (improvement by a factor of 22.4). This is justified by the additional computational costs associated with the PDSA scheme: the two SP N resolutions (T spn = 47.8 s) and the communications required to exchange data between subdomains (T comm PDSA = 35.7 s). Considering the run times when the PDSA scheme is used, we obtain that the performance of the sweep operation reaches a 2.7 Tflop/s on four computing nodes of athos, corresponding to 65.2% of the theoretical peak performance of the four nodes. Experiments with other partitionings yield similar trends. In particular, using 768 cores, the performance of the sweep operation is 8.3 Tflop/s, corresponding to 24.72% of the peak of the corresponding 32 nodes.

In Table 5.7, we report a performance comparison between Domino and the radiation transport code Denovo on the 2-group problem. We used 64 computing nodes of the athos supercomputer, distributed into a grid of 4 × 4 × 4 processes. The corresponding global computation time for the whole solver is 0.8 min. As a comparison, the Denovo code solves a slightly different version of this 2-group benchmark in 2.05 min on the Jaguar XT5 supercomputer [START_REF] Gregory G Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF] (18688 compute nodes, each with dual 2.6 GHz AMD 6-core Istanbul processor). One should note that this machine differs from athos and as a consequence the performance comparison is not very precise. In addition, the axial meshes used in Domino and Denovo are slightly different (756 vs 700). Nonetheless, we hope that this comparison provides a correct trend on how Domino compares to Denovo.

8-group PWR core model

Table 5.8 presents performance results of a S 12 8-group 3D PWR k eff computations using 64 computing nodes of the athos cluster partitioned into (4,4,4). As in the 2-group benchmark presented previously, ρ max PDSA with this partitioning is similar to ρ max DSA . The convergence on this benchmark is reached in 65 external iterations, and the obtained eigenvalue is k eff = 1.009408. This number of external iterations is similar to that was obtained for a run with a single subdomain. The total computation time is 128. [START_REF] Mckee | Reflections on the Memory Wall[END_REF] (4, 4, 4). As in Figure 5.9b, the drop points in Figure 5.13b correspond to the iterations where the Chebyshev acceleration method is applied. The detailed computing times for this benchmark are reported in ). The obtained eigenvalue is k eff = 1.008358. As in the case of 8-group, we did not observe any increase on the number of external iterations as compared to a run with a single domain. This is a remarkable result, highlighting the perfect efficiency of the PDSA method on representative benchmarks of our target applications.

▶ In this chapter, we studied the accuracy and the performances of our massively parallel approach for solving the neutron transport equation according to the discrete ordinates method. We first integrated our task-based implementation of the sweep with PaRSEC in the Domino solver, then we integrated a new acceleration method to speed-up the convergence of the scattering iterations in strongly diffusive media as PWR cores. All the experiments were carried-out using Domino solver, and we used three classes of benchmarks.

We then considered one of the kobayashi benchmarks series in order to assess the accuracy of the discretization schemes and to study the performances of the source iterations. We considered the case with 50% of scattering. We found that for some regions, the discrepancies between the flux obtained from Domino and the reference solution are particularly large (rel. error of 2.4) due to ray-effects. However, one should note that for our target applications, which are PWR core simulations, the ray-effects are not present. Indeed, the PWR core corresponds to a strongly diffusive medium where the neutron sources are uniformly distributed, enabling to dramatically mitigate the ray effects. The runtime performances of the Domino solver compares favorably with the Denovo on the considered kobayashi problem.

The second benchmark that we studied is the takeda benchmark where the control rods are inserted. On this benchmark with a high scattering ratio, we assessed both the accuracy of the solution and the performances of the PDSA scheme. We highlighted the correlation between the convergence of PDSA scheme with an indicator, ρ max PDSA , depending on the partitioning of the domain. The lower bound of this indicator is the spectral radius of the classical DSA method. Therefore, we have found that the convergence of the PDSA scheme improves for partitionings with lower values of this indicator. According to our study, we found that even though ρ max PDSA is defined for 1D cases, it enables us to conservatively predict the convergence of 3D heterogeneous problems.

Finally, we have studied the performances of Domino on 3D PWR core models. For these benchmarks, the convergence of PDSA scheme is optimal and requires roughly the same number of iterations as a global DSA scheme to reach the convergence. Thereby, the performance of the Domino solver is very satisfactory. For instance, on a 2-group k eff computation, the performance of the sweep operation reaches 8.3 Tflop/s using 768 cores of the athos platform. This performance corresponds to 24.72% of the the theoretical peak of the considered nodes.

Conclusion and Future Work Conclusion

The goal of the research presented in this thesis was to study, and to propose a suitable solution, to the challenges posed by the use of hierarchical massively parallel computers, for solving the neutron transport equation according to the discrete ordinates method (S N ). The addressed challenges were manifold: design of efficient algorithms capable to handle the SIMD paradigm on modern CPUs; efficient utilization of multicore-based clusters by means of emerging taskbased models on top of generic runtime systems; and the use of efficient numerical methods.

From this perspective, we have designed theoretical performance models of the sweep operation intended to bound its performances both on a single CPU and on hierarchical multicorebased architectures. These performance models have been used to establish a strategy that we have adopted to meet the above mentioned challenges.

At the processor core level, we have shown that the SIMD capabilities of modern processors impose some constraints on the design of efficient computational kernels. In the case of the sweep operation, we have studied different strategies for the vectorization of the computations. We have theoretically proved that the efficiency of the vectorization over the spatial variable is limited by the padding that must be used to ensure data alignment. However, the vectorization over the angular variable allows to reduce the overhead due to the paddings. In order to explicitly take advantage of the vector units, we have studied and validated the use of generic programming models in C++ for exploiting directly the SIMD units corresponding to the target architecture. We have shown that this strategy allows to preserve the modularity of the code and to enhance its performance portability. On a Westmere processor, the performance of the angular vectorization reaches a maximum of 13.6 Gflop/s on a single CPU, corresponding to 62.9% of the theoretical peak performance of that CPU.

From this efficient single threaded implementation of the sweep kernel, we addressed the challenges posed by the parallel execution of the sweep on multi-processor architectures. We first built an accurate parallel sweep simulator in order to explore the different parallel spatial decomposition of the transport sweep. Furthermore, we used our sweep simulator to justify the need for a task-based implementation of the sweep operation in order to maximize its performances on multicore-based architectures. Then, we compared different emerging task-based models on top of generic runtime systems (Intel TBB, StarPU, PaRSEC). As a result, the PaRSEC framework based on parametrized DAG model, produced the most efficient implementation for the Cartesian transport sweep. This PaRSEC based implementation helped us to show that our performance model accurately predicts optimal partitionings. Using optimal partitioning, the performance of the sweep operation reaches 6.1 Tflop/s of 768 cores of the ivanoe supercomputer, which corresponds to 33.9% of the theoretical peak performance of this set of computational resources. 107 Finally, we addressed the challenge of converging the scatter iterations in highly diffusive media such as the PWR cores. We have implemented and studied the convergence of a new acceleration scheme (PDSA) that naturally suits our Hybrid parallel implementation. The efficiency of the PDSA scheme have been investigated on the reference takeda benchmark, according to a convergence indicator defined by the PDSA scheme. This indicator depends on the optical thickness of the subdomains and has a lower bound equal to the spectral radius of the classical DSA scheme. We have shown that when this indicator is similar to the spectral radius of the classic DSA scheme, as in the case of PWR cores, the number of external iterations required for convergence is the same as that of DSA.

The combination of all these techniques have enabled us to develop a massively parallel version of the Domino solver. It is capable of tackling the challenges posed by the neutron transport simulations and compares favorably with state-of-the-art solvers such as Denovo. For a typical 26-group PWR calculations involving 1.02 × 10 12 DoFs, the time to solution required by the Domino solver is 46 min using 1536 cores. Consequently, this Domino solver can be used by nuclear power plant operators such as EDF for improving the efficiency and safety of nuclear power plants.

Future work

The work presented in this thesis has raised several open questions for future research.

One such question is related to the improvement of the accuracy of the spatial discretization scheme by extending the DD0 scheme to high-order diamond differencing schemes (DD1 and DD2), or using a Discontinuous Galerkin (DG) scheme [START_REF] Cockburn | The Development of Discontinuous Galerkin Methods[END_REF][START_REF] Reed | Triangular Mesh Methods for the Neutron Transport Equation[END_REF][START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation[END_REF]. The main advantage of using these high-order discretization schemes is the possibility to obtain a comparable, or even better, level of accuracy with a coarser spatial mesh. However, for ensuring the stability of the PDSA scheme, the spatial discretization scheme must be consistent with the RTk finite element used for solving SP N equations. This scheme consistency requirement is ensured for diamond differencing scheme as shown in [START_REF] Hébert | High order diamond differencing schemes[END_REF], but there exists no equivalent result for the DG scheme.

Another question is oriented towards the design of efficient kernels for accelerator-based architectures. Early works on this topic such as in [START_REF] Gong | GPU Accelerated Simulations of 3D Deterministic Particle Transport Using Discrete Ordinates Method[END_REF][START_REF] Jamelot | High performance 3D neutron transport on petascale and hybrid architectures within APOLLO3 code[END_REF] require the use of quadrature formula featuring larger numbers of angular directions in order to achieve acceptable performance on GPU-based clusters. Recently, a new multilevel decomposition in energy has been introduced in [START_REF] Thomas M Evans | Three-dimensional discrete ordinates reactor assembly calculations on GPUs[END_REF], which is based on the use of Krylov methods for solving multigroup equations. This decomposition enables to extract a larger number of parallel tasks for problems having large number of energy groups, which is a requirement for maximizing the occupation of GPUs. However, for problems with few angles and energy groups, a solution for maximizing the efficiency would consists of merging the angular and spatial variables.

In Chapter 3, we have presented our Hybrid model and the Hybrid-Async simulator that can be used to predict the computation time of the sweep operation. One of the parameters used in these models is T grind : the computation time per cell, per group, and per angular direction. This parameter was obtained through experimental measurements. A challenge raised by this process is the sensitivity of the model predictions to T grind . Therefore, it is of great interest to accurately estimate this parameter according to machine-dependent constants such as sizes and hierarchy of caches, clock rate, instruction and the microarchitecture of the target CPU.

Furthermore, our performance models can be extended to describe the behavior of the whole S N algorithm including the computation time required by the SP N resolutions used by the PDSA scheme. To achieve this, it is required to develop a performance model of the SP N method as in [START_REF] Pinchedez | Calcul parallèle pour les équations de diffusion et de transport homogènes en neutronique[END_REF].

In a previous study [START_REF] Courau | DOMINO: A Fast 3D Cartesian Discrete Ordiantes Solver for Reference PWR Simulations and SPN Validations[END_REF], on shared memory systems, we found that initializing the external iterations of the S N algorithm by a flux obtained from a SP N calculation reduces the number of external iterations required by the transport solver by a factor up to 2.5. Such a strategy can therefore significantly reduce the global computation time. To efficiently take advantage of this result, it is required to have a distributed implementation of the SP N solver as in [START_REF] Barrault | Efficient Parallel Resolution of the Simplified Transport Equations in Mixed-Dual Formulation[END_REF][START_REF] Jamelot | Fast non-overlapping Schwarz domain decomposition methods for solving the neutron diffusion equation[END_REF][START_REF] Lathuilière | Méthode de décomposition de domaine pour les équations du transport simplifié en neutronique[END_REF].

Finally, in this thesis we have studied the stationary form of the neutron transport equation. Hence, this work provides an entry point towards time-dependent neutron transport simulation as 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]. However, time integration requires the storage of the angular flux. As we have shown in Chapter 2, the storage of the currents decreases the arithmetic intensity of the sweep kernel. Therefore, efficient time-dependent neutron transport simulation represents a big challenge for future research.

1 . 1

 11 Analysis of nuclear reactor cores . . . . . . . . . . . . . . . . . . . . . 1.1.1 Basic concepts of nuclear reactor physics . . . . . . . . . . . . . . . . . . 1.1.2 Boltzmann transport equation . . . . . . . . . . . . . . . . . . . . . . . 1.2 Multigroup formulation of the transport equation . . . . . . . . . . . 1.3 Discrete ordinates method on Cartesian meshes . . . . . . . . . . . . 1.3.1 Angular discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Spatial discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Cartesian transport sweep operation . . . . . . . . . . . . . . . . . . . . 1.4 Modern parallel computers and performance evaluation . . . . . . . 1.4.1 Architecture and design of parallel computers . . . . . . . . . . . . . . . 1.4.2 Metrics for performance evaluation . . . . . . . . . . . . . . . . . . . . .

  quadrature in 2D. Directions are grouped in quadrants. (b) In each direction, cells have two incoming components of the flux (Here from the left and bottom faces: ψL and ψB), and generates two outgoing components of the flux (Here on the right and top faces: ψR and ψT ).

Figure 1 . 2 :

 12 Figure 1.2: Illustration of the sweep operation over a 6 × 6 2D spatial grid for a single direction.

Figure 2 . 1 :

 21 Figure 2.1: Sweep of a single spatial cell (in DD0) in 2D, with a single angular direction.

Figure 2 . 2 :

 22 Figure 2.2: Sweep over a 5 × 4 2D spatial domain. It requires 20 scalar steps to go through the whole domain in scalar mode.

Figure 2 .

 2 Figure2.2 shows a sweep moving from the bottom-left corner of the 2D grid to the upperright, where all the cells on the same front can be processed at once. Our aim here is to use vector instructions to parallelize the processing of all the cells belonging to the same front. However, the number of cells on the fronts is not constant. Hence, we cannot directly use vector instructions as they operate on fixed-size registers. For this reason, we must add padding cells to the fronts 0, 1, 2 and 5, 6, 7 in order to have fixed-size vectors that fit the SIMD width of the target processor. This strategy is illustrated on Figure2.3. It shows that as the sweep moves

Figure 2 . 3 .Figure 2 . 3 :

 2323 Figure2.2 shows a sweep moving from the bottom-left corner of the 2D grid to the upperright, where all the cells on the same front can be processed at once. Our aim here is to use vector instructions to parallelize the processing of all the cells belonging to the same front. However, the number of cells on the fronts is not constant. Hence, we cannot directly use vector instructions as they operate on fixed-size registers. For this reason, we must add padding cells to the fronts 0, 1, 2 and 5, 6, 7 in order to have fixed-size vectors that fit the SIMD width of the target processor. This strategy is illustrated on Figure2.3. It shows that as the sweep moves

Listing 2 . 2 :

 22 The SIMD implementation of the sweep algorithm using Eigen. // BlockArray is a type representing an array of packed data typedef Eigen :: Array < RealType , packSize , 1 > BlockArray ; /* BlockArrayView maps a BlockArray expression to an existing array of data */ typedef Eigen :: Map < BlockArray , Eigen :: Aligned > BlockArrayView ; typedef Eigen :: Map < const BlockArray , Eigen :: Aligned > C on s t Bl o c kA r r ay V i ew ; const BlockArray sijk ( BlockArray :: Constant ( source [ k ][ j ][ i ])); const int nblocks = dir ec tionPerOct ant / packSize ; for ( int b =0; b < nblocks ; b ++){ const int dir = b * packSize ; BlockArray denom ( sigmaIJK ); BlockArray epsX = C on s tBlockArrayVi e w (& omegaX [ dir ])* twoInvStepX ; BlockArray epsY = C on s tBlockArrayVi e w (& omegaY [ dir ])* twoInvStepY ; BlockArray epsZ = C on s tBlockArrayVi e w (& omegaZ [ dir ])* twoInvStepZ ; BlockArrayView psiX (& psiXd [ dir ]); BlockArrayView psiY (& psiYd [ dir ]); BlockArrayView psiZ (& psiZd [ dir ]); BlockArray psiOut00 ( sijk ); psiOut += epsX * psiX + epsY * psiY + epsZ * psiZ ; BlockArray denom ( sigmaIJK ); denom += epsX ; denom += epsY ; denom += epsZ ; psiOut /= denom ; // Update of the scalar flux phi += psiOut * Co n s tB lockArrayVi ew (& weight [ dir ]); } The neutron currents are not stored

Figure 2 . 4 :

 24 Figure 2.4: Single precision arithmetic intensity as a function of the MacroCell size (n x = n y = n z ) and angular quadrature.

Figure 2 . 5 :

 25 Figure 2.5: Roofline model for an Intel Xeon E78837 processor and measured performances of the sweep operation as a function of the arithmetic intensity I a . The variation of I a is obtained by varying size of MacroCells as in Figure 2.4b.

Figure 2 . 6 :

 26 Figure 2.6: Roofline model for an Intel Xeon E52697 V2 processor and measured performances of the sweep operation as a function of the arithmetic intensity I a . The variation of I a is obtained by varying size of MacroCells as in Figure 2.4b.

3 . 1

 31 Preliminary definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Literature review on performance models of the sweep . . . . . . . . 3.2.1 Flat models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Hybrid models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 New performance model of a parallel sweep . . . . . . . . . . . . . . 3.3.1 Computation steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Communication steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Asynchronous simulator of the sweep . . . . . . . . . . . . . . . . . . 3.4.1 General presentation of the simulation algorithm . . . . . . . . . . . . . 3.4.2 Communication costs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Comparative study of the performance models . . . . . . . . . . . . . 3.5.1 Parameters of the models . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Evaluation of the performance models . . . . . . . . . . . . . . . . . . . 3.5.2.1 Optimizing process grid and task-granularity . . . . . . . . . . 3.5.2.2 Performance comparisons of the models . . . . . . . . . . . . . 3.5.2.3 Impact of the scheduling overhead in the Hybrid-Async model 39

  4): 1. a new performance model of the sweep, targeting distributed multicore systems; 2. a simulator of an asynchronous implementation of the sweep.

  a) Illustration of the blocked data distribution. MacroCells of similar colors belong to the same process. (b) Communication pattern of the sweep for a 6 × 6 MacroCell grid, over a 3 × 3 process grid (Ax = Ay = 7; Sx = Sy = 6; Lx = Ly = 2; F = 11).

Figure 3 . 1 :

 31 Figure 3.1: 2D blocked data distribution and communication pattern of the simulateneous sweep of all quadrants.

  where (a o , b o , c o ) represents the coordinates of the MacroCell (a, b, c) relative to the octant o

Figure 3 . 2 :

 32 Figure 3.2: Illustration of the formula giving the number of computation steps, for a single quadrant sweep, over a 2D spatial grid of 6 × 6 MacroCells. This sweep is performed on a dual-core computing node.

Figure 3 . 3 :

 33 Figure 3.3: Verification of the computation steps formula for a concurrent sweep of all quadrants over a 6 × 6 MacroCell grid, using nine dual-core computing nodes distributed into a 3 × 3 process grid. Tasks of similar colors are on the same node and follow distribution on Figure 3.1.

7 ▷ 8 ▷

 78 Update the clock of the communication thread t (p,q) c Perform communications and schedule ready tasks 9 clock+ = T task ; 10 return clock;

Figure 3 . 4 :

 34 Figure 3.4: Dependencies for the execution of a task in 2D.

1 ,▷

 1 Communication overhead (o c ) Out: Sweep execution time (clock) clock = 0.0; ▷ Global clock RemainingTasks = 4 × S x × S y ; while RemainingTasks > 0 do for w = N (p, q) ∈ N do ▷ Execute n ready tasks on w: T q,a,b • • • , T q,a,b n (n ≤ N ) RemainingTasks -= n; Update the clock of the communication thread T c t

Figure 3 . 5 :

 35 Figure 3.5: Performance of the network interconnect of the ivanoe platform. Measurements are obtained using NetPIPE utility, compiled with OpenMPI 1.6.5. The two vertical lines define the range to which belong messages sent through the network for the benchmarks considered in this study.

Figure 3 . 6 :

 36 Figure 3.6: Average performance of the sweep operation per MacroCell, as a function of the MacroCell size. This experiment was carried out using a sequential implementation of the sweep operation on the bigmem computing node.

Figure 3 . 7 :

 37 Figure 3.7: Performance of the Hybrid model using the small test case on ivanoe. This figure illustrates how the Hybrid model integrates the overlap of communications by computations.

Figure 3 . 8 :

 38 Figure 3.8: Comparison of Adams et al., Hybrid and Hybrid-Async models using the small test case using parameters of the ivanoe platform.

Figure 3 .

 3 Figure 3.11 presents a sensitivity study, of the Hybrid-Async model, to the scheduling overhead.This study is carried out with the small test case, using 64 nodes of the ivanoe platform. We fixed the task granularity and the process grid partitioning to the optimal values obtained from a run without any scheduling overhead (see Table3.4). We made the following observations: the performance of the sweep operation decreases as the scheduling overhead increases, with a large drop of performance (from 2.1 Tflop/s to 1.6 Tflop/s) when using δ = 57.43 µs. This value coincides with the computation time of a single task for the considered task granularity.

Figure 3 . 9 :

 39 Figure 3.9: Performance of the Hybrid model using the big test case on ivanoe.

Figure 3 . 10 :

 310 Figure 3.10: Comparison of Adams et al., Hybrid and Hybrid-Async using the big test case on ivanoe.

Figure 3 . 11 :

 311 Figure 3.11: Impact of the scheduling overhead on the small test case, using 64 computing nodes. Process grid partitioning and MacroCell sizes are fixed to the optimal values obtained from the run without any overhead (see Table3.4). The vertical line corresponds to the computation time of a single task.

Listing 4 . 1 :

 41 JDF file of the 2D sweep for one quadrant ComputePhi (a , b ) /* Execution space */ a = 1 .. ncx b = 1 .. ncy /* Parallel partitioning */ : mcg (a , b ) /* Parameters */ RW PSIX <-( a != 1) ? PSIX ComputePhi (a -1 , b ) : psi_x ( b ) -> ( a != ncx ) ? PSIX ComputePhi ( a +1 , b ) : psi_x ( b ) RW PSIY <-( b != 1) ? PSIY ComputePhi (a , b -1) : psi_y ( a ) -> ( b != ncy ) ? PSIY ComputePhi (a , b +1) : psi_y ( a ) RW MCG <-mcg (a , b ) -> mcg (a , b ) /* Priority of this task */ ; priority (a , b ) BODY { computePhi_CPU ( MCG , PSIX , PSIY ); } END

Listing 4 . 2 :

 42 Blocked data distribution of the MacroCell grid over the P × Q grid of processes. // Rank of the process owning MacroCell (a , b ) int rank_of ( Parameters & param , int a , int b ) { int lp = a / ( param . ncx / param . P ); int lq = b / ( param . ncy / param . Q ); return lq * param . P + lp ; } // Address of the object MacroCell (a , b ) void * data_of ( Parameters & param , int a , int b , DataType & mcgData ) { int aa = a % ( param . ncx / param . P ); int bb = b % ( param . ncy / param . Q ); return mcgData [ bb ][ aa ]; }

  (a) LFQ (b) PBQ -PlaneZ (c) PBQ -Front

Figure 4 . 1 :

 41 Figure 4.1: Animation snapshots showing the behavior of various scheduling strategies on the sweep progress throughout the spatial domain for one single direction. Data is distributed over a 2 × 2 × 2 process grid. Threads of similar colors belong to the same node. Brightest colors are used to show what MacroCells are processed during the snapshots.

Listing 4 . 3 :

 43 Basic submission of the sweep task-graph with StarPU. // ncx , ncy and ncz : number of MacroCells along the x , y and z axes const int nfront = ncx + ncy + ncz -2; // An octant o is defined by ( fx , fy , fz ) such that o = fx +2* fy +4* fz for ( fx =0; fx <2; fx ++) { int xinc = 1 -2* fx ; int xbeg = fx * ( ncx -1); int xend = ((1+ fx )%2)* ncx + fx * xinc ; for ( fy =0; fy <2; fy ++) { int yinc = 1 -2* fy ; int ybeg = fy * ( ncy -1); int yend = ((1+ fy )%2)* ncy + fy * yinc ; for ( fz =0; fz <2; fz ++) { int zinc = 1 -2* fz ; int zbeg = fz * ( ncz -1); int zend = ((1+ fz )%2)* ncz + fz * zinc ; for ( x = xbeg ; x != xend ; x += xinc ) { for ( y = ybeg ; y != yend ; y += yinc ) { for ( z = zbeg ; z != zend ; z += zinc ) { void * data_mcg = NULL ; if ( rank == mcg -> rank_of ( mcg , x , y , z ) ){ data_mcg = mcg -> data_of ( mcg , x , y , z ); } /* Create and submit the task T o,x,y,z to StarPU */ st ar pu _t as k_insert ( cl_solve , STARPU_VALUE , & data_mcg , sizeof ( void *) , STARPU_VALUE , & CS , sizeof

Figure 4 . 2 :

 42 Figure 4.2: Single-core performances of the PaRSEC, StarPU and Intel TBB implementations of the sweep operation, averaged by the number of MacroCells: 25 (there are five MacroCells along each dimension). This experiment was conducted on a computing node of the athos platform (Intel Xeon E5-2600 V2 processors). The angular discretization is a S 16 Level Symmetric quadrature (288 directions).

Figure 4 . 3 :

 43 Figure 4.3: Performance of the network interconnect of the ivanoe and athos platforms. Measurements are obtained using NetPIPE utility, compiled with OpenMPI 1.6.5. The two vertical lines define the range to which messages sent through the network belong to, for small and big test cases, according to task granularity.

Figure 4 . 4 :

 44 Figure 4.4: Performances of the implementations of the sweep operation using big test case (see Table3.2), on a single 24-cores computing node of the athos supercomputer. All measurements were obtained using the same task-granularity of 20 × 20 × 20 cells per task.

  (o c =0) Hybrid-Async (o c =1.8e-10) PaRSEC Hybrid Model (o c =1.8e-10)

Figure 4 . 5 :

 45 Figure 4.5: Sensitivity of Hybrid and Hybrid-Async performance models to data distributions using 48 nodes (576 cores) of the ivanoe cluster, and the big test case. Predicted performances by the Hybrid model and Hybrid-Async simulator are obtained using a scheduling overhead value of δ = 0 s.

Figure 4 . 6 : 4 . 5 . 3 . 3 Figure 4 . 7

 46453347 Figure 4.6:Comparison of Hybrid and Flat approaches using small test case on the ivanoe cluster. We used 64 nodes and 6 computations cores per node. The performance measurements were obtained using optimal partitioning, according to predictions of Hybrid and Adams et al. models.

Figure 4 . 7 :

 47 Figure 4.7: Performance comparison of the PaRSEC implementation of the sweep against predicted performances by Hybrid model and Hybrid-Async simulator. This experiment was conducted on the ivanoe cluster, using the small test case. Using δ = 8.1 • 10 -6 s, o c = 4.2 • 10 -10 s/Byte, the prediction of the Hybrid-Async complies with actual measurements.

Figure 4 . 8 :

 48 Figure 4.8: Performance of the PaRSEC implementation of the sweep using big test case on the ivanoe cluster.

1

 1 Acceleration of source iterations (SI) . . . . . . . . . . . . . . . . . . 5.1.1 Diffusion Synthetic Acceleration (DSA) method . . . . . . . . . . . . . . 5.1.2 Piecewise DSA method (PDSA) . . . . . . . . . . . . . . . . . . . . . . 5.2 Validation and performances of Domino . . . . . . . . . . . . . . . . 5.2.1 Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Validation and performances of the source iterations scheme . . . . . . . 5.2.3 Efficiency of the PDSA scheme . . . . . . . . . . . . . . . . . . . . . . . 5.2.4 Full-core 3D PWR calculations . . . . . . . . . . . . . . . . . . . . . . .

Figure 5 . 1 :

 51 Figure 5.1: RT0 finite element in 2D: 5 DoFs (4 for the current and 1 for the scalar flux).

Figure 5 . 2 :

 52 Figure 5.2: Correspondance between flux and current DoFs.

Figure 5 . 3 :

 53 Figure 5.3: Illustration of the PDSA method on a domain split in two. The first step consists of solving two diffusion problems in parallel on D 1 and D 2 , with Neumann boundary conditions. The second step solves also two diffusion problems with non-homogeneous Dirichlet boundary conditions: null flux boundary conditions on the external boundary of the domain and an average value of the flux at the inner interface.

Figure 5 . 4 :

 54 Figure 5.4: Illustration of the communication pattern in PDSA method on a domain split in two. Two point-to-point communications are needed to exchange flux at the interface between the two subdomains.

Figure 5 . 5 :

 55 Figure 5.5: Configuration of the kobayashi benchmark [67] (Problem 1).

Figure 5 . 6 :

 56 Figure 5.6: Core configuration the takeda benchmarks [116].

  1Cii -Relative error on the flux

Figure 5 . 8 :

 58 Figure 5.8: Comparison of the neutron fluxes for the Kobayashi problem 1ii (50% of scattering), obtained from Domino, with reference values obtained from GMVP Monte-Carlo code [67].

  Convergence on the source term.

Figure 5 . 9 :

 59 Figure 5.9: Convergence of the inverse power algorithm, without the acceleration, using takeda benchmark.

Figure 5 . 10 :

 510 Figure 5.10: Convergence of the PDSA scheme using a partitioning of (1, 1, 1).

Figure 5 . 11 :Figure 5 . 12 :

 511512 Figure 5.11: Convergence of the PDSA scheme as a function of ρ max PDSA .

Figure 5 .Figure 5 . 13 :

 5513 Figure 5.13 presents the convergence of a k eff computation on the 26-group 3D PWR core model, using the PDSA scheme. We used 64 computing nodes of the athos cluster partitioned into

Fission source computation 4 ▷ External iterations: inverse power algorithm 5 while Non convergence do 6

  g .ϕ g ; ▷

Multigroup iterations: Gauss-Seidel 8 while Non convergence do 9 for g ∈ 1, G do 10

	11	▷ External sources Q ext = S g + ∑ Σ g ′ →g s	.ϕ g ′ ;
		g ′ ̸ =g

12 ▷ Scattering iterations 13 while Non convergence do 14 Q = Q ext + Σ g→g s ϕ g ; 15 15

2.1 Review of the SIMD paradigm . . . . . . . . . . . . . . . . . . . . . . 20 2

  .1.1 General presentation of the SIMD paradigm . . . . . . . . . . . . . . . . 20 2.1.2 On the way for vectorization . . . . . . . . . . . . . . . . . . . . . . . . 23 2.1.3 Generic programming and tools for the vectorization . . . . . . . . . . .

23 2.2 Arithmetic intensity of the sweep kernel . . . . . . . . . . . . . . . . 24 2

  .2.1 Memory traffic and flops per cell and per direction . . . . . . . . . . . . 24 2.2.2 General formula of the arithmetic intensity of the sweep . . . . . . . . .

26 2.3 Vectorization over spatial domain . . . . . . . . . . . . . . . . . . . . . 28

	2.3.1 The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
	2.3.2 Maximum theoretical speed-up . . . . . . . . . . . . . . . . . . . . . . . 30
	2

.4 Vectorization over the angular variable . . . . . . . . . . . . . . . . . 30 2.4.1 The algorithm

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.4.2 Study of the arithmetic intensity of the sweep . . . . . . . . . . . . . . . 31 2.4.3 Actual performances vs Roofline model . . . . . . . . . . . . . . . . . . 34 19

Table 2 .

 2 .1 shows the evolution of the SIMD width for Intel processors. The number of elements that fit into a SIMD register of a given architecture can also be under-1: SIMD width as a function of the processor architecture.

	Architecture Launch date SIMD width	Elements per register	ISA
			(bits)	single precision double precision
	Westmere	01/2010	128	4	2	sse
	Ivy Bridge	04/2012	256	8	4	avx
	Skylake	08/2015	512	16	8	avx-512

  .2 we give the values of I c for our test machines: the bigmem computing node (Intel Xeon E78837) and a node of the athos platform (Intel Xeon E52697 V2). We distinguish the cases of a single-core,

			bigmem			athos
		Stream	Theoretical		Stream	Theoretical
		Bandwidth	Peak	I c	Bandwidth	Peak	I c
		(GB/s)	(GFlop/s)	(F/B)	(GB/s)	(GFlop/s)	(F/B)
	single-core	5.8	21.2	3.6	12.6	43.2	3.4
	single-socket	20.3	170.2	8.3	34.4	518.4	15.0
	full node	76.3	680.9	8.9	67.5	1036.8	15.3
	Table 2.2: Critical arithmetic intensities of Intel Xeon E78837 and Intel Xeon E52697 V2 pro-
	cessors.

Table 2 .

 2 .4). Except in the case

		N dir	Directions per octant (M )	Sustainable speed-up
			w/o	w/
			padding	padding
				sse	avx	sse	avx
	S 2	8	1	4	8	1	1
	S 4	24	3	4	8	3	3
	S 8	80	10	12	16	3.3	5
	S 12	168	21	24	24	3.5	7
	S 16	288	36	36	40	4	7.2

4: Impact of the padding on the vectorization efficiency in single precision. The sustainable speed-up gives the maximum attainable speed-up taking into account the padding.

Table 2 .

 2

	)

5: Illustration of the speed-up brought by the angular vectorization of the sweep, as a function of the MacroCell size (n x = n y = n z ), and for different angular quadrature orders. The performance measurements were carried-out on an Intel Xeon E78837 processor.

Table 2 .

 2 4) because of the padding we used. This slow speed-up is justified by the fact that for a MacroCell of size 1 × 1 × 1, as mentioned in section 2.2.2, the sweep kernel is memory bound and thus the computational gain brought by the vectorization is overshadowed by a slower DRAM bandwidth. The same observation applies to the cases of a sweep on a single cell with S 8 , S 12 and S 16 quadratures. As expected, the speed-up brought by the angular vectorization increases with the MacroCell size, or equivalently with the arithmetic intensity, and reaches 2.4, 3.78, 4.23, 4.44 respectively for the four angular quadratures, when the MacroCell size is 64 × 64 × 64 cells. One should note that the speed-up of 3.78 in the case of S 8 with n x = 64 is larger than 3.3 which was expected (see

Table 2

 2 AVX We performed the same study as in the previous paragraph, changing only the target architecture to use an Intel E52697 V2 processor and the results are reported on Table2.6.

		S 4		S 8			S 12			S 16
	n x	T scalar T avx	(×)	T scalar T avx	(×)	T scalar T avx	(×)	T scalar T avx	(×)
		(ms) (ms)		(ms) (ms)		(ms)	(ms)		(ms)	(ms)
	1	0.038 0.043	0.8	0.027 0.032 0.84	0.027 0.032 0.84	0.028 0.032 0.87
	2	0.031 0.035	0.8	0.032 0.034 0.94	0.037 0.035 1.05	0.041 0.035 1.17
	4	0.043 0.043	1.0	0.066 0.045 1.46	0.105 0.049 2.14	0.140 0.059 2.37
	8	0.144 0.108 1.33	0.342 0.140 2.44	0.638 0.171 3.73	0.934 0.239 3.90
	16	1.062 0.761 1.39	2.635 0.985 2.67	5.004 1.243 4.02	7.362 1.779 4.13
	32	7.335 4.806 1.52	19.90 6.671 2.98	38.76 8.684 4.46	57.67 13.13 4.39
	64	59.28 37.00 1.62	157.3 51.85 3.03	308.3 68.09 4.52	458.5 105.9 4.32

.4). Indeed, by vectorizing the computations, the number of integer calculations required for loop indexing, for instance, are reduced in the same time, which can then lead to the observed superlinear speed-up.

Table 2 .

 2 6: Illustration of the speed-up brought by the angular vectorization of the sweep, as a function of the MacroCell size (n x = n y = n z ), and for different angular quadrature orders. The performance measurements were carried-out on an Intel Xeon E52697 V2 processor.

Table 3 .

 3

	Adams et al. Hybrid Hybrid-Async

1: Notations used in the performance models of the sweep operation.

Table 3 .

 3 .2). Experiments are carried out using settings for ivanoe platform.

	small	big

2: Characteristics of the test cases used. GFlops is the number of floating point operations required for one complete sweep. We count 25 floating point operations per spatial cell per angular direction (see section 1.3.3).

Table 3 .

 3 3: Optimal domain partitioning and aggregation factors for the Adams et al. model using small test case on the ivanoe platform.

				2 40 60 15	384	126.5 84.3 337.5
	24	6	2	2 20 60 15	768	126.5 42.2 168.7
	48	6	4	2 20 30 10	2304	42.2 28.1 84.3
	96	8	6	2 15 20 10	4608	28.1 21.1 42.2
	192	12 8	2 10 15 10	9216	21.1 14.0 21.1
	384	24 8	2	5	15	5	36864 10.5	3.5	10.5
	576	24 12 2	5	10	5	55296	7.0	3.5	7.0
	768	24 8	4	5	15	5	36864 10.5	3.5	10.5
		Partitioning Aggregation	N tasks Message size (KB)
	N cores P Q R A x A y A z		s x	s y	s z
	12	1 1	1	5	5	5 110592 3.5 3.5	3.5
	24	2 1	1	5	5	5 110592 3.5 3.5	3.5
	48	2 2	1	5	5	5 110592 3.5 3.5	3.5
	96	2 2	2	6	5	5	92160 3.5 4.2	4.2
	192	4 2	2	10 10	5	27648 7.0 7.0	14.0
	384	4 4	2	10 10	5	27648 7.0 7.0	14.0
	576	6 4	2	10 10	5	27648 7.0 7.0	14.0
	768	8 4	2	15 10	5	18432 7.0 10.5	21.1

Table 3 .

 3 4: Optimal domain partitioning and aggregation factors for the Hybrid model using small test case on the ivanoe platform.predicted performances from Hybrid-Async simulator is 1.41 times higher than Hybrid model, and 1.05 times than Adams et al. model. This is due to the fact that the asynchronous approach is more suitable to leverage more parallelism compared to Adams et al. and Hybrid models.

Table 4 .

 4 .1). The scheduling and communication overheads are obtained from a fit of the Hybrid-Async model to real measurements. The network performance is given by a NetPIPE benchmark on the target computer. Figure4.3 presents the network latency of the ivanoe and athos computers.

	small	big

1: Characteristics of the test cases used. GFlops is the number of floating point operations required for a single complete sweep. We count 25 floating point operations per spatial cell per angular direction (see section 1.3.3).

Table 5 .

 5 1: One-group cross-sections and source strength S for the kobayashi benchmark.

Table 5 .

 5

		kobayashi	takeda		PWR
		(Problem 1)	(Model 1)
	N G	1	2	2	8	26
	N x	100	120		2 × 289
	N y	100	120		2 × 289
	N z	100	120	9 × 84	2 × 84	2 × 70
	N dir	288	80	168	80	288
	N dof	864.0 ×10 6	829.44 ×10 6	254.58 ×10 9	107.76 ×10 9	1.05 ×10 12
	Flops	7.2 ×10 9	6.91 ×10 9	2.12 ×10 12	898.0 ×10 9	8.75 ×10 12
	A x	20	10		20
	A y	20	10		20
	A z	20	10		20
	ϵ k eff ϵ ψ	--	10 -6 10 -5	10 -6 10 -6	10 -6 10 -5	10 -5 10 -5
	I g	1	1	1	5	4

[START_REF] Michael D Mccool | Structured Parallel Programming: Patterns for Efficient Computation[END_REF]

for the 26-group (resp. 2-group and 8-group) along the z axis. The spatial mesh is then refined by 2 × 2 × 2 for the 8-group and 26-group, and by 2 × 2 × 9 for the 2-group. The larger spatial mesh for the 2-group Figure 5.7: Radial view of a PWR 900 MW model [25]. 2: Description of benchmarks and calculation parameters.

Table 5 .

 5 4: Discrepancies on the scalar flux, between Domino and reference solution extracted from

			|δϕ 2 | (%)
	Core	0.44	0.16
	Control rod	0.33	0.20
	Reflector	0.01	0.24

Table 5 .

 5 5: Eigenvalue and computation time for the Model 1 of the takeda benchmark, on a 24-cores computing node of the athos platform. For this experiment, we fixed the number of SP N iterations to N SPN = 4.

	eff	k eff	N outer	Time (s)
		(Ref. Monte Carlo)	w/ PDSA	w/o PDSA	w/ PDSA	w/o PDSA
	0.962437	0.9624 ± 0.0005	14	281	2.43	27.46

Table 5 .

 5 ). 6: Solution times for a S 12 2-group 3D PWR k eff computation on the athos platform.

		Partitioning	(2, 2, 1)	(2, 2, 2)	(4, 2, 2)	(4, 4, 2)
		N cores		96	192	384	768
	ρ max PDSA			0.21	0.21	0.21	0.21
		N outer		2116	2120	2116	2125
	w/o PDSA	T sweep	(s)	6270.2	3831.0	2609.3	2289.8
		T total	(s)	7724.4	4543.6	2974.2	2479.0
		% sweep		81.1	84.3	87.7	92.3
		N outer		92	95	97	81
		T sweep	(s)	279.9	144.6	77.0	41.0
	w/ PDSA	T spn	(s)	47.8	24.6	12.4	5.1
		T comm PDSA	(s)	35.7	30.4	15.6	7.9
		T total	(s)	451.6	245.7	129.2	65.0
		% sweep		61.9	58.8	59.5	63.0
	Perf. (Tflop/s)	sweep Domino		2.7 0.8	3.5 1.6	5.3 3.2	8.3 5.3
	speed-up		17.1	18.5	23.0	45.7

Table 5 . 7 :

 57 s of which 91.53 s comes from the sweep operation, Domino-Denovo comparison for a S 12 2-group 3D PWR k eff computation.

	Denovo [29]	Domino
	PWR 2g	PWR 2g

(P, Q, R) ρ max

Table 5 .

 5 8: Solution times for a S 12 8-group 3D PWR k eff computation. illustrating that the sweep operation is still dominant (71% of the total time).

26-group PWR core model

Table 5

 5 .9. The convergence is reached in 126 outer iterations, for a global solver (P, Q, R) ρ max

		PDSA	N outer	T sweep (s)	T spn (s)	T comm PDSA (s)	T total (s)
	(4, 4, 4)	0.14	126	2226.42	56.56	147.2	2763.52

Table 5 .

 5 9: Solution times for a S 16 26-group 3D PWR k eff computation. time of 2763.52 s (46 min

HiePACS: High-End Parallel Algorithms for Challenging Numerical Simulations, https://team.inria.fr/ hiepacs/.

330 • C in the case of PWR 900 MW. The coolant is maintained in a liquid state thanks to a high pressure (10 MPa) maintained in the core by the pressurizer.

This process is known as thermalization or moderation of the neutrons.

There is an up-scattering if a neutron of energy g is scattered into a neutron of g ′ > g.

Thermal groups are those having lowest energy. Fast groups are those having highest energy.

The number of transistors per chip doubles approximately every two years

Power density is the amount of power (time rate of energy transfer) per unit volume. https://en.wikipedia. org/wiki/Power_density

http://www.gotw.ca/publications/concurrency-ddj.htm

A full list of latency and throughput of instructions is available at http://www.agner.org/optimize/ instruction_tables.pdf

According to an Intel communication at https://gcc.gnu.org/wiki/cauldron2014?action=AttachFile&do= view&target=Cauldron14_AVX-512_Vector_ISA_Kirill_Yukhin_20140711.pdf, processors that will implement the avx-512 ISA are: Intel Knights Landing co-processor and processors of Xeon series.

It should be noted that Stream does not account for caches. It targets DRAM bandwidth.

https://code.google.com/p/likwid/

http://icl.cs.utk.edu/papi/

http://research.microsoft.com/en-us/um/people/minka/software/lightspeed/

It is expected that the penalty cost for an unaligned load/store operation is going to decrease for next generations of processors.

A comprehensive study on processor microarchitecture is available from http://agner.org/optimize/ microarchitecture.pdf.

http://www.lanl.gov/asc/

http://wwwc3.lanl.gov/pal/software/sweep3d/sweep3d_readme.html

Partisn is a Cartesian S N neutron transport code developed at LANL.

Note here that, as previously said, we do not use the full number of cores per node for task computations; rather a single core is entirely dedicated for communications.

A dependency D of a task T is a task that has to be processed before the processing of T.

http://www.mpi-forum.org/

In the future, we will add another level of parallelism by the mean virtual processes to take advantage of "highly" NUMA nodes.

Let a 1D domain composed of n regions, of different properties; Σ i , li, for i = 1,

2, . . . , n being respectively the total cross-section and length of each region. Then the optical thickness of the whole domain is defined byτ = Σ 1 l1 + Σ 2 l2 + • • • + Σ n ln.

Remerciements

https://software.intel.com/en-us/code-modernization-enablement 2

// ncx , ncy and ncz : number of MacroCells along the x , y and z axes const int nfront = ncx + ncy + ncz -2;

// An octant o is defined by ( fx , fy , fz ) such that o = fx +2* fy +4* fz for ( int f =0; f < nfront ; ++ f ){ for ( int z =0; z < min ( ncz , f +1); ++ z ){ for ( int y = max (0 ,f -ncz ); y < min ( ncy , f +1 -z ); ++ y ){ const int x =f -y -z ;

, STARPU_RW , mcg -> getPsi ( mcg , 0 , fx , fy , fz , 0 , y , z ) , STARPU_RW , mcg -> getPsi ( mcg , 1 , fx , fy , fz , x , 0 , z ) , STARPU_RW , mcg -> getPsi ( mcg , 2 , fx , fy , fz , x , y , 0 ) , STARPU_EXECUTE_ON_NODE , mcg -> rank_of ( mcg , x , y , z ) , 0 );

Experimental platforms

We consider three computing platforms based on x86 multicore processors, as presented on Table A bigmem is a NUMA node featuring four octo-core Intel Xeon E7-8837 processors running at 2.66 GHz. Each CPU core on each of these processors supports Intel SSE 1 extensions of the x86 Instruction Set Architecture (ISA), allowing to perform 4 (resp. 2) floating point arithmetic operations (add/mul/sub) in single (resp. double) precision, in a single clock. The theoretical peak performance of this computing node is 680.96 GFlop/s in single precision (SP).

ivanoe is a distributed memory computer. Each computing node of this cluster is a NUMA node featuring two hexa-core Intel Xeon X5670 processors running at 2.93 GHz. Each of the CPU cores supports Intel SSE extensions, as for the bigmem computing node. Each computing has a theoretical peak performance of 281.28 GFlop/s in single precision. The network topology of this cluster is a fat-tree switch fabric of QDR 2 InfiniBand TM (IB) type. This cluster is 1 Streaming SIMD Extensions (SSE) is an SIMD instruction set extension to the x86 architecture, designed by Intel and introduced in 1999 in their Pentium III series processors (source https://en.wikipedia.org/wiki/ Streaming_SIMD_Extensions).

2 Quad Data Rate (QDR) is a communication signaling technique wherein data are transmitted at four points 111 equipped with a Mellanox TM IB card whose theoretical effective throughput is 40 Gbits/s. athos is another distributed memory computer, whose network topology is a fat-tree switch fabric of FDR 1 InfiniBand TM type, and equipped with a Mellanox TM IB card, achieving a theoretical effective throughput of 56 Gbits/s. Each computing node of this cluster is a NUMA node featuring two twelve-core Intel Xeon E5-2697 v2 processors. Each CPU core supports Intel AVX 2 extensions of the x86 ISA, allowing to perform 8 (resp. 4) floating point operations in single (resp. double) precision, in a single clock. The theoretical peak performance of each computing node is 1036.8 GFlop/s in single precision.

in the clock cycle (source https://en.wikipedia.org/wiki/Quad_data_rate).

1 Fourteen Data Rate (FDR)