Keywords: heterogeneous architectures for Monte Carlo neutron transport codes Manycore architectures, Heterogeneous architectures, Particles transport 3P Portability, Performance and Productivity. 102 COMA Cache-Only Memory Access Model. 12 DRAM Dynamic Random-Access Memory. 15 DRAM Static Random-Access Memory. 15 EPP Explicit Productivity Portability. 122, 129 HB History-based Method. 77 IPP Implicit Productivity Portability. 121

Titre : Évaluation de mod èles de programmation pour les architectures manycore et / ou h ét érog ènes pour les codes de transport neutronique Monte Carlo Mots cl és : Architectures manycore, Architectures h ét érog ènes, Transport des particules R ésum é:Dans cette th èse nous nous proposons d' évaluer les diff érents mod èles de programmation disponibles pour adresser les architectures de type manycore et / ou h ét érog ènes dans le cadre des codes de transport Monte Carlo. On consid èrera dans un premier temps un cas test d'application simple mais repr ésentatif pour couvrir un éventail assez large de solutions et les comparer en terme de performance, de portabilit é de la performance, de facilit é de mise en oeuvre et de maintenabilit é. Les architectures cibles sont les CPU 'classique', Intel Xeon Phi et GPU. Les mod èles de programmation les plus pertinents seront ensuite mis en place dans un code de transport Monte Carlo.

Evaluation of performance portability for different implementations of slabAllNuclides across a set of platforms via

the generic model F and the absolute performance metric, particle tracking rate, the higher the better. . . . . . . . 5.6 Evaluation of performance portability for different implementations of slabAllNuclides across a set of platforms via the generic model F and the relative performance metric, application efficiency, the higher the better. . . . . . . . .

Introduction 1.1 Motivations

Monte Carlo neutron transport simulation is a stochastic method which is widely used in the nuclear field for reference calculations. Instead of introducing mathematical and physical approximations to simulate a real-world system, Monte Carlo method aggregates particles behaviors through statistical sampling which significantly leads to a huge amount of computational cost. In order to alleviate this cost, the use of supercomputers for Monte Carlo simulation has become a trend for many researchers and developers.

However, as more and more modern architectures (multi-core/manycore, heterogeneous architectures) emerge in a rather fast rhythm, it is non-trivial to develop full portable applications on all these architectures, furthermore, to optimize and maintain the application continuously so as to achieve preferable performances on them.

From one side, the development and life time of an application are much greater than the time for the development of computing architectures. Thus, it makes little sense to design an application for a given architecture since after the development of the application there is a good chance that the architecture has already evolved.

From the other side, it is necessary for applications to be scalable and maintainable by as many people as possible, allowing those people who are not familiar with HPC to develop the codes capable of delivering good performance.

Therefore, developing applications of high performance portability becomes a hotspot to address in the field of HPC.

To meet the challenges mentioned above and gain better performance portability, some high-level programming languages such as OpenACC, OpenMP have been proposed to allow the application performing well on a wide range of architectures without large development and maintenance effort. Most of these programming models is dedicated to accelerator-based architectures with the massive increase of multi-level parallelism which represents the underlying trend in computational architectures. In this thesis, we are interested in the development of portable Monte Carlo neutron transport code targeting to hybrid architectures (CPU + GPU) with the utilization of hybrid programming models.

Chapter 2

Parallel Computing

The thesis work is based on the development of parallel code for Monte Carlo neutron transport simulation. Therefore, we begin with the literature of thesis by giving an introduction of parallel computing from both the hardware and software perspectives. It mainly consists of the basic information about parallel computing and some hardware architectures, programming models as well as profiling tools that are used in research work.

General Terminology

Parallel computing is a type of computation in which many calculations or the execution of processes are carried out simultaneously [START_REF] Gottlieb | Highly Parallel Computing[END_REF]. Instead of dividing a problem into a group of instructions and executing them sequentially, parallel computing firstly breaks the problem into several sub-problems and carries out them concurrently.

Nowadays, the demand of parallel computing has reached to a very high level. With the blowout of big data, numerical simulation needs a great amount of processing time which is far beyond the capacity of any single compute resource. Thus, more and more parallel computing are used in a wide range of scientific and industrial domains for the purpose of exploring the potential of underlying parallel hardware, making better use of dispersed compute resource, and eventually saving the cost of time and money. In what follows some general parallel terminology will be introduced to make better acknowledge of parallel computing.

Von Neumann Architecture

The von Neumann architecture is the fundamental design of modern parallel computers proposed by the Hungarian mathematician and physicist John von Neumann in his 1945 paper [START_REF] Neumann | First Draft of a Report on the EDVAC[END_REF]. This stored-program design allows instruction and data being stored in electronic memory. Its main components are illustrated in Figure 2.1.

• Central Processing Unit (CPU): An electronic circuitry which executes instructions of a computer program.

It contains a processing unit and a control unit.

⇤ Processing Unit (PU): It is comprised of an arithmetic logic unit (ALU) and a variety of registers. ALU allows for arithmetic and bitwise operations on input data and returns the result as output. Registers are fast storage areas which provide the fasted way to access data.

⇤ Control Unit (CU): A component of CPU that fetches and decodes instructions and directs the operations of the processor. It is comprised of a current instruction register (CIR) and a program counter (PC). The CIR is responsible for holding the current instruction during processing. The PC contains the address of the next instruction to be executed and increments itself after fetching an instruction.

• Memory Unit: It consists of random access memory (RAM) which stores data and instructions.

• Input/Output: Input and output mechanisms. One of the key features of the von Neumann architecture is that both data and instructions are designed to share common bus to communicate between central processing unit and memory unit. Thus, data operation and instruction fetch cannot be executed simultaneously. This drawback is referred to as the von Neumann bottleneck and usually limits the performance of the system significantly [START_REF] Markgraf | The Von Neumann Bottleneck[END_REF]. Two strategies are mostly adopted by modern parallel computers to mitigate this bottleneck at memory hierarchy level, which are respectively adding internal memory layer such as cache to communicate with central processing unit and providing separate storage and pathways for data and instructions.

Flynn's Taxonomy

Flynn's Taxonomy is classification of parallel computing machines proposed by Michael J. Flynn in 1966 [START_REF] Flynn | Some Computer Organizations and Their Effectiveness[END_REF]. This classification is based on distinguishing the multiplicity (single/multiple) of instruction and data streams. As a result, there are four sub-categories according to Flynn's taxonomy:

• SISD: single instruction stream, single data stream

• SIMD: single instruction stream, multiple data stream

• MISD: multiple instruction stream, single data stream

• MIMD: multiple instruction stream, multiple data stream

SIMD organization

SIMD organization allows multiple processing units handling the same instruction broadcast from one control unit.

Meanwhile, all processing units take different data streams to accomplish arithmetic operations simultaneously. The details of SIMD organization is depicted in Figure 2 Some parallel computers designed via SIMD organization are IBM 9000, Cray X-MP, ILLIAC-IV and Thinking Machines CM-2. Furthermore, SIMD organization is widely employed in most parallel computers as specific instructions and execution units to achieve SIMD vector operation such as vector processing unit (VPU) and graphics processing unit (GPU). MISD machine operates multiple instruction streams and single data stream with multiple control units and processing units. Each control unit handles one instruction stream and sends it to the corresponding processing unit. All processing units share a single data stream to complete execution as shown in Figure 2.5. There are few examples of MISD machines that have ever existed.

MISD organization

MIMD organization

MIMD organization is the most popular for a parallel computing machine. Most modern supercomputers, networked grids and all multiprocessor systems fall into this classification, such as IBM POWER5, Intel IA32, Cray XT3, etc. In contrast to MISD organization, in which all processing units are fed with only a single data stream, each processing unit in MIMD organization can handle a different data stream. Since there is no data and instruction dependencies among processors, executions on different processors can be asynchronous. The diagram of MIMD organization is denoted in Figure 2 Furthermore, MIMD contains a subcategory called as SPMD (single program, multiple data) [START_REF] Darema | The SPMD Model: Past, Present and Future[END_REF]. It allows a single program to be executed simultaneously by multiple processors with different input, which usually refers message passing programming (section 2.4.1.1) on distributed memory architectures (section 2.2.2).

Parallel Speedup

Amdahl's law

Amdahl's law [START_REF] Amdahl | Validity of the Single Processor Approach to Achieving Large Scale Computing Capabilities[END_REF] is a formula indicating the potential speedup of an algorithm which can be expected by parallel computing, as illustrated in Equation 2.1:

Speedup(P, N)= 1 P N +1 P (2.1)
where:

• P : The proportion of run time of a given algorithm that can be parallelized.

• N : The number of processors performing the program.

• P 2 [0, 1]:

⇤ if P =0, then Speedup =1, meaning that the algorithm has no parallel fraction and cannot be parallelized to achieve performance speedup.

⇤ if P =1 , then Speedup = N , signifying that the optimal speedup of algorithm is N. With the increase of number of processors, the speedup can be infinite in theory.

• N 2 [1, +1):

⇤ if N =1, then Speedup =1, the code is executed in serial mode.

⇤ if N ! +1, then Speedup = 1 1 P
. This is a derived version of Amdahl's law, assuming that the execution time of parallel fraction can be reduced to zero with multiple processors. This theoretical achievable speedup is considered as the upper performance limit of parallelization of the given algorithm.

Amdahl's law assumes the parallelizable part of an algorithm P is a constant. However, this assumption does not always fit into real-world problems. In practice, the parallel fraction P often varies with the change of problem size. Considering of this factor, we can have another derived version of Amdahl's law as shown in Equation 2.2:

Speedup(P (W ),N)= 1 P (W ) N +1 P (W ) (2.2)
where:

• W : The workload of a given problem.

• P (W ): The parallelizable proportion of execution time of a given problem with workload W .

Unfortunately, P (W ) is not a predictable function varying via specific rules. To overcome this shortcome, we can consider in another way. Instead of fixing the problem size and comparing the execution time, it is possible to firstly define a fixed execution time and then to calculate the speedup by comparing the achieved workload of program in parallel and serial modes.

Gustafson's law

Gustafson's law was proposed by John L. Gustafson and his colleague Edwin H. Barsis in the article "Reevaluating Amdahl's Law" in 1988 [START_REF] Gustafson | Reevaluating Amdahl's Law[END_REF]. As shown in Equation 2.3, Gustafson's law assumes that the total parallel workload varies linearly with the number of processors. It estimates the theoretical attainable speedup of a program by comparing the workload of the parallel version of the code W p with N processors to the workload of the serial version of the code W s at a fixed execution time T .

W s =(1 P )W s + PW s ;

W p =(1 P )W s + NPW s ; Speedup(N )= W p W s =1 P + NP; (2.3)
where:

• W s : The workload of the code in serial mode at a fixed execution time T .

• W p : The workload of the code in parallel mode with N processors at a fixed execution time T .

• P : The parallelizable fraction of the code at a fixed execution time T .

Scalability

In the context of parallel computing, scalability refers to the property of a system to handle a growing amount of work by adding resources to the system [START_REF] Bondi | Characteristics of Scalability and their Impact on Performance[END_REF]. There are two common types of scaling:

• Strong scaling: The strong scaling indicates the scalability tested under the condition that the time to solution varies with the number of processors for a fixed total problem size.

• Weak scaling: The weak scaling describes the scalability evaluated under the condition that the time to solution varies with the number of processors for a fixed problem size per processor.

There are many factors which influence the scalability of a parallel program. For example, the algorithm used in the program often contains some inherent limits which lead to performance degradation with the increase of compute resources. The memory capacity, memory bandwidth or network bandwidth of a parallel computing machine may also be the bottleneck of scalability.

Memory Architectures

Parallel computing machines can be separated into three types in the context of memory: shared memory architecture, distributed memory architecture and hybrid memory architecture. Processors in shared memory architecture share the global memory modules while processors in distributed memory architecture have their own local memory and communicate with each other via a network. Hybrid memory architecture is the mix of shared and distributed memory architectures. These memory architectures have their own advantages and disadvantages, we will dig into them in the following part.

Shared Memory

Shared memory is a global address space which can be accessed by all processors. The key advantage of shared memory architecture is that multiple processors have fast access to data in memory location, leading to high memory bandwidth. Meanwhile, such fast access to memory incurs two disadvantages. First, the resource contention and false sharing usually happen when several processors try to access the same or nearby memory locations. When adding more processors to shared memory architecture, the performance degradation will become more and more serious. Second, since several processors may handle data load/store operations simultaneously, data coherence must be kept to ensure that different processors always use the "correct" data that has just been updated by other processor. Such obligation for data coherence is quite a bottleneck to performance. The shared memory architecture can further be divided into three models which are based on different manners of memory access, as listed below:

• Uniform Memory Access Model (UMA): All processors in UMA are identical. They have equal access and equal access time to shared memory, as shown in Figure 2.7 (a).

• Non Uniform Memory Access Model (NUMA): Global memory is not accessed uniformly by all processors.

Some processors attached to the same local memory have equal access time. However, the access to a non-local memory across link for these processors is slower. The NUMA model is described in Figure 2.7 (b).

• Cache-Only Memory Access Model (COMA): As illustrated in Figure 2.7 (c), COMA model is a specific NUMA model where its local memories are used as cache (attraction memory) instead of main memory.

The difference between main memory and attraction memory is that when a data is accessed by non-local processors, main memory keeps the space of data and gives a copy of data to processors to operate. On the contrary, attraction memory migrates the data to wherever it is needed.

Distributed Memory

In distributed memory architecture, processors have their own local memory instead of sharing a global address space. The separate memories are connected by message passing interconnection network such as Ethernet.

When a processor needs to access remote data in another processor, it uses the network link to communicate with the remote processor. The diagram is depicted in Figure 2.8: Generally speaking, the key advantage of distributed memory architecture is that the problems of memory conflicts (for example memory contention, race condition, cache coherence) can be alleviated since each processor has its own private memory. However, the data access to remote memory consumes longer time than the access to local memory, which may degrade significantly the performance if it happens too much inter-communications during the execution of a program. Another shortcome of distributed memory is that programmers need to explicitly define the mechanism of communications and synchronizations of processors for a given code, which augments much programming complexity.

Hybrid Memory

Hybrid memory architecture, also called distributed-shared memory architecture, is widely employed by modern high performance computers. As Figure 2.9 shows, hybrid memory architecture consists of a group of shared memory machines which are connected to each other via inter-connection network. A parallel computing machine designed in hybrid memory architecture scales well with the increase of shared memory components, which makes it easier to exploit the computing resources at the cost of increased programming complexity.

Figure 2.9: Hybrid Memory Architecture.

Computing Architectures

The current computer market is full of different hardware architectures produced by many vendors such as INTEL, NVIDIA. These architectures are designed for solving different problems and each product has its advantages and defects. In order to explore maximal computing capacity of a targeting architecture and learn about the performance bottlenecks on this architecture for a given code, it is required to get better knowledge of the targeting architecture.

Therefore, we will express in detail several computing architectures used in the thesis work.

Hardware Characteristics

We start from highlighting some features which are widely used to characterize hardware capacity.

• Latency: refers to the time it takes for the execution of an operation. Its conventional unit of measurement is microseconds.

• Frequency: indicates the number of electronic pulses generated by a processor to synchronize operations of its components per second. It can also be called "clock rate" with the units of measurement MHz or GHz.

• Bandwidth: signifies the amount of data which can be handled per unit of time. gigabytes/sec, megabytes/sec, or bytes/cycle are widely used as units of measurement.

• Throughput: denotes the amount of floating-point operations which can be processed per unit of time. GFLOPS is a common unit of measurement.

In a word, the evolution of computing architectures is mainly to design a hardware with lower latency, higher frequency, higher bandwidth and higher throughput.

CPU Architectures

.10: Modern CPU Architecture [90].

CPU is the electronic circuitry within a computer that executes a computer program by performing the basic arithmetic, logical, control and input/output(IO) operations specified by instructions [START_REF] Weik | A Third Survey of Domestic Electronic Digital Computing Systems[END_REF]. It is a general purpose processor designed to perform a bunch of operations integrating the arithmetic and control operations while its performance of these tasks is not good enough. Unlike the traditional structure of CPU as we illustrated in section 2.1.1, where its components refer to processing unit and control unit, modern CPUs include memory units such as cache and DRAM as depicted in Figure 2.10:

• Cache: An intermediate memory layer made of SRAM which has higher bandwidth and lower latency compared with main memory.

• DRAM: Dynamic random-access memory is primarily used as main memory. It is comprised of an array of memory cells in which a transistor-capacitor circuit is used to store one bit of data based on its state of charge.

Due to capacitor leakage, DRAM needs to perform refresh process, which is a key difference to SRAM and makes DRAM much slower than SRAM.

• DRAM: Static random-access memory uses multiple transistors in a latching circuit to store each bit of data.

Since it makes no use of transistors, SRAM does not need to be refreshed and it has higher efficiency on data transfer comparing to DRAM.

Intel Broadwell

.11: Broadwell Die Diagram [4].

Intel Broadwell is the fifth model generation of Intel Processor released in 2015 with 14nm die shrink of Haswell microarchitecture [START_REF] Nalamalpu | Broadwell : A family of IA 14nm processors[END_REF]. Based on 14nm process technology, Broadwell reduces 49% feature neutral area and gains better power efficiency due to improved Vmin (minimum operating voltage). Comparing to Haswell microarchitecture, Broadwell introduces some new instruction set extensions including RDSEED, PREFETCHW, SMAP and Intel ADX. There are several instructions that are improved for lower latency in Broadwell such as FP MUL instructions and CLMUL instructions. Broadwell also enlarges instruction queue, scheduler and second-level TLB for higher throughput.

Intel Skylake

Skylake is the 6 th generation of Intel microarchitecture succeeding the Broadwell microarchitecture with 14nm process (client) or 14nm+ process (server) technologies. The modified 14nm+ process has promoted transistor channel strain which leads to higher frequency of processor (100 ⇠ 300 MHz) and better performance. Unlike Skylake (client) which adopts ring topology on-die interconnect firstly introduced in Sandy Bridge architecture [START_REF] Yuffe | A Fully Integrated Multi-CPU, GPU and Memory Controller 32nm Processor[END_REF], Skylake (server)

implements a new networking topology called mesh interconnect architecture that mitigates the scaling bottlenecks of latency and bandwidth as the number of cores increases in ring interconnect architecture. These two types of Skylake machines are both used in the thesis work.

Note that Skylake introduces a number of new instructions in which the advanced vector extension using 512bit SIMD techniques (AVX-512) [START_REF] Reinders | AVX-512 Instructions[END_REF] is one of the most important feature contributing to parallel computing. The AVX-512 instruction set in Skylake is composed of five separate sets:

• AVX512F: The fundamental of the 512-bit SIMD instruction extensions which is required for any implementations of AVX-512 extensions on microarchitectures. It is the extensions to AVX/AVX2 instruction sets using EVEX prefix so as to support 512-bit registers, arithmetic operations, bitwise operations, data management, parameter broadcasting, etc.

• AVX512CD: AVX512CD instruction set is designed specifically for conflict detection. Loops without address conflict can be detected and then vectorized to improve performance.

• AVX512BW: BW refers to bytes (8-bit) and words (16-bit). This instruction set extends AVX-512 extensions to cover 8-bit and 16-bit integer operations [START_REF] Reinders | Additional AVX-512 Instructions[END_REF].

• AVX512DQ: DQ represents doubleword (32-bit) and quadword (64-bit). AVX512DQ instruction set introduces some new 32-bit and 64-bit operations such as conversion and transcendental support.

• AVX512VL: One of an additional AVVX-512 instruction set dedicated to the compatibility of AVX-512 operations on registers with different vector length. For instance, XMM (128-bit SSE) registers and YMM (256-bit AVX) registers.

IBM POWER9

IBM POWER9 is a 14nm microarchitecture based on Power ISA (instruction set architecture) v3.0 [START_REF]Power ISA Version 3.0B[END_REF] firstly released in 2017. The key component of POWER9 architecture is the 4-way or 8-way simultaneous multithreading core, as SMT4 core and SMT8 core. Each super-slice combines 2 slices together where each slice consists of a 64-bit vector and scalar unit (VSU) coupled with load/store unit (LSU). Each slice is able to handle an address generation (AGEN) operation and a 64bit scalar operation including either a computational operation (fixed-point/floating-point) or a load/store operation at any given cycle. A super slice can be considered as a vector processing unit issuing two AGEN operations and a 128-bit vector operation every cycle. Besides, a POWER9 SMT core also contains 4 double-precision floating-point units (DFUs), which allows for issuing double-precision floating-point multi-add/divide/square root instructions every cycle. This design increases POWER9 CPU's double-precision throughput.

The POWER9 processor has a 10MB eDRAM-based L3 cache chunk per SMT8 core (or shared between 2 SMT4 cores), as denoted in Figure 2.13. Different L3 cache chunks are linked by fabric interconnect to allow cores sharing one L3 cache chunk accessing data residing in other L3 chunks. The POWER9 L2 cache is similar to L3 cache as each 512KB L2 cache is associated with a SMT8 core or two SMT4 cores. Moreover, the POWER9 processor has a 32KB L1 instruction cache in conjunction with a 32KB L1 data cache per SMT4 core.

GPU Architectures

GPU (graphics processing unit) is originally a specialized electronic circuit dedicated to rapidly handle image processing computations for output to a display device. Back then, it was widely used in personal computers, workstations or game consoles, only dealing with graphical operations. However, with more and more emergence of programmable GPUs and programming interfaces, GPU's highly parallel architecture has been found more performance-efficient than CPU architecture performing compute-intensive and data parallel calculations. Thus, a new software concept called "general-purpose computing on graphics processing unit" has been proposed and largely developed which allows GPUs handling general-purpose applications in different domains such as geometric computing, bioinformatics, computational finance, scientific computing and so on [START_REF] Chitty | A Data Parallel Approach to Genetic Programming using Programmable Graphics Hardware[END_REF].

As denoted in Figure 2.14, GPU contains thousands of GPU cores settled massively in parallel so as to hide latency and maximize throughput. Comparing GPU core with CPU core, GPU core is designed for data-parallel computations with simple control logic whereas CPU core is rather heavy-weight, aimed at handling efficiently complex control logic such as an algorithm composed of multiple conditional operations (if...else).

Regarding to GPU caches, state-of-the-art GPUs are being designed with sizable hardware-managed multilevel caches, such as L1 cache, L2 cache and texture cache [START_REF]A Survey of Techniques for Managingand Leveraging Caches in GPUs[END_REF]. The cache memory has been enlarged with the evolution of modern GPUs. For example, the Nvidia Tesla Volta GPU [START_REF] Nvidia | NVIDIA TESLA V100 GPU ARCHITECTURE[END_REF] has 6MB last level cache while the latest released Nvidia Ampere Tensore Core GPU [START_REF] Nvidia | NVIDIA A100 Tensor Core GPU Architecture[END_REF] has 40MB last level cache. Note that the last level cache size of Nvidia Ampere architecture is competitive with that of Intel Skylake architecture containing up to around 30 cores. But if we consider the average cache memory which can be used for each core, Nvidia Ampere has 5KB last level cache per core and that of Intel Skylake is 1.375MB per core. We can say that GPU is composed of tens of hundreds of cores with relatively small caches. On the contrary, CPU is composed of tens of cores with relatively large caches.

Furthermore, modern GPUs consume significant amount of power as the cost of performance optimization.

For example, the Nvidia Tesla Pascal GPU [START_REF] Nvidia | [END_REF] has a TDP (thermal design power) level of 300W for maximum Figure 2.14: Modern GPU Architecture [90], where "green square" refers to GPU core, "yellow rectangle" is GPU control unit, "orange rectangle" signifies GPU cache or shared memory.

performance. In contrast, Intel Xeon Gold 6148 has a maximum power consumption of 150W. Since more and more

GPUs have been widely used in high-performance computing (HPC) machines, the optimization of GPU power efficiency has become a research hotspot to make supercomputers in future more friendly in economy and energy.

There are some research work that has been carried out to optimize power efficiency of GPUs from a hardware perspective, such as saving static energy in L1 and L2 caches [START_REF] Wang | Run-time Power-gating in Caches of GPUs for Leakage Energy Savings[END_REF], using filter-cache to reduce accesses to instruction cache [START_REF] Lashgar | Inter-warp Instruction Temporal Locality in Deep-Multithreaded GPUs[END_REF] and so on. From a developer perspective, enhancing power efficiency of GPUs for a specific application is mainly to optimize codes so as to make maximal use of computing resources.

In the following part of this section, we will introduce several GPU microarchitectures related to our research work.

Intel Gen9.5

Intel Gen9.5 is the 14nm+ microarchitecture of Intel's integrated graphics processing unit (IGPU) used in several Intel's CPU microarchitecture such as Kaby Lake, Coffee Lake, etc. Unlike dedicated graphics cards which have their own RAM, IGPUs share the system RAM with the CPU.

Nvidia Pascal

Nvidia Pascal microarchitecture was firstly introduced in 2016 with the release of Tesla P100 card. As the successor to Nvidia Maxwell microarchitecture, GP100 architecture provides some improvements such as enhanced peak throughput (For Tesla P100, 5.3 TFLOPS in double-precision, 10.6 TFLOPS in single-precision). The number of FP32 (single-precision floating-point) cores and FP64 (double-precision floating-point) cores in GP100 architecture are respectively 3840 and 1920. Thus, the ratio of FP32 units to FP64 units for GP100 is 2. By comparison, the ratio of FP32 units to FP64 units in Maxwell architecture is 32. Furthermore, GP100 can handle FP16 (half-precision floating-point) operations which deliver great speedups for many deep learning algorithms.

With respect to memory, the overall memory size of GP100 is enlarged including register, L2 cache and shared memory. HBM2, as the second generation of High Bandwidth Memory, is firstly adopted in GP100 architecture which offers three times higher memory bandwidth compared with the Maxwell GM200 GPU, in which GDDR5 (graphics double data rate type 5) memory is used. GP100 also use a new high-speed interconnect technology called NVLink to increase performance for both CPU-to-GPU and GPU-to-GPU inter-communications as a substitution of PCI Express (PCIe) [START_REF] Nvlink | Pascal and Stacked Memory: Feeding the Appetite for Big Data[END_REF]. The NVLink implementation in GP100 (NVLink 1.0) delivers up to 160GB/s bidirectional bandwidth with maximum 4 links. Besides, GP100 adds two hardware features to improve unified memory functionality: 49-bit virtual addressing support and page faulting mechanism. Due to these features, GP100 unified memory offers a system-wide virtual address space and it can be accessed simultaneously by CPU and GPU.

A full GP100 is composed of 60 streaming multiprocessors (SMs), a 4096KB L2 cache, eight 512-bit memory controllers, four 4096-bit HBM2 stacks and 4 NVLinks. The number of SMs may vary for different GP100 products.

Chapter 2

Parallel Computing 21

For example, the Tesla P100 only enables 56 SMs. Each SM contains an instruction cache, a texture/L1 cache, 4 texture units, a 64KB shared memory and 2 processing blocks. Each processing block is comprised of 32 FP32 cores, 16 FP64 cores, an instruction buffer, a warp scheduler, 2 dispatch units and a 128KB register file, as illustrated in Figure 2.16. is unfriendly to handle thread-aware algorithms with data dependency among threads in a warp.

Nvidia Volta

Nvidia Volta is the microarchitecture succeeding Pascal architecture firstly released in 2017 with the Tesla V100 product. A full GV100 consists of 84 SMs, a 6144KB L2 cache, eight 512-bit memory controllers, four 4096-bit HBM2 stacks and 6 NVLinks. Note that the Tesla V100 uses 80 SMs.

As we can see in Figure 2.17 Figure 2.17: Nvidia Volta GV100 Streaming Multiprocessor Unit [START_REF] Nvidia | NVIDIA TESLA V100 GPU ARCHITECTURE[END_REF].

Comparing to GP100 architecture, GV100 has many advances from both hardware and software perspectives.

Since L1 data cache is merged with shared memory, L1 cache in GV100 benefits from the low latency and high bandwidth of shared memory. Therefore the performance gap between a code with and without the use of shared memory is significantly narrowed. GV100 supports 6 NVLinks 2.0 delivering in total 300 GB/s bidirectional bandwidth which is around 2⇥ higher than the interconnect bandwidth of GP100. The HBM2 memory used in GV100 can yield maximum 900 GB/s memory bandwidth which is higher than 700 GB/s in GP100. Besides, GV100 add a new "independent thread scheduling" feature to SIMT model, where each thread within a warp maintains its own program counter and call stack so as to keep an independent execution state. Due to this new feature, statements from if branch can be interleaved with statements from else branches and vice versa, making GV100 capable of carrying out thread-aware algorithms where data-dependent executions are required among threads in a warp. With respect to unified memory in GV100, it is improved by the introduction of a new Access Counter and the Address Translation Services (ATS) support over NVLink.

Heterogeneous Architecture

Heterogeneous architecture refers to a hybrid architecture combined through memory bus (PCIe, NVLink) with

CPUs as the host and some processing units as the devices, including GPUs, DSPs (digital signal processors) [START_REF] Frantz | Digital signal processor trends[END_REF],

FPGAs (field-programmable gate arrays) [START_REF] Cong | High-Level Synthesis for FPGAs: From Prototyping to Deployment[END_REF] and so on. A host is responsible for initializing the environment and launching the application. A device, also called as an accelerator, is dedicated to handle the compute-intensive parts of the application so as to enhance the performance of application. The calculations running on a heterogeneous architecture is called heterogeneous computing, which is a popular trend in HPC domain. The most widely used heterogeneous architecture is abstracted as CPUs + GPUs, including many different combinations such as Intel Broadwell CPUs + 2 Nvidia P100 and IBM Power9 CPUs + 4 Nvidia V100.

Parallel Programming Models

A parallel programming model refers to an abstraction of underlying hardware architecture that allows for the expression of both algorithms and data structures. It is different from programming languages with the invocation of an abstract execution model which cannot be understood by programming languages. In other words, a parallel programming model exists independently of the choice of the programming language. It bridges the gap between hardware and software.

There are mainly two classes of parallel programming models, which are individually called process interaction and problem decomposition. With regards to process interaction, it denotes that the parallel processes are able to communicate with each other. The most common forms of process interaction are shared-memory and messagepassing. There is also implicit interaction in which the process interaction is invisible to the programmer, such as domain-specific languages and functional programming languages.

As for problem decomposition, it specifies the way in which the child processes of a program are formulated.

Task parallelism is one of the problem decomposition which focuses on processes. Data parallelism describes a way of problem decomposition by performing operations in parallel on a data set. There is also implicit parallelism which reveals nothing to the programmer, such as automatic parallelization support in compilers.

Low-Level APIs

CUDA

CUDA is a general-purposed parallel computing platform and programming model created by NVIDIA [START_REF][END_REF] that leverages Nvidia GPUs to solve complex computational problems. It comes with a software environment that allows developers to use C, C++, Fortran, Python or other languages. This accessibility makes it easy for developers to manipulate Nvidia GPUs in parallel programming.

The basic CUDA processing flow follows the pattern: 1. Copy data from host (CPU memory) to device (GPU memory); 2. Invoke kernel to do operations on device data; 3. Copy back data from device to host. Since CUDA programming model is primarily asynchronous, during the execution of kernel, the host CPU can also execute other tasks simultaneously.

CUDA programming model provides three key abstractions: a hierarchy of thread groups, a hierarchy of memory groups and barrier synchronization. These three abstractions works together to explore maximal computing power of accelerators.

When a kernel is launched from the host side, a large number of device threads are generated and each thread executes the statements specified by the kernel [START_REF] Cheng | Professional CUDA C Programming[END_REF]. All threads invoked by this kernel composes a grid. A grid is consisted of many thread blocks. Each thread block contains a group of threads which can cooperate with each other by using block-local synchronization and block-local shared memory. To get use of shared memory, variables in kernel function need to be decorated with __shared__ attribute. Shared memory is visible among threads in a thread block. In other words, it is a tool for inter-thread communication.

The __syncthreads() function must be used to access to shared memory without producing potential data hazard.

Registers are specified with the highest bandwidth and the lowest latency since they are thread private. Variables declared in a kernel function without any other type qualifiers is generally stored in a register residing on each thread. However, register are scarce resources, there is a limitation of registers per thread depending on different GPU devices. Thus, if too much registers are used, there will be few thread blocks executing concurrently, leading to an extreme low occupancy and bad performance.

CUDA programming model provides two levels of barrier synchronization: system-level and block-level. System-Figure 2.19: CUDA Memory Hierarchy [START_REF] Cheng | Professional CUDA C Programming[END_REF].

level synchronization refers to the synchronization between host side and device side. cudaDeviceSynchronize() function can be used to block the host process until all CUDA operations have been completed. Block-level synchronization aims to barrier all threads in a thread block until they reach the same point in execution. It can be achieved by using the function __syncthreads(), as we've mentioned above.

Besides, CUDA provides asynchronous multi-streams to enable concurrency of CUDA operations such as kernel launches and memory copies. Operations within the same stream are ordered and unable to overlap. By comparison, operations marked with different streams can overlap to achieve concurrency [START_REF] Luitjens | CUDA STREAMS[END_REF].

OpenCL

OpenCL is an open industry standard and framework for programming a heterogeneous collection of CPUs, GPUs, and other processors or hardware accelerators. Unlike CUDA only supporting Nvidia GPUs, OpenCL targets to a wide range of accelerators which are considered as OpenCL devices. OpenCL device is, in fact, an abstract term describing a physical architecture. Each OpenCL devices is made up of many compute units which can be further divided into a group of processing elements.

With regard to execution model, two terms "work-group" and "work-item" are used to describe OpenCL kernel where computations occur. As an example, if we choose a Nvidia Pascal card as OpenCL device, then its compute unit is equal to CUDA streaming processor and its processing element is mapped into CUDA thread.

Directive-based Programming Models

Apart from the basic programming models introduced above, which are relatively at low-level, making them difficult to learn and to maintain, there are some directive-based programming models which are much more simple and easy to use, such as OpenMP, OpenACC, XMP [3][12], etc. In what follows we introduce several directive-based programming models that we use in our thesis.

OpenMP

OpenMP is an application programming interface that supports portable parallel programming in C/C++ and Fortran across a wide range of architectures including multicore nodes and chips, NUMA systems, GPUs, FPGAs, etc. The most useful component of OpenMP API is called compiler directives. The structure of compiler directive is shown as follow:

Fortran:!$OMP PARALLEL LOOP PRIVATE(A)
C/C++: #pragma omp parallel loop private(a)

Besides compiler directives, OpenMP also provides runtime library routines to set and query the number of threads in parallel region or even initialize the parallel loop. Furthermore, certain environment variables can be set manually to either specify loop schedule or to set number of threads, etc. Note that other programming models such as Intel Cilk+ [1] and TBB [START_REF] Kukanoc | The Foundations for Scalable Multi-core Software in Intel Threading Building Blocks[END_REF] also support data and task parallelism as OpenMP.

The execution model of OpenMP is called as the fork-join model. When an initial thread encounters a parallel region defined implicitly or explicitly by OpenMP directives, it generates a team of threads to accomplish the work inside this region concurrently. At the end of the parallel region, all members of the team join together and the program resumes being as a single thread of execution.

OpenMP 4.0 starts to provide a set of directives to perform offload computing. Such functionality is achieved by using the directive target. If a region is described by target, the data and code in this region may be offloaded to an available device. It should be noted that the host device may offload target regions to multiple target devices and OpenMP provides an interface to explicitly choose the target device. If there is no available target device, all target regions will be executed on the host device. In general, The OpenMP API defines two types target which are respectively dedicated to data mapping and code parallelization. The first one is composed of target data, target enter/exit data, target update and declare target which specify that variables or functions are mapped to or from a device. The second one consists teams, distribute (simd), distribute parallel for (simd). There are already some implementations of OpenMP API allowing for target directive. For instance, the IBM XL C/C++ V13.1.5 and XL Fortran V15.1.5 compilers are one of the first compilers that provide support for Nvidia GPU offloading using OpenMP 4.5 specification.

The OpenMP API provides a relaxed-consistency, shared-memory model where all OpenMP threads have access to a common place to store and retrieve data and each thread is allowed to have its own temporary view of the memory which allows the thread to cache variables and thereby to avoid going to memory for reference to a variable [START_REF][END_REF]. However, such relaxed-consistency model may cause data races and thus lead to incorrect result of the program. In order to avoid the memory consistency issue, OpenMP provides flush operations to enforce consistency between the temporary view and memory which can be specified using the flush directive.

OpenACC

OpenACC is a user-driven performance-portable accelerator programming model supporting implicit offload computing by uses of compiler directives, library routines, and environment variables in C/C++ and Fortran [START_REF]OpenACC-Standard[END_REF]. In contrast to CUDA, it is much easier for OpenACC to port CPU-based original code to multiple architectures without any significant structural changes, which leads to lower portability and more complexity for code maintenance.

Tetsuya Hoshino and his group studied the performance implications of OpenACC for two microbenchmarks and one real-world CFD application. The evaluations indicate that the current OpenACC compilers (PGI, HMPP, Cray) achieve approximately 50% of performance of the CUDA versions [START_REF] Hoshino | CUDA vs OpenACC: Performance Case Studies with Kernel Benchmarks and a Memory-Bound CFD Application[END_REF].

As a OpenMP-like directive-based programming model, OpenACC gives a basic directive format shown as follow:

Fortran:!$acc directive-name [clause-list] C/C++: #pragma acc directive-name [clause-list]
The compiler directives are used to accomplish data transfers between host and device, kernel execution. data or enter data/exit data are adopted to achieve structured and unstructured data mapping. parallel loop or kernel directives invoke parallel region executed on devices. OpenACC also provides a set of runtime library routines as OpenMP. These routines can accomplish a large amount of works such as getting the number of devices or retrieving the device pointer by a host pointer.

The execution model of OpenACC is similar to that of OpenMP, which is host-directed execution with attached parallel accelerators, such as GPUs [START_REF]OpenACC Application Programming Interface[END_REF]. In the case of treating the multicore CPU as a device, the initial host thread may initiate a team of threads and perform parallel execution for code regions decorated by parallel loop and kernel directives.

OpenACC exposes coarse-grained, fine-grained and vector parallelism via gang, worker, and vector abstractions. In terms of CUDA terminology, a number of gangs is equal to a set of thread blocks. Each gang (thread block) has one or more workers, referring to as warps. Vector parallelism is for SIMD or vector operations within a worker (warp), which equals to CUDA threads in a warp.

OpenACC also supports concurrent activity queues for a targeting device where the host thread may enqueue operations of data transfer or kernel execution onto multiple device queues. The host thread may continue execution as the enqueued operations are handled asynchronously with the utilization of async directive clause or wait for all operations in a queue to complete with the use of wait async directive.

The memory model of OpenACC manages implicitly data movement between the host and device based on specific OpenPACC compiler directives. If a device shares memory with the host thread, new copies of the data will be created for the device and there is no data movement occurring. If a device has a discrete memory with the host thread, the space will be allocated in device memory and the data will be copied between the host memory and device memory. Moreover, OpenACC does not provide an interface for explicit management of shared memory or hardware caches. On the contrary, they are managed implicitly by the compiler with hints from the programmer in the form of directives [START_REF]OpenACC Application Programming Interface[END_REF].

Libraries

There are many programming libraries which have been developed and used for parallel computing such as Kokkos, SYCL, RAJA, etc. We intend to introduce Kokkos and SYCL which have be implemented in our thesis work in detail.

Kokkos

Kokkos is a C++ based programming model for writing performance portable applications targeting complex node architectures with N-level memory hierarchies and multiple types of execution resources. It currently can use OpenMP, Pthreads and CUDA as backend programming models. The principal abstractions of Kokkos are parallel execution and data structure [START_REF] Edwards | Kokkos: Enabling Manycore Performance Portability through Polymorphic Memory Access Patterns[END_REF].

Concerning about parallel execution, Kokkos uses "execution space" term to describe the targeting architecture where code can actually be executed. It also provides three execution patterns which are called respectively parallel_for, parallel_reduce and parallel_scan to invoke parallelism. "Execution policy" is adopted to determine the manner how a parallel function is executed, such as setting the number of threads of kernel.

Kokkos provides a set of abstractions to manage data. Firstly, it uses "memory space" term to specify physical location of data and certain default memory access patterns. Secondly, a term named "memory layout" is used to express the mapping style of data array. Finally, Kokkos uses "memory trait" term to tell compiler how to access a data structure. There are different traits like atomic access, random access, etc.

The

Others

RAJA [START_REF] Hornung | ASC Tri-lab Co-design Level 2 Milestone Report[END_REF] is a collection of C++ software abstractions, being developed at Lawrence Livermore National Laboratory (LLNL), that enables architecture portability for HPC applications. StarPU is a task programming library for heterogeneous multi-core architectures which offers either the StarPU's high level GCC plugin pragmas, StarPU's rich C/C++ API, or OpenMP pragmas to handle task scheduling and data transfer on available CPUs and GPUs [11].

ArrayFire is a high performance software library for parallel computing with an easy-to-use API. Its array based function set makes parallel programming more accessible [START_REF] Yalamanchili | ArrayFire -A high performance software library for parallel computing with an easy-to-use API[END_REF]. OCCA (oca-rina) is an open-source library which aims to provide a unified API for interacting with backend device APIs (e.g. OpenMP, CUDA, OpenCL) targeting to different architectures (CPUs, GPUs, FPGAs) [START_REF] Medina | OCCA: A unified approach to multi-threading languages[END_REF]. DLB (Dynamic Load Balancing) library is dedicated to improve the load balance of the outer level of parallelism (e.g. MPI) by redistributing the computational resources at the inner level of parallelism (e.g. OpenMP) with LeWI (Lend When Idle) algorithm [2].

Comparison of Implemented Programming Models

Profiling Tools for Performance Analysis

Performance analysis is significant for parallel programming to mitigate bottlenecks. There are many profiling tools which are widely used for programmers to optimize their codes. In what follows several tools related to our thesis work will be introduced.

nvprof

nvprof is one of NVIDIA profiling tools that enables users to collect and view profiling data from the command-line.

As a light-weight command-line profiler, nvprof is very handy to profile CUDA applications for quick checks such as collecting a summary of run time consumption for all the kernels and memory transfer operations.

Besides, nvprof also offers users a way to perform more detailed analysis collecting a list of available events and metrics during kernel execution where an event refers to a hardware counter describing a specific countable activity, action or occurrence on a CUDA device such as active_cycles_pm (number of cycles a multiprocessor has at least one active warp) and a metric is a property of a CUDA application characterized by one or more event values such as achieved_occupancy (ratio of the average active warps per active cycle to the maximum number of warps supported on a multiprocessor).

Furthermore, nvprof is capable of remote profiling by enabling the command-line option -output-profile so as to obtain a data file which may be later imported into nvprof and the NVIDIA Visual Profiler (nvvp) for analysis.

NVIDIA Visual Profiler

The NVIDIA Visual Profiler (nvvp) offers different views including the Timeline View, the Analysis View, the Source-Disassembly View and so on to allow users to analyze potential performance bottlenecks of their CUDA applications and thus find a way to mitigate those bottlenecks and eventually optimize the performance [START_REF] Nvidia | Profiler User's Guide[END_REF].

The Timeline View consists of a set of timeline rows where each row shows the start and end times of the activities corresponding to the type of this row which may represent a process for the profiled application, a CPU thread performing either a CUDA driver or CUDA runtime API call as well as OpenACC and OpenMP directives, an estimate of the compute utilization of a device over time, a context summarizing all memcpys of a type performed during the profiling (device-to-host, host-to-device, device-to-device, and peer-to-peer), a CUDA default or created stream undergoing a series of memcpys and kernel executions.

The Analysis View provides users a straightforward analysis which may be performed in guided or unguided mode exploring some detailed performance information of a CUDA application which indicates the kernel performance limiter, the kernel latency, the PC sampling as well as the memory statistics during the execution.

The Source-Disassembly View shows the hotspots and profiling data analyzed and aggregated for a kernel at the source and assembly instruction level where the hotspots indicate the lines of source code that may generate the performance bottlenecks and they are colored according to different level of importance (low, medium or high).

In addition to the results displayed in the timeline, analysis and source-disassembly views, users may also collect the specific event and metrics values of a memory copy or kernel execution in the GPU Details View as what they can do with the command-line profiler, nvprof. Moreover, there are some other views such as CPU Details/Source View, OpenACC Details View, OpenMP Details View which are provided in nvvp to allow users to perform a complete performance analysis.

Chapter 3

Monte Carlo Neutron Transport Simulation

In this chapter, we will elaborate some aspects of nuclear reactor physics which are considered in Monte Carlo neutron transport simulation along with the context of Monte Carlo transport codes for the HPC.

Nuclear Reactor Physics

A nuclear reactor is a system which generates heat from nuclear fission where a self-propagating nuclear chain reaction happens due to the collision of a neutron with an atomic nucleus. Each collision produces energy and releases neutrons which may be used to trigger more collisions and free more neutrons to form the chain reaction.

Based on different characteristics as the types of fuel or coolant and so on, there are many types of nuclear reactor such as Pressurized Water Reactor (PWR), Boiling Water Reactor (BWR), Advanced Gas-cooled Reactor (AGR), Fast Breeder Reactor (FBR), etc. The main components of fuel are UO 2 and PuO 2 that consists uranium ( 233 U, 235 U) and plutonium ( 239 Pu) nucleus. The coolants of reactors include H 2 O, D 2 O, CO 2 , Na and so on.

Nuclear reactor physics is the physics of neutron fission chain reacting systems. It encompasses those applications of nuclear physics and radiation transport and interaction with matter that determine the behavior of nuclear reactors [START_REF] Stacey | Nuclear Reactor Physics -Second Edition[END_REF]. As a relatively mature discipline which has originated in the middle of the twentieth century and evolved for tens of years, nuclear reactor physics has now become the fundament directing the uses of nuclear energy for electricity production, neutron physics research, radiation therapy and so on. Nowadays, nuclear energy provides 10% of the world's electricity from about 440 nuclear reactors. In France, there is about 75% of its electricity derived from nuclear power plants.

The study of nuclear reactor physics is intended to understand and predict the behavior of nuclear reactor which is determined by the transport of neutrons and their interaction with nuclides composing the materials within a reactor. In other words, it can be summarized as the study of neutron transport, calculating macroscopic physical quantities such as the density of particles, the reaction rates, the heat power. Mathematically, this process can be abstracted as to solve the neutron transport equation derived from the Boltzmann transport equation that is firstly proposed by Ludwig Boltzmann in 1872 [START_REF] Stewart | An Introduction to the Theory of the Boltzmann Equation[END_REF].

Neutron Transport Simulation

Neutron Transport Equation

The neutron transport equation is a time-dependent linear equation which indicates a balance statement (neutron production equals to neutron losses) in the phase space element drdΩdE for the angular neutron density (r , Ω,E,t) at time t and kinetic energy E, with the position vector r and the unit vector in travel direction Ω (equals

to ṽ (E) |ṽ (E)|
, where ṽ (E) signifies the neutron velocity vector), as described in Equation 3.1. It should be noted that we are interested in the stationary form of the following equation without consideration of the first term @n/@t in the thesis work.

∂n/∂t z }| { 1 |ṽ (E)| @ @t (r , Ω,E,t)+ streaming z }| { Ω •r (r , Ω,E,t)+ collision z }| { Σ t (r , E ) (r , Ω,E,t) = Z 1 0 dE 0 Z 4π dΩ 0 Σ s ⇣ r , E 0 ! E, Ω0 ! Ω⌘ ⇣ r , Ω0 ,E 0 ,t ⌘ | {z } scattering + (E) 4⇡ Z 1 0 dE 0 Z 4π dΩ 0 ⌫ (E 0 )Σ f (r , E 0 ) ⇣ r , Ω0 ,E 0 ,t ⌘ | {z } fission + S ext (r , Ω,E,t) | {z } external source (3.1)
On the left hand side:

• @n/@t: The partial derivative of the number of neutrons with respect to time t inside the phase space element.

• streaming: The time variation of neutrons in the phase space element caused by travel without nuclear reactions.

• collision: The number of neutrons inside the phase space element that encounters collisions (radiative capture + scattering) at time t.

On the right hand side:

• scattering: The number of scattered neutrons at time t moved inside the space area from other areas.

• fission: The prompt fission neutrons produced inside the space area due to fission at time t where (E) refers to the energy spectrum of the fission neutrons and ⌫(E 0 ) means the average number of neutrons generated per fission.

• external source: The neutrons provided by external source that is independent of the neutron flux at time t inside the phase space element.

With reference to Σ t , Σ s and Σ f denoted in Equation 3.1, they are respectively the macroscopic total cross section, the macroscopic scattering cross section and the macroscopic fission cross section which characterize the probability of a given reaction that a neutron may experience per unit path length traveled in the material. The details of cross section will be elaborated in section 3.2.2.

Note that the scattering, radiative capture and fission mentioned above are different types of nuclear reactions that happens during the bombardment-driven process. Scattering describes the interaction without changing the nature of the target nuclide. It can be further divided into two types: elastic scattering and inelastic scattering. In an elastic scattering process, the kinetic energy and momentum of the "incident neutron-target nucleus" system is conserved and there is no energy transferred into nuclear excitation. By contrast, in an inelastic scattering reaction, the momentum of the system is preserved while the kinetic energy of the system is not conserved due to the fact that some energy of the incident neutron is transferred into nuclear excitation. Such energy transfer leaves the target nucleus in a short-lived energy state and ends up with the emission of one or more gamma rays turning the nucleus from an excited state to the ground state. In a radiative capture process, the target nucleus forms a compound nucleus in an excited state of energy by absorbing the incident neutron and is de-excited by emitting gamma rays.

As for fission, the target nucleus captures the incident neutron and splits into lighter nuclei along with the production of several neutrons and the emission of gamma rays, as well as the release of a large amount of energy.

Solutions to the neutron transport equation can be divided into fixed source problems and eigenvalue problems according to whether an external neutron source exists [START_REF] Wu | Fusion Neutronics[END_REF]. Eigenvalue problems can be further classified into (k, , , ↵)-eigenvalue problems [122][14][140] [START_REF] Carreño | The Solution of the Lambda Modes Problem Using Block Iterative Eigensolvers[END_REF] based on different ways to reach criticality of a steady-state multiplying system in the absence of external source as described in Equation 3.2.

M • (r , Ω,E)=F • (r , Ω,E) (3.2)
where:

• M : The transport operator ( Ω •r+Σ t (r , E ) R 1 0 dE 0 R 4π dΩ 0 Σ s ⇣ r , E 0 ! E, Ω0 ! Ω⌘ ).
• F : The fission operator (

χ(E) 4π R 1 0 dE 0 R 4π dΩ 0 ⌫ (E 0 )Σ f (r , E 0 )).
If we perform a k-eigenvalue calculation for Equation 3.2, the problem is reduced to determine the effective multiplication factor k ef f and its corresponding eigenfunction ef f (r , Ω,E).

M • (r , Ω,E)= F k ef f • (r , Ω,E) (3.3)
The effective multiplication factor indicates the ratio of the neutron production to the neutron loss which is used to characterize the subcritical, critical and supercritical states of a multiplying system listed as follow:

• k ef f < 1: Subcritical state, the neutron population of system decreases in time.

• k ef f =1 : Critical state, the chain reaction is self-sustaining as the neutron population remains unchanged in time.

• k ef f > 1: Supercritical state, the neutron population increases in an exponential way with time.

Nuclear Cross Section

Nuclear cross sections are the key ingredients of the neutron transport equation as they represent the probability of the neutron to interact with the crossed material.

Macroscopic Cross Section

The macroscopic cross section characterizes the probability of the neutron to interact with the material per unit path length traveled in the material. It is the sum of the microscopic cross section of all the nuclides present in the material weighted by their atomic density, as described in Equation 3.4.

Σ r (E, T)= X i N i r,i (E, T) (3.4)
where:

• r: The nuclear reaction type.

• i: The index of a nuclide that resides in the material.

• E: The energy of the incident neutron.

• T : The temperature of the material.

• N i : The atomic density of the nuclide i.

• r,i (E, T): The r-type microscopic cross section of the nuclide i at a given energy E and temperature T .

Microscopic Cross Section

The microscopic cross section of a nucleus is denoted in a simple form as r with the units of measurement barn (10 28 m 2 ) or 10 24 cm 2 where r refers to one sort of nuclear reactions such as elastic scattering ( s ), inelastic scattering ( i ), radiative capture ( γ ) and fission ( f ). The total microscopic total cross section t is used to characterize the probability that any sort of nuclear reaction will occur for a nucleus as described in Equation 3.5.

t = s + i + γ + f | {z } absorption + ... (3.5)
Figure 3.1 depicts the cross sections of several reaction types for 238 U varying with the incident energy of the neutron. As we can see, the total cross section is typically dominated by the elastic scattering cross section plus absorption cross section. Furthermore, the energy-dependent cross section can be generally divided into three energy ranges including low-energy region, resonance region and high-energy region [START_REF] Herman | ENDF-6 Formats Manual[END_REF]. The full range of cross section for a nuclide may be comprised of one or more regions. For example, 238 U contains all three regions while 1 H is only represented in the way of low-energy region.

Figure 3.1: Comparison of different types of cross sections for 238 U via JANIS 4.0 [START_REF] Soppera | JANIS: NEA JAva-based Nuclear Data Information System[END_REF] with the database set to ENDF/B-VII.0 library [START_REF] Chadwick | ENDF/B-VII.0: Next Generation Evaluated Nuclear Data Library for Nuclear Science and Technology[END_REF].

In the low-energy region (typically E<1eV), the incident energy is well below the resonance and the 1 behavior ( the velocity of incident neutron) is observed instead for both the scattering cross section and the absorption cross section when the effects of temperature and Doppler broadening is taken into account (mainly discussed in section 3.2.2.3). According to this behavior, the cross sections are tabulated as smooth functions of

1 or 1 p E .
The resonance region is one characterized by large fluctuations of several orders of magnitude. It can be divided into two regions which are respectively the resolved and unresolved resonance region. In the resolved resonance region, resonances can be seen and distinguished using methodologies such as the R-matrix formalism [START_REF] Chadwick | ENDF/B-VII.0: Next Generation Evaluated Nuclear Data Library for Nuclear Science and Technology[END_REF]. In the unresolved resonance region, there are so many resonances that they cannot be distinguished with each other and a statistical representation is adopted.

In the high-energy region, the cross sections become smooth again such as in the low-energy region. Moreover, we can find from the Figure 3.1 that at high energy the reactions involving neutron emission such as the inelastic scattering and fission become more significant comparing to ones without neutron emission such as the elastic scattering and radiative capture.

Effects of Temperature and Doppler Broadening

The effect of temperature on cross sections is caused by the thermal motion of atoms which adds the dependence of cross sections on the relative energy between the incident neutron and the target nucleus. Such effect leads to different phenomena across energy regions. In the high-energy region, the effect of temperature is negligible since the thermal motion of target nucleus is much less dominant than the motion of incident neutron. The target nucleus is considered at rest as it behaves at 0 Kelvin.

Mathematically, the Doppler broadening equation, as described in Equation 3.6, has been proposed by Cullen and Weisbin to calculate Doppler broadened cross sections along with the SIGMA1 algorithm which performs the convolution in this equation by piecewise-linear exact integration [START_REF] Cullen | Exact Doppler Broadening of Tabulated Cross Sections[END_REF].

(y, T )= 1 y 2 ✓ 1 ⇡ ◆ 1/2 Z 1 0 x 2 (x, T 0 )[e (x y) 2 e (x+y) 2 ]dx (3.6)
where

y 2 = ↵E, x 2 = ↵E r , ↵ = M k B (T T 0 )
:

• E: The incident neutron energy.

• E r : The relative energy between incident neutron and target nucleus.

• M : The target mass.

• k B : The Maxwell-Boltzmann constant equaling to 1.38064852

⇥ 10 23 m 2 • kg • s 2 • K 1 .
• T 0 : The base temperature of cross sections to be broadened.

In addition to the SIGMA1 method, there are many methods that introduces approximations to the Doppler broadening equation. For instance, the Gauss-Legendre quadrature and Gauss-Hermite quadrature approaches approximate the integral in Equation 3.6 to a form of Gaussian quadrature which contributes to higher computational efficiency comparing to the SIGMA1 method due to the avoidance of complementary error function evalua-

tions [107][76][41][65].
The Doppler broadening algorithms mentioned above are widely used in cross section processing codes such as NJOY [START_REF] Macfarlane | The NJOY Nuclear Data Processing System, Version 2016[END_REF], AMPX [START_REF] Dunn | AMPX-2000: a Cross-Section Processing System for Generating Nuclear Data for Criticality Safety Applications[END_REF], PREPRO [START_REF] Cullen | PREPRO 2019: 2019 ENDF/B Pre-processing Codes[END_REF] to prepare cross sections at the desired temperatures from 0 Kelvin cross sections. These codes also build 0 Kelvin cross sections from more elementary parameters such as resonances parameters.

Note that the simulation performed in our thesis work is based on thermohydraulic coupling which takes the effect of Doppler broadening into account for temperature feedback. As the temperature varies during the execution of simulation, a large number of temperature-dependent cross section data is required and the on-the-fly Doppler broadening method becomes essential as the conventional pretabulated method needs too much memory to store broadened cross section data.

Isotopic Depletion Calculations and Thermohydraulic Feedback

One of the main challenges for Monte-Carlo neutron transport code in reactor physics is the realization, at the fuel pin scale, of an irradiation cycle at a fuel nuclear reactor core. To perform this kind of calculation in realistic conditions, it is necessary to couple a neutron transport code with:

• An isotopic depletion solver to take into account the creation and disappearance of nuclides in the material under irradiation. This is governed by the Boltzmann-Bateman non-linear coupled system [START_REF] Bateman | The Solution of a System of Differential Equations Occurring in the Theory of Radio-active Transformations[END_REF]. This kind of calculation induces a large number of nuclides in the fuel material.

• The coupling with a thermohydraulic solver to take into account the temperature and/or density change in the fuel and the coolant. This kind of calculation needs a large number of temperature to be considered. To avoid a prohibitive memory foot print (about 1GB per temperature for the cross section), on-the-fly Doppler broadening techniques become mandatory.

Computational Methods

Two major approaches are conventionally used to solve Equation 3.1, including the deterministic method and the Monte Carlo method. The deterministic methods consist in applying numerical discretizations to the variables of the neutron transport equation. The Monte Carlo method is to transform the integral form of the neutron transport equation into the average of N values of particle histories collected through a large number of iterations [START_REF] Lux | Monte Carlo Particle Transport Methods[END_REF] which is given as Equation 3.7:

X ⇡ 1 N N X n=1 1 X k=1 x nk ! (3.7)
where X is a macroscopic quantity of interest such as a reaction rate, x nk is the observation of this quantity at the very position of k th iteration for particle n. According to the central limit theorem and the law of large numbers, X converges to the exact value of the integral when N !1.

Deterministic Method

The deterministic method introduces discretizations and approximations to simplify and obtain numerical solutions of the neutron transport equation. The discretizations are typically made with respect to the variables such as energy, angular, and spatial [START_REF] Wagner | MCNP: Multigroup/Adjoint Capabilities[END_REF][134] [START_REF] Vaz | Neutron Transport Simulation (selected topics)[END_REF].

• Energy discretization: The multigroup cross section library is used in substitution of the continuous-energy cross section data where the energy range is divided into multiple groups and the cross sections in each group become a constant.

• Angular discretization: The integral of angular is transformed to a sum of weighted neutron flux function following discretized direction vector.

• Spatial discretization: The geometry is discretized into a mesh of points and the integral of space becomes a sum of neutron flux at each point.

The introduction of discretizations makes the deterministic method a compute-efficient solution to neutron transport simulation at the expense of precision. Moreover, due to the limits of discretizations under complex conditions, the deterministic method may be unable to handle problems with specific features.

Monte Carlo Method

The Monte Carlo method is a stochastic approach to obtain the statistical properties of real-world problems [START_REF] Sawilowsky | You think you've got trivials?[END_REF].

The key feature of Monte Carlo method is the use of random sampling and statistical analysis. For example, to solve a mathematical equation, unlike conventional deterministic method which introduces numerical discretizations to achieve approximations, Monte Carlo method aggregates the results of a large amount of samples generated from a probability distribution which is constructed by (pseudo-)random numbers to simulate the reality corresponding to the equation.

Due to the fact that Monte Carlo method solves the equations abstracted from real-world problems with minimal mathematical approximations, Monte Carlo method is widely used to simulate complex systems with many coupled degrees of freedom in many domains such as statistical physics [START_REF] Binder | Applications of the Monte Carlo Methods in Statistical Physics[END_REF], computational biology [START_REF] Ojeda | Monte Carlo Simulations of Proteins in Cages: Influence of Confinement on the Stability of Intermediate States[END_REF], artificial intelligence [START_REF] Browne | A Survey of Monte Carlo Tree Search Methods[END_REF] and so on. When it is applied to solve the neutron transport problem as described in Equation 3.7, a neutron transport problem will be turned into the tracking of a large number of particle histories from their birth to their death in conjunction with the aggregated estimates obtained by tallying.

A history (random walk process) of neutron n consists of a sequence of iterations (I n0 ,I n1 ,I n2 ,...,I nKn ) in which the neutron moves a distance from its initial position to next position and then interact with the material it is flying through (Figure 3.3). K n is the total number of iterations contained in the history of neutron n and after the final iteration the neutron is desactived by leakage, absorption or energy cut. Each iteration I nk,k2Kn can be denoted by a triple ( Pnk , Dnk ,C nk ) where Pnk is the position of neutron n at the beginning of k th iteration, Dnk is the displacement of neutron n at k th iteration. C nk is one type of interactions for neutron n at k th iteration including scattering, absorption as well as none of collisions if the neutron moves across boundaries.

In the context of nuclear reactor physics, especially taking isotopic depletion into account, the most computational work during each iteration is the calculations of cross sections which are required for random samples of the distance of movement, the nuclide to be collided, the sort of interaction and the update of kinetic state of the neutron.

As for the tallying part, the observation x nk given in Equation 3.7 of a random variable X varies with the different choices of estimators. For example, in order to tally the reaction rate for reaction r within a volume without considering weights, we can use several estimators listed as follows:

• Analog estimator: observations of reaction r are accumulated. x nk equals to 1 if reaction r happens, 0 otherwise.

• Collision estimator: observations resulting in any type of collisions are scored with the ratio of macroscopic cross section for reaction r to total macroscopic cross section. x nk becomes Σ r Σ t if there is a collision, 0 otherwise.

• Track-length estimator: observations resulting in any type of collisions are scored with the multiplication of corresponding distance of movement and macroscopic cross section for reaction r. x nk equals to | Dnk |Σ r for a colliding iteration, 0 otherwise.

Furthermore, the error estimation of a tally is basically achieved by the use of variance by which we know the statistical error of the corresponding estimate. The comparison among estimates can be handled by the adoption of variance of the mean of estimates from which we can analyze the order of number of particles that makes the errors of several estimators tending to an arbitrary low level.

Overall, the very few approximations introduced make the Monte Carlo method a precise approach which is capable of handling a large number of neutron transport problems under complex conditions and performing reference calculations. However, it incurs a much higher computational cost comparing to the deterministic method widely used in the industry. For this reason, researchers and developers of many Monte Carlo transport solvers have turned to solutions by porting Monte Carlo codes to modern architectures so as to reduce time to solution and/or model more complex systems.

HPC and Monte Carlo Neutron Transport Simulation

The development of Monte Carlo neutron transport simulation in HPC mainly focus on the parallelization of Monte Carlo codes so as to benefit from high throughput of parallel computing machines and eventually yield improving performance. In what follows an explicit introduction of the state-of-the-art in the development of Monte Carlo neutron transport in HPC will be elaborated.

State-of-the-art

The general workflow of a serial Monte Carlo neutron transport simulation is illustrated in The distributed memory parallelism reached in Monte Carlo transport simulation needs to arrive at a trade-off between memory footprint on each rank and communication workloads across ranks. For example, if the geometries of Monte Carlo simulation are divided into separate domains across ranks, communications among ranks will be significantly augmented because each rank needs to address the domain-crossing issue by sending and receiving the data of out-domain and in-domain particles. By contrast, if distributed memory parallelism is only responsible for dispatching particles, there will be much less inter-communications among ranks at the expense of larger memory footprints since each rank is required to hold entire data of geometries [START_REF] Bleile | Thin-Threads: An Approach for History-Based Monte Carlo on GPUs[END_REF].

Similarly, there is always a compromise between memory footprint and achievable bandwidth for shared memory process of Monte Carlo transport simulation. For instance, the data of nuclides is usually stored as one copy and used for all multithreading processes on a shared memory system. Such strategy highly reduces memory footprints but degrades the efficiency of memory access due to the avoidance of data races and extra overheads of non-local memory access via links and so on. The corresponding performance penalty varies a lot with hardware properties such as cache hierarchy, memory access model (NUMA or UMA). It's up to developers to find a way the most adaptive to the coarse-grained parallelism of Monte Carlo transport simulation on a specific working machine taking memory footprint and bandwidth into consideration.

In addition to coarse-grained data parallelism, the parallelization at fine-grained level has also been developed

for Monte Carlo transport simulation and it has become a major aspect to dig into as modern architectures have added more features of vectorization such as SIMD and SIMT techniques. Typically, after the parallelization of coarse granularity, each processing thread is assigned to track histories of a batch of particles. To further parallelize the Monte Carlo transport codes, one may either vectorize the inner-loops in particle tracking kernel injected by cross section calculations or divide particles within a batch into vectors based on different events and execute a vector of particles in parallel once a time. The first strategy in which particles are tracked in sequence is actually called as history-based method. The second strategy where particles of a sorted event type are simulated in parallel is instead event-based method.

Parallelism

The parallelism of history-based Monte Carlo neutron transport codes used in our thesis is showed in Algorithm 1. At the distributed memory level, each MPI rank is in charge of performing an independent replication of the simulation (the first loop at line 1). The inter-communications among MPI ranks occurs only at the end of each simulation to compute the mean and variance of all the tallies. However, each MPI rank duplicates the global data such as the cross section data, the isotopic description of the materials, the description of the geometry and above all the tallies data structure, which lead to large memory footprints. For depletion calculation, tally data structure and materials description may represent several tens of gigabytes. Hence, for these kind of calculation one MPI rank is typically attached to one node.

The second loop (line 3) is executed in sequence to handle all particles in the current process batch by batch.

The third loop (line 5) is achieved by the shared memory level of parallelism as this loop contains the largest trip count (between 10 5 and 10 6 ). All threads are responsible for tracking the histories of a number of particles assigned to themselves in parallel. It is obvious that with the increase of the number of threads, the scalability will always degrade due to the common obstacles occurred at shared memory level such as memory contention, NUMA systems. Besides, since in our case we use the on-the-fly Doppler broadening method which contributes to a large amount of computational cost, the workloads of each thread may be significantly unbalanced. In order to mitigate this load balancing issue, we adopt the strategy of dynamic scheduling to assign the untracked particles to idle threads.

The while loop (line 6) simulates a particle from its birth to its death. This type of tracking approach is called history-based because we won't go to the next particle until the disappearance of the current particle.

History-based and Event-based Methods

Most of the MC transport codes are based on the history-based method, but in the 80's to make the best use of vector processors the event-based algorithm was developed which is brought up to exploit GPUs and vectorization.

Algorithm 2 shows the general procedure of history-based Monte Carlo transport simulation at the processing thread level. The outermost for loop (line 1) is executed in serial containing a while loop (line 2) tracking a particle from its birth to death. In each iteration, cross sections are firstly calculated (line 3-6) to sample distance as well as collided nuclide and interaction type (line 7). The displacement and collision take place afterward as illustrated from line 8 to line 12. The vectorizable parts reside in macroscopic cross section calculations including an outer loop for total cross section sum-up, an inner loop of energy lookup and an inner loop of integral in addition if on-the-fly Doppler broadening method is adopted. The details of these cross section algorithms are expressed in section 3.4.4.

Algorithm 2: History-based Monte Carlo transport algorithm. As an alternative solution, the event-based method was proposed by regrouping particles according to their different event types and undertaking the simulation of banked particles in parallel [25]. The general procedure of event-based method is described in Algorithm 3. The outermost while loop (line 1) ensures that all particles dispatched to a given processing thread will be tracked until they are desactived. The first inner for loop (line 2) is devoted to fill identified particles to vectors of different event types. The second inner for loop (line 5) is responsible for processing particles of an event type vector by vector. The classification of event types may be of fine or coarse granularity which is totally up to developers to define.

For example, Bleile developed an event-based Monte Carlo transport algorithm where only three event types are used: collision processing, material interface crossing, and cell boundary crossing [START_REF] Bleile | Algorithmic Improvements for Portable Event-based Monte Carlo Transport Using the Nvidia Thrust Library[END_REF][21]. Bergmann and Vujić proposed an event-based Monte Carlo transport algorithm adapting to GPUs in which more than a dozen of event types are defined [START_REF] Bergmann | Algorithmic Choices in WARP-A Framework for Continuous Energy Monte Carlo Neutron Transport in General 3D Geometries on GPUs[END_REF]. Fine-grained event types introduce less branch divergence for each event kernel while it incurs extra workload due to shuffle operations so as to regroup particles into different vectors. On the contrary, coarse-grained event types reduce particle regrouping need but make branch divergence a major point to degrade kernel performance. If we define three events including migration, scattering and absorption, the event-based Monte Carlo transport algorithm can be unfolded as shown in Algorithm 4.

Initially, all particles on the processing thread are stored in the vector of migration. In each iteration of while loop (line 2), three for loops are performed successively to process particles of different event type. The loop on the migration vector (line 3) takes responsibility of calculating cross sections, sampling distance and moving particles to collision sites. If a particle is going to undertake a collision, it will be moved to the vector of scattering or absorption and be processed in one of the following two loops. The loop of scattering (line 13) carries out scattering interactions for all particles in the vector and then shuffles them back to the vector of migration to prepare for next iteration. The loop of absorption (line 17) processes all particles of absorption, kills these particles afterward and adds newly born particles (production of fission) to the vector of migration. All three inner for loops are candidates for vectorization using SIMD or SIMT techniques.

Overall, compared with history-based method, event-based approach is better suited for vectorization on modern computers to benefit from data-level parallelism. At the same time, it introduces extra workload of consolidating surviving particles during the simulation and requires restructuring the control flow and redesigning data structure, which makes it difficult to implement with full-physics capabilities.

Recently, a significant number of studies have explored using Intel MIC and Nvidia GPU for Monte Carlo transport solvers and micro-benchmarks [START_REF] Bergmann | Algorithmic Choices in WARP-A Framework for Continuous Energy Monte Carlo Neutron Transport in General 3D Geometries on GPUs[END_REF][95][55] [START_REF] Hamilton | Continuous-energy Monte Carlo neutron transport on GPUs in the Shift code[END_REF]. Both history-based method and event-based method are im-plemented and the results are quite informative and promising. They have proved that history-based method is far more straightforward to implement than event-based method and event-based method may outperform historybased method with specific tuning strategies such as remapping data references, use of intrinsic functions and design of trivial kernel. Romano has developed a simple model to theoretically estimate the efficiency of eventbased method and observed that the good performance is yielded when the ratio of the particle bank size (e.g. migration vector size) to the vector length (e.g. AVX-512 SIMD-vector is able to store 8 double) is over 20 [START_REF] Romano | Limits on the Efficiency of Event-based Algorithms for Monte Carlo Neutron Transport[END_REF].

Algorithm 4: (Migration, scattering, absorption)-based Monte Carlo transport algorithm.

1 initially, store all particles in the vector of migration;

2 while any particles are alive on the processing thread do

Cross Section Computing Algorithms

For a given material comprised of a set of nuclides, each nuclide has its own energy grid stored in a continuous array mapping to corresponding cross section data stored in another array of the same length. The macroscopic cross sections are summed up as described in Equation 3.4. As for microscopic cross section calculations, there are mainly two approaches, which are called respectively pre-tabulated method and on-the-fly Doppler broadening method.

With reference to pre-tabulated method, an energy lookup is employed to find the nearest energy lower-upper bound in cross section tables of given temperatures and the required cross section is retrieved by performing an interpolation based on this bound. The cross sections are prepared for all the temperatures occurring in the simulated system and then loaded into memory. This kind of method induces a lot of memory operations for little computation and is therefore memory-bound.

With reference to on-the-fly Doppler broadening method, instead of storing all cross section data in memory, only the cross section data of a baseline temperature is stored and the temperature dependence is calculated on-the-fly during the simulation. This kind of method consumes less memory but is much more compute intensive.

Macroscopic Cross Section Lookup Algorithm

As shown in Algorithm 5, the macroscopic cross section of a given material accumulates with the multiplications of microscopic cross sections and atomic densities of the nuclides composing the material. This algorithm has well vectorization potential since the loop size is large enough (from tens of to hundreds of nuclides) and the workload assigned to each processing unit is equal to others which gets rid of the performance penalty caused by imbalance.

Bergmann [START_REF] Bergmann | Algorithmic Choices in WARP-A Framework for Continuous Energy Monte Carlo Neutron Transport in General 3D Geometries on GPUs[END_REF] implemented a macroscopic cross section kernel where in addition to the total macroscopic cross section, other values related to it are computed as well, including interaction distance, collided nuclide and so on.

These extra calculations are equally distributed to all processing units and thus the workload balance is nicely maintained.

Algorithm 5: Macroscopic cross section lookup.

(E)= (E k )+ E E k E k+1 E k [ (E k+1 ) (E k )] (E)= (E k+1 ) E k+1 E E k+1 E k [ (E k+1 ) (E k )] (3.8)
The key nature of such linearly-interpolable cross section data is that energy grid sizes highly depend on their corresponding nuclides. Some nuclides have hundreds of energy points (e.g. 1 With respect to workload imbalance, some methods such as unionized energy grid [START_REF] Leppänen | Two Practical Methods for Unionized Energy Grid Construction in Continuous-Energy Monte Carlo Neutron Transport Calculation[END_REF], fractional cascading [START_REF] Lund | Using Fractional Cascading to Accelerate Cross Section Lookups in Monte Carlo Neutron Transport Calculations[END_REF], hash map [START_REF] Brown | New Hash-based Energy Lookup Algorithm for Monte Carlo Codes[END_REF][125], N-ary map [START_REF] Wang | Optimization of Monte Carlo Neutron Transport Simulations with Emerging Architectures[END_REF] have been proposed to address this issue by reconstructing energy grids and supplementing additional mapping tables to effectively reduce iteration path. The basic logic of these methods is either to expand energy tables of various lengths to a same length or to compress searching ranges of long energy tables into narrower ones.

With respect to search schemes, binary search and linear search are two major approaches commonly used to perform energy lookups. The computational complexity of binary search is O(log 2 N ) while linear search is O(N ),

N is the number of elements in table. Comparing to linear search, binary search has optimized computational complexity which gains much better performance with a large problem size. However, the random memory access pattern of binary search leads to bad memory coalescing and high cache misses, giving it few chances to be vectorized with SIMD or SIMT optimizations [START_REF] Siegel | Improved Cache Performance in Monte Carlo Transport Calculations Using Energy Banding[END_REF].

Previous work done by Wang [START_REF] Wang | Competing Energy Lookup Algorithms in Monte Carlo Neutron Transport Calculations and their Optimization on CPU and Intel MIC Architectures[END_REF] has tested the search schemes, reconstruction and mapping techniques mentioned above with 512-bit SIMD optimizations. The results shows that the threshold table length for choosing between vectorized linear search and binary search is 200. Such behavior is informative and indicates that with uses of mapping techniques such as hashing and N-ary mapping, it may contribute to performance improvement if the vectorized linear search is adopted for energy lookups through narrow bins. Besides, the vectorized binary search, also called as N-ary search, yield worse performance compared with binary search by reason of higher cache misses that are caused by extra non-coalesced memory access. As for reconstruction and mapping techniques, unionized energy grid method contributes to the best performance improvement whereas N-ary mapping the worst.

In general, there are few opportunities to exploit vectorization for Monte Carlo algorithms based on pre-tabulated cross section where latency-bound energy lookups take up most of computing time. As a solution, the computeintensive FLOP work between frequent memory loads introduced by on-the-fly Doppler broadening may mitigate the latency-bound bottleneck mainly induced by the binary search in the pre-tabulated approach [START_REF] Tramm | Memory bottlenecks and memory contention in multi-core monte carlo transport codes[END_REF].

Algorithm 6: Basic pre-tabulated algorithm with binary search scheme.

Input: a given energy E, the energy grid eg[N ] and cross section grid csg[N ] of a nuclide with grid size equaling to N Result: calculated microscopic cross section (E) of the nuclide

1 k = getLowerBound(E, eg); /* 0 6 k<N 1 and eg[k] <E6 eg[k + 1] */ 2 (E)=csg[k]+ E eg[k] eg[k + 1] eg[k] (csg[k + 1] csg[k]); /* linear-linear interpolation */
Now we only consider a basic pre-tabulated method using binary search scheme without the introduction of reconstructed cross section data and mapping techniques as a comparison of on-the-fly Doppler broadening algorithm expressed in section 3.4.4.3. The general procedure of this method is illustrated in Algorithm 6 where getLowerBound() is a binary search function retrieving the last element which is lower than the input energy in the input energy grid.

On-the-fly Doppler Broadening Microscopic Cross Section Calculation Algorithm

On-the-fly Doppler broadening algorithms calculate temperature-dependent cross sections during the Monte Carlo transport simulation which reduces memory footprints at the expense of much higher computational cost. In addition to the SIGMA1 on-the-fly Doppler broadening method which performs piecewise-linear exact integration, some other methods have also been proposed to reconstruct cross sections by various multipole representations and regression models [START_REF] Hwang | A rigorous pole representation of multilevel cross sections and its practical applications[END_REF][66] [START_REF] Yesilyurt | On-the-fly Doppler broadening for Monte Carlo Codes[END_REF].

Algorithm 7: Basic SIGMA1 on-the-fly Doppler broadening algorithm.

Input: a given energy E and temperature T , the energy grid eg[N ] and cross section grid csg[N ] of a nuclide with grid size equaling to N , the elements in csg are of a baseline temperature T 0

Result: calculated microscopic cross section (E) of the nuclide A basic SIGMA1 on-the-fly Doppler broadening approach can be abstracted as a process to retrieve lower-upper bound of integral range along with another process of integral computation. In Algorithm 7, lb and ub are calculated by the same binary search function getLowerBound() described in Algorithm 6 with inputs E 1 , E 2 and eg. The compute_integral() function is an integral computation kernel approximated from Equation 3.6 which only considers calculations inside the range 4 6 x y<4 since the term e (x y) 2 becomes negligible when |x y| is large [START_REF] Cullen | Exact Doppler Broadening of Tabulated Cross Sections[END_REF].

The PATMOS Monte Carlo Prototype

There are many simulation codes developed to solve Monte Carlo transport problems, such as MCNP [START_REF] Team | MCNP-A General N-Particle Transport Code[END_REF],

TRIPOLI [27], Serpent [START_REF] Leppänen | Serpent-A Continuous-energy Monte Carlo Reactor Physics Burnup Calculation Code[END_REF], RMC [START_REF] Wang | RMC-A Monte Carlo Code for Reactor Core Analysis[END_REF], OpenMC [START_REF] Romano | The OpenMC Monte Carlo Particle Transport Code[END_REF], Shift [START_REF] Pandya | Implementation, Capabilities, and Benchmarking of Shift, A Massively Parallel Monte Carlo Radiation Transport Code[END_REF] as well as the code used in our thesis, PAT-MOS [START_REF] Brun | Patmos: A Prototype Monte Carlo Transport Code to Test High Performance Architectures[END_REF].

PATMOS is a prototype of Monte Carlo neutron transport under development at CEA dedicated to the testing of algorithms for high-performance computations on modern architectures [START_REF] Brun | Patmos: A Prototype Monte Carlo Transport Code to Test High Performance Architectures[END_REF]. One of the goals is to perform pinby-pin full core depletion calculations for large nuclear power reactors with realistic temperature fields. Numerical results have been verified by comparisons with TRIPOLI-4 R [27] and MCNP5 [START_REF] Team | MCNP-A General N-Particle Transport Code[END_REF].

The physics of PATMOS is simplified with two types of particles (mono-kinetic pseudo-particles and neutrons).

Four types of physical interactions including elastic scattering, discrete inelastic scattering, absorption and simplified fission have been implemented. The scoring part is encapsulated into a scorer class which deals with tally computation during the simulation and gathers statistical results afterwards.

PATMOS relies on a hybrid parallelism based on MPI for distributed memory parallelism and OpenMP, C++ native threads or Intel TBB for shared memory parallelism. It is entirely written in C++, with a heavy use of polymorphism in order to always allow the choice between competing algorithms such as the mix of nuclides with pre-computed Doppler-broadened cross sections and on-the-fly Doppler broadening.

Several energy lookup algorithms have been implemented in PATMOS with much vectorization effort to improve performance especially for pre-tabulated method dominated by energy lookup process. The SIGMA1 on-the-fly Doppler broadening method used in NJOY as the reference method [START_REF] Macfarlane | The NJOY Nuclear Data Processing System, Version 2016[END_REF] has also been implemented in PATMOS with some effort to reorganize the main loop of the SIGMA1 algorithm which in turn yields the reconstructed main loop a better potential for vectorization. This on-the-fly Doppler broadening algorithm has been ported to Nvidia

GPUs as a rewritten CUDA kernel which allows CPU threads offloading the compute-intensive workloads of cross section calculations to multiple GPUs and copying back results for further simulation.

In this chapter, we give a brief introduction about nuclear reactor physics, neutron transport simulation, as well as Monte Carlo method. The state-of-the-art progress of HPC and Monte Carlo neutron transport codes are discussed with respect to parallelism, particle tracking methods, cross section computing algorithms.

Chapter 4

Portable Implementation of on-the-fly

Doppler Broadening in PATMOS

Chapter 2 and Chapter 3 have presented the context of this thesis work. In this chapter, the implementation of the SIGMA1 Doppler Broadening algorithm in PATMOS with the chosen programming models are described. Then, the performance results of these different implementations are analyzed on a benchmark executed on several platforms [START_REF] Chang | Portable Monte Carlo Transport Performance Evaluation in the PATMOS Prototype[END_REF]. At last, performance profiling and optimization are presented.

Implementation of SIGMA1 on-the-fly Doppler Broadening Algorithm

This section presents the SIGMA1 Doppler broadening algorithm which allows, from a cross section at a base temperature (typically 0 Kelvin), to compute the cross section at a given temperature and a given energy.

As we have illustrated in Algorithm 7, our SIGMA1 on-the-fly Doppler broadening algorithm mainly consists in computing the indexes corresponding to the lower and upper bounds of the interval of integration (function getLowerBound()) and then compute the integral of Equation 3.6 (function compute_integral()) in which the cross section has the following form:

(x)= i + s i (x 2 x 2 i )= i+1 s i (x 2 i+1 x 2 ) (4.1)
where s i = σi+1 σi

x 2 i+1 x 2 i .
According to Equation 3.6 and Equation 4.1, we have:

(y, T )= ⇤ (y, T ) ⇤ ( y, T ) ⇤ (y, T )= 1 y 2 1 π 1/2 R 1 0 x 2 (x, T 0 ) e (x y) 2 dx ⇡ 1 y 2 1 π 1/2 P i R xi+1 xi x 2 (x, T 0 ) e (x y) 2 dx ⇡ 1 y 2 1 π 1/2 P i h R xi+1 xi x 2 i s i x 2 i e (x y) 2 dx + R xi+1 xi x 4 s i e (x y) 2 dx i (4.2)
Assume z = x y, we have:

⇤ (y, T )= P i 1 π 1/2 R xi+1 y xi y ⇣ 1 y 2 z 2 + 2 y z +1 ⌘ i s i x 2 i e z 2 dz
+ P i 1 π 1/2 R xi+1 y xi y ⇣ 1 y 2 z 4 + 4 y z 3 +6z 2 +4yz + y 2 ⌘ s i e z 2 dz = P i A (x i y) ⇥ i s i x 2 i ⇤ + B (x i y) s i P i A (x i+1 y) ⇥ i+1 s i x 2 i+1 ⇤ + B (x i+1 y) s i = A (x 0 y) ⇥ 0 s 0 x 2 0 ⇤ + B (x 0 y) s 0 + P N i=1 [s i s i 1 ] B (x i y) A (x i y) x 2 i A (x N +1 y) ⇥ N +1 s N x 2 N +1 ⇤ + B (x N +1 y) s N (4.3)
where

A (a)= 1 y 2 F 2 (a)+ 2 y F 1 (a)+F 0 (a) B (a)= 1 y 2 F 4 (a)+ 4 y F 3 (a)+6F 2 (a)+4yF 1 (a)+y 2 F 0 (a) (4.4) with F n (a)=( 1 π ) 1/2 R 1 a
z n e z 2 dz, which can be expanded as follows: The total cross section calculations account for up to 95% of total run time where the most computational cost comes from the functions erf c and exp which take up respectively 49% and 18% of total run time. Besides, comparing to the run time percentage of compute_integral function which gains around 79% of total run time, the function getLowerBound is quite negligible with only 2% run time percentage.

F 0 (a)= 1 2 erf c(a) F 1 (a)= 1 2 p ⇡ exp( a 2 ) F 2 (a)= 1 4 erf c(a)+ 1 2 p ⇡ exp( a 2 )a F 3 (a)= 1 2 p ⇡ exp( a 2 )(1 + a 2 ) F 4 (a)= 3 8 erf c(a)+ 1 2 p ⇡ exp( a 2 )( 3 2 a + a 3 ) (4.
Thus, we adopt a heterogeneous offloading strategy for our portable implementations of Monte Carlo code where only the cross section calculation is offloaded on the accelerator. The rest of the calculation remains on the host. Each accelerator can be considered as a cross section server that can be queried during the particle tracking process, by one or several host threads, to compute cross sections at a given energy and temperature. To achieve this heterogeneous offloading strategy, the main development effort is to rewrite codes which are dedicated to porting nuclide data and microscopic cross section calculations to devices via implementations of several programming models. The next section gives more details about the difference.

Implementations in the Different Programming Models

PATMOS allows two levels of parallelism via MPI + Multi-thread libraries (OpenMP, C++ threads). Each MPI rank performs history tracking for a batch of particles and the average scores of simulation are calculated via inter-node communications. In the shared memory parallellism, each CPU thread is in charge of tracking a group of particles. It allows all threads to share the non-mutable data and scores which are concurrently updated through atomic memory operations.

In order to address architectures where multiple accelerators are associated with a host upon a single node, we may use either "Multiple Threads, Single Accelerator (MTSA)" or "Multiple Threads, Multiple Accelerators (MTMA)"

strategies. The main difference between them is that MTSA requires one MPI process using a single accelerator while MTMA allows multiple threads of a MPI rank to target to multiple accelerators. MTMA avoids extra launch latency and memory footprint on a single node, which makes it a memory-friendly way for our heterogeneous offloading strategy.

Since the offloading part of our MC simulation introduces no inter-node communication, we do not consider MPI in this thesis work but only OpenMP thread + {X}, where {X} can be any languages or libraries which are capable of parallel programming on modern accelerators. In this thesis work, we have considered CUDA, OpenACC, OpenMP offload, Kokkos, SYCL or none at all (called OpenMP thread implementation in the following).

Code Architecture

The rewritten work consists in developing a set of derived classes inheriting from three base classes that include Nuclide, NeutronMedia and NeutronMediaNavigator.

The base class Nuclide engages in storing basic information related to a given nuclide such as its atomic weight (AW), energy and cross section grids of a particular base temperature as well as its table length At last, the base class NeutronMediaNavigator handles the iterations of particle history tracking which include mainly the total cross section calculations, the random samples of distance, collided nuclide, interaction type as well as kinetic state of the tracked particle after intereaction. The core part, total cross section calculation, is put in the function computeSigma() which differs a lot between its derived classes. The attribute device_id is retrieved by the thread_id in cooperation with the total number of devices nb_devices, as described in the formula device_id =t h r e a d _ i d%n b _ d e v i c e swhere % is the operator to obtain the remainder.

OpenMP Thread Implementation

The OpenMP thread version does not induce nuclide data transfer between host and devices, thus the single point to address is the OpenMP SIMD directive that we adopt to vectorize the for loop in the function compute_integral(),

as

CUDA Implementation

With respect to other portable implementations which allow for offloading microscopic total cross section calculations to accelerators, the rewritten parts are designed to accomplish the nuclide data transfer between host and devcies as well as the microscopic cross section calculations on devices. The data management of a particular nuclide is Generally, from the thread hierarchy perspective, our CUDA implementation makes use of two levels of work units warp-thread to achieve the task-level and data-level parallelism of the SIGMA1 on-the-fly Doppler broadening algorithm. From the memory hierarchy point of view, the programmable memory including registers and global memory are adopted to accomplish data management such as load/store and warp shuffle operations. Shared memory is not considered for the initial design as other programming models such as OpenACC and OpenMP offload do not yet offer programming interface to shared memory.

OpenACC Implementation

The OpenACC implementation is quite similar to the CUDA version as the CUDA runtime instructions are replaced by OpenACC directives to achieve nuclide data transfer between host and devices and offload microscopic cross section calculations to accelerators.

In order to handle the data management of memory allocation, movement and removal between host and devices, the directives #pragma acc {enter/exit} data {copyin/copyout/delete} are adopted as shown in Algorithm 11 (at lines 4, 10, 14, 18). The OpenACC runtime function acc_set_device_num() is used to set current device which may be CPU host, Nvidia GPU and so on. Moreover, the device pointers that are associated with the related host pointers by the directive #pragma acc enter data copyin can be retrieved by another OpenACC runtime function acc_deviceptr().

Concerning the parallelization of SIGMA1 on-the-fly Doppler broadening method implemented in the function computeSigma() of class OpenACCNeutronMediaNavigator, our OpenACC version adopts the directive #pragma acc parallel loop gang to parallelize the outer-loop and the directive #pragma acc loop vector the inner-loop.

The number of parallel gangs for the outer-loop is set to N by clause num_gangs(N) which corresponds to the gridsize N in the CUDA implementation. The length of vector unit for the inner-loop is defined by clause vector_length [START_REF] Chang | Performance Portability Analysis of Different Programming Models for High Performance Monte Carlo Neutron Transport[END_REF] which can be interpreted as 32 CUDA threads within a warp under the condition that the device type is set to Nvidia.

the clause reduction(operator:vars) on a loop to perform one type of operation for one or more variables. In our case, the variable sigma_tmp is involved in a sum reduction that spans the inner-loop with the vector clause and the aggregated value of this variable undergoing the sum reduction will be updated at the exit of the inner-loop.

OpenMP Offload Implementation

Like the OpenACC version, the OpenMP offload implementation also relies on the utilization of directives and runtime functions to execute some portion of the code on the device. However, there are several differences between them which need to be highlighted.

Firstly, although OpenMP target specification offers the directives to manage data transfer between host and devices including #pragma omp target {enter/exit} data, as far as we know, it doesn't give any interface to access the device addresses allocated by these directives of data management as what the function acc_deviceptr() does in OpenACC specification. Since these device pointers are required to be stored in the data structure that will be transferred between host and device during each call of the function computeSigma(), we directly use some OpenMP runtime functions dedicated to performing data management such as omp_target_alloc(), omp_target_memcpy() and omp_target_free() where omp_target_alloc() allocates memory in a device data environment and returns the corresponding device address (at lines 3, 4, 9 in Algorithm 12).

Secondly, the directives adopted in the function computeSigma() which are responsible for data transfer and kernel launch, #pragma omp target {data map/teams distribute have no clause specifying the asynchronous stream. Despite the fact that OpenMP offload does offer a couple of directives #pragma omp task, #pragma omp taskwait as well as a clause widely used to specify a set of directives denoted as nowait, such techniques only allow the corresponding code region being executed asynchronously with respect to the current host thread. Thus, OpenMP offload is not able to provide CUDA asynchronous multistreaming support as OpenACC which achieves this functionality by using the clause async(stream).

It should be noted that the map-type of the directive #pragma omp target data map that we use in our implementation is tofrom or to. The map-type tofrom signifies the mapping of data from host to device at the start of the region and the update of data from device to host at the end of the region while the map-type to refers to barely the mapping of variable to the device data environment. Furthermore, the outer-loop parallelized by the directive #pragma omp target teams distribute are specified by the clauses num_teams(N) and thread_limit(32) which define respectively N parallel teams and the upper limit of 32 threads per team.

As for the parallelization of the inner-loop in the function compute_integral(), the directive #pragma omp parallel for is used along with the clause reduction(+:sigma_tmp) so as to perform a sum reduction for the variable sigma_tmp. Besides, the clause device(idev) utilized in directives #pragma omp target [data] creates the data environment on the device of ID idev and the clause if(target [data]:device_enabled) makes the target region in the host data environment when the scalar expression represented by device_enabled evaluates to zero. The variable device_enabled is set at runtime depending on the variable nb_devices that is returned by the OpenMP API runtime function omp_get_num_devices() as device_enabled = 0, if nb_devices = 0; device_enabled = 1, otherwise.

Algorithm 12: Pseudo-code of OpenMP offload implementation. Firstly, since Kokkos implementation does not support multiple devices, our Kokkos implementation is not able to be processed on a heterogeneous architecture of a CPU associated to multiple GPUs.

In our Kokkos implementation, the data management between host and device is operated by the use of the class Kokkos::View and the Kokkos runtime function Kokkos::deep_copy(). A Kokkos view is an array of zero or more dimensions specified by a type of entries, an execution space, a memory space as well as a memory layout.

In our scenario, we utilize views of a single dimension allocated in the default execution space such as Serial, Cuda and manually set different memory spaces including Kokkos::CudaSpace and Kokkos::HostSpace with their default memory layouts such as Kokkos::LayoutLeft. The Kokkos function deep_copy() is used to manage data placement between views in different memory spaces. In order to avoid the limitations that only those views with an identical memory layout can be performed by deep_copy() correctly, the Kokkos::View::HostMirror is used as a type of view mapped to the host memory space with a compatible memory layout comparing to its original type Kokkos::View.

It should also be noted that the execution space OpenMP is not taken into account since our programming model explicitly makes use of OpenMP thread directives to accomplish task-level parallelism and Kokkos is only responsible for serial execution inside a host memory space. Moreover, like OpenACC and OpenMP offload, we use the data() method defined in Kokkos View class to obtain the device pointer of a given view conserving nuclide data which is put in other views whose entries are of a C++ struct and these views will be used in the function computeSigma().

As for the function deleteN ucOnDevices() of parent class Nuclide, it is not explicitly implemented in the child class KokkosNuclide due to the fact that all views produced in our Kokkos implementation are stored in either a C++ struct or a C++ function and Kokkos provides a reference-counting mechanism to automatically manage deallocation of Views as once the reference count of a view reaches zero, this view will be deallocated from its memory space (constructor and copy increase the reference count, destructor and assignment decrease the reference count).

Alternatively, to ensure that before the call of Kokkos::finalize(), reference counts of all views are returned to zero, we implement the function deleteN ucDataOnDevice() of the class KokkosNeutronMedia to forcefully call the destructor of the class KokkosNuclide in which views of nuclide data are conserved.

The computational kernel is launched by calling the function Kokkos::parallel_for() parameterized by a Kokkos execution policy setting its league size and team size along with a user-defined functor encapsulating the sigma1DopplerBroadening() function. In our case, the Kokkos::TeamPolicy is set to N leagues (handling cross section calculations of N nuclides in parallel) and Kokkos::AUTO threads per team (Cuda execution space equaling to a multiple of 32 and Serial equaling to the corresponding vector length of a given host processor). Moreover, Kokkos adopts a mechanism to verify if a host thread generated by Kokkos already exists in a thread pool, which may sometimes cause runtime error with the utilization of our programming model OpenMP thread + Kokkos::Serial.

For the purpose of avoiding this bug, we manually remove the verification part implemented in the Kokkos source file Kokkos_HostThreadTeam.cpp.

The compute_integral() function is implemented in the functor as a KOKKOS_INLINE_FUNCTION where the inner loop is parallelized by the Kokkos::parallel_reduce() function which accepts a parameter of an instance of Kokkos::ThreadVectorRange() defining the loop range and a parameter of KOKKOS_LAMBDA function encapsulating the calculation part as illustrated in Algorithm 13 (at lines [START_REF] Bondi | Characteristics of Scalability and their Impact on Performance[END_REF][START_REF] Brown | New Hash-based Energy Lookup Algorithm for Monte Carlo Codes[END_REF][25]. Although Kokkos provides an interface to pro-gramme "scratched pad" memory shared by threads in a team, we don't use explicitly this functionality in order to make our Kokkos version in accordance with other implementations.

It should be noted that Kokkos provides indeed the interoperability of CUDA runtime API functions. By instantiating an object of the class Kokkos::Cuda with its constructor fed with a parameter of the class cudaStream_t which is initialized by the CUDA runtime function cudaStreamCreate(), the kernels and memory operations invoked by different host threads can be overlapped as an achievement of CUDA asynchronous multistreaming in the latest version of Kokkos. Nevertheless, such technique has not been implemented in our basic Kokkos version since for all versions implemented via high-level programming languages or libraries, we intend to hide low-level APIs completely and offer users a pure high-level interface to programme.

SYCL Implementation

As shown in Algorithm 14 (at lines 1-8), the memory copies of nuclide data are invoked as memory copy kernels in shared by all work-items within the group, we declare these variables of the type cl::sycl::private_memory<T> as a substitution of the type T for the purpose of explicitly declaring them in private memory and avoiding the use of shared memory. In this way, the SYCL version can be in accordance with other implementations.

After the execution of the computational kernel, a member function update_host() of SYCL handler for memory operation is enqueued to the SYCL queue to copy the calculated data back to host.

Note that each call of the function submit() invokes a single kernel or memory operation by convention. It should be avoided to put multiple kernels and memory operations in the same scope of a function submit(). to the class cl::sycl::group so as to be invoked in the parallel_work_for_group scope [START_REF]SYCL Specification[END_REF]. Therefore, it is not suitable to our coding strategy of hierarchical parallelism and we hope such functionality can be implemented in the future release.

Additionally, SYCL allows the host to target to multiple devices. However, unlike CUDA, OpenACC and OpenMP offload which offer users the interface to explicitly associate a particular device with a host thread, the function select_device() of the class device_selector set up the device based on the highest score gained by all available SYCL devices and it will randomly select a device if more than one device obtains the same high score. As for the kernel execution order in our SYCL implementation, since all command groups are submitted to different queues where each queue is owned by a particular host thread, the order of execution is determined by the SYCL runtime.

Moreover, because all command group objects use almost the same nuclide data as the contents of requirement sets, they depend on each other and cannot be overlapped as concurrent kernels.

Implementations of Particle Tracking Methods

Heterogeneous Offloading History-based Method

Algorithm 15: Heterogeneous offloading history-based Monte Carlo neutron transport algorithm. Algorithm 15 shows the procedure of history-based method implemented in PATMOS adopting our heterogeneous offloading strategy. Microscopic cross section lookup using the SIGMA1 on-the-fly Doppler broadening algorithm is the only part offloaded to accelerators. Once the required data for cross section lookups are transferred from host to device, calculations on device begin and the host summarizes macroscopic total cross section after the results being transferred back. It is obvious that our design brings in too many back-and-forth data transfers between host and device. The size of data movement for each calculation (a group of nuclides in one material) is quite small but the large number of memcpy calls induces many launch overheads which may degrade performance in an overwhelming way.

Heterogeneous Offloading Pseudo Event-based Method

As one of general guidelines for mitigating the performance bottleneck due to host-device data transfers is to batch many small transfers into one large transfer. Thus, we developed another tracking algorithm in which one host thread treats several particles at the same time so as to request the GPU to compute the microscopic cross sections for all of these particles at the same time, in one bigger kernel. This is achieved by banking multiple particles into one group and offloading microscopic cross section lookups for all these particles. In this way, the number of data transfers can be reduced and the amount of work for each kernel is increased. To fulfill this tuning strategy, our history-based method is redesigned to a new "pseudo event-based" approach.

The details of this new method are described in Algorithm 16. The overall procedure of history tracking is reorganized into three events including a for loop of data banking, a kernel of microscopic cross section calculations and another for loop of macroscopic cross section calculations along with the distance sampling, particle movement and interaction.

Numerical Results

In this section, the numerical results of a benchmark implemented in PATMOS via different programming languages and libraries as well as two particle tracking approaches (history-based and pseudo event-based methods) executed across a set of architectures will be elaborated.

Description of the benchmark

The benchmark, called slabAllNuclides, consists in a fixed source Monte Carlo simulation of a slab discretized in 10,000 volumes. Each volume is made with the same material containing 388 nuclides, which corresponds to the number of nuclides available in the ENDFB-VII.0 nuclear data library. The number of nuclides in the material has a big impact on the macroscopic cross section calculation time. That is why we chose a material composition representative of what is encountered in depletion calculation. Moreover, 1 H and 238 U are predominant so as to obtain a neutronic behavior representative of a pressurized water reactor. The temperature of the material is 900

Kelvin. We simulate 10 batches of 500,000 particles.

In the following, all the tests are performed with the SIGMA1 on-the-fly Doppler broadening algorithm. The base temperature from which the cross section are Doppler broadened is 0 Kelvin. Performances of the different implementations (programming models and particle tracking method) are evaluated by comparing the simulation time (the run time per cycle) or the particle tracking rate (the number of particles treated per second).

Compiler Support

Since the programming languages or libraries implemented in the research (OpenMP thread, CUDA, OpenACC, OpenMP offload, Kokkos and SYCL) do not cover all architectures and each one has different degrees of maturity achieved by different compilers, it is of significance that we make a global summary of the compiler support for these programming models.

From Table 4 [START_REF] Software | ComputeCpp[END_REF]. Except for SYCL implementations mentioned above, there are another two incomplete LLVM-based implementations, triSYCL [68] and sycl-gtx [START_REF] Zuzek | sycl-gtx[END_REF] which both build on OpenCL backend to provide offload support for OpenCL devices.

Execution Environment

All intra-node tests of slabAllNuclides benchmark expressed in this chapter are carried out on five machines which are listed in Table 4.3: where HT refers to hyper-threading, SMT refers to simultaneous multi-threading.

In summary, all five machines listed above compose a set of architectures including x86 (Broadwell, Skylake, Kaby Lake G), OpenPower (Power9), Nvidia GPU (P100, V100) and Intel GPU (Intel HD Graphics 630). The first four are computing center machines while the last one is a desktop machine.

The basic configurations of these architectures are illustrated in Table 4.4. Note that the cores and threads of Nvidia GPUs listed in Table 4.4 signify respectively the streaming multiprocessors (SMs) and CUDA threads. As for the Intel Gen9.5 architecture, the number of cores refers to its number of execution units.

The compilers installed on GridCL include GCC 7.

Performance Analysis in Host Mode

Due to the fact that our implementations in PATMOS adopts a heterogeneous offloading strategy, distributing particle tracking tasks to a group of CPU threads where the microscopic cross section calculations required in all tasks may be offloaded to devices or performed on host, the performance analysis should be classified into two modes: host and offload. We start from the study of the effects of thread binding on performance in host mode.

Thread Affinity Analysis

The test is carried out on the machine GridCL where the SYCL version is not able to be executed successfully.

The benchmark is compiled by Intel C++ 17.0 with the auto-vectorization functionality enabled. The configuration of GridCL can be checked from Table 4.3 and Table 4.4.

OpenMP mainly provides three levels of thread affinity, sockets, cores and threads so as to bind a thread to a hardware socket, core or thread. In order to deliver better performance for a given number of threads, it is a priority to bind these threads to as much as hardware cores rather than hardware sockets or threads since the latter two choices may end up with many cores bound by more than one thread. Therefore, we prefer to set the thread affinity level to cores by using the command export OMP_PLACES=cores instead of sockets and threads. As for the thread affinity policy, we consider two cases, close and spread which bind the threads respectively by a close and sparse distribution across the available places.

From Table 4.5, we find that the spread policy outperforms the close policy. The performance gap becomes greater when the number of threads is less than the total number of cores on the machine (In this case, GridCL contains two sockets of 20 cores). The reason for this performance gap is that binding all threads to the cores in a single socket as what the close policy does gives each thread limited bandwidth. By contrast, binding all threads over two sockets gives each thread a higher available bandwidth and mitigates the memory-bound issue. The performance is evaluated by the simulation time (s/cycle), the lower the better.

Note that the performance of slabAllNuclides using the pseudo event-based method is similar to that using the history-based method in host mode. The spread policy is also better than the close policy. Thus, we set the thread affinity in host mode to cores/spread for all tests undergone in the following.

Scaling Analysis

Now, we carry out a series of tests for scaling analysis (strong scaling) on the machine GridCL. For host mode, the scaling performance is mainly influenced by overheads to avoid memory contention when shared resources such as nuclide data are accessed by multiple threads simultaneously. of slabAllNuclides benchmarks using the history-based method (denoted as HB). In terms of simulation time, all curves are in a form of 1 x as the run time per cycle decreases with the increase of threads. For example, the OpenMP thread version gains a factor of approximate 28x performance speedup with 40 threads comparing to that with a single thread. In terms of strong scalability, all versions resemble in each other. With the increase of threads from 10 threads to 20 threads and finally to 40 threads, their scalability drops from around 90% to 80%, and then to 70%. It can be concluded that we obtain 10% drop of strong scalability as the number of threads doubles.

Furthermore, when the hyper-threading is enabled, the scalability drops from 70% to around 45% as the number of threads varies from 40 to 80, which corresponds to the normal behavior of hyper-threading.

Figure 4.4: Scaling performance of the slabAllNuclides using the history-based method in host mode (Simulation time, the lower the better; Strong scalability, the higher the better).

Figure 4.5: Scaling performance of the slabAllNuclides using the pseudo event-based method in host mode (Simulation time, the lower the better; Strong scalability, the higher the better).

Similarly, Figure 4.5 shows the simulation time and strong scalability of slabAllNuclides benchmarks using the pseudo event-based method (denoted as PEB). Note that the bank size of pseudo event-based method is set to 100 as there is no need to perform a bank size analysis in host mode. From the perspective of performance, the simulation time of different implementations of slabAllNuclides using the pseudo event-based method is generally a little higher than the one using the history-based method. This is caused by the data banking and additional loop overheads of the pseudo event-based method. From the perspective of strong scalability, the scalability of all versions decreases in the similar way of the history-based method.

Vectorization

The performance gaps between different programming models are caused by different compilers. For instance, by using the options -Wa,-adhlng and -fopt-info with GCC compiler to generate assembly code along with the source code and show optimization information.

OpenMP directive enabled

LOOP BEGIN at sigma1xse.cxx [START_REF] Nvidia | NVIDIA TESLA V100 GPU ARCHITECTURE[END_REF]3) remark

LOOP END

We find that the for loop in the function compute_integral() is vectorized with a vector length 2 for holding two 64-bit double-precision floating point variables simultaneously as only the SSE (XMM0-XMM15) registers are used in the instructions such as vmulsd %xmm4, %xmm7, %xmm0. In summary, the OpenMP SIMD directive and GCC compiler options both fail to make the GCC compiler fully exploiting the hardware feature (Skylake processor with AVX512 extensions).

The Intel C++ compiler manages to parallelize the for loop in the function compute_integral() with a 256bit vectorization (vector length 4 for holding four variables of double type simultaneously) achieved by the use of AVX2 (YMM0-YMM15) registers, as shown in the optimization report listed above which is enabled by adding the -qopt-report=5 into compiler options. It indicates that the OpenMP SIMD directive also has no effect on Intel C++ compiler for the utilization of AVX2 extensions on GridCL as it acts for the GCC compiler. In order to accomplish a higher level of vectorization, the compiler options for auto-vectorization must be added. However, all compilers, including GCC, Intel C++, as well as PGI and LLVM Clang, fail to vectorize the inner-loop with AVX512 extensions on GridCL.

Performance Analysis in Offload Mode

Thread Affinity Analysis

Due to the GPU topology, the CPU affinity may affect the performance. By using the command nvidia-smi topo -m, we find that on GridCL, a single GPU is connected closely to a socket including 20 cores, 40 threads via PCIe.

A comparison of the performances of slabAllNuclides adopting different thread affinity policies in offload mode tested on GridCL shows that for slabAllNuclides using the pseudo event-based method PEB, both the OpenMP cores/spread and cores/close policies introduces few performance gap between each other, as illustrated in Table 4.6. The maximal ratio of the performance gap to the greater performance for a given number of threads is 4.57% achieved by 10 threads while the average percentage of performance gap equals to 1.56%.

On the contrary, the performances of slabAllNuclides HB adopting different thread affinity policies differentiate quite a lot with each other. From Table 4.6 we find that the maximal percentage of performance gap reaches up to 24.06% and the average ratio of performance gap to the higher performance is 12.66%. It indicates that comparing to pseudo event-based method, history-based method is much more influenced by the choice of thread affinity policy related to GPU topology. Note that the tests listed in Table 4.6 are performed with a single device by setting the environment variable CUDA_VISIBLE_DEVICES to a single number, 0 or 1. We can conclude that the cores/close policy works better than the cores/spread policy for slabAllNuclides HB in offload mode targeting to a single device.

As for multi-threads targeting to multi-GPUs, the effect of CPU affinity is similar to that of a single device, as described in Table 4.7. With respect to history-based method, the cores/close policy is more suitable than the cores/spread policy, improving around 16.33% performance in average. With respect to pseudo event-based method, the cores/spread policy outperforms a little the cores/close with a maximal 8.94% performance improvement and an average percentage of performance gap equaling to 2.50%.

Overall, in what follows we set the thread affinity to cores/spread for all tests carried out in offload mode to accord with the thread affinity set for tests executed in host mode, as we discussed in the section 4.4.4. where HB and PEB are respectively the abbreviations of history-based and pseudo event-based methods. The performance is evaluated by the simulation time (s/cycle), the lower the better. where HB and PEB are respectively the abbreviations of history-based and pseudo event-based methods. The performance is evaluated by the simulation time (s/cycle), the lower the better.

Bank Size Analysis

Moreover, when it comes to the pseudo event-based method, the bank size is also a key factor that influences the performance as well as the optimal thread number. In what follows we perform a bank size analysis for the CUDA We find that too small and too large bank size cannot maximize performance. For the region of bank size ranging from 1 to 100, with the increase of bank size, the performance of benchmark improves significantly, gaining a factor of approximately 20x speedup. When the bank size surpasses 100, most of performances executed by different thread numbers varies little and becomes stable, until the bank size exceeds a rather large number (in our case, 5 ⇥ 10 4 ), as the green region shown in Figure 4.6. After this green region, the performance of benchmark drops down which is due to the CPU idle threads caused by the lack of tasks to be distributed to all available threads (number of banks is smaller than number of threads). For example, when the bank size equals to 5 ⇥ 10 5 , there is only a single bank to be assigned to a CPU thread and all performances executed by multiple threads tend to the same value.

It is obvious that choosing a bank size residing in the green region is a better choice to exploit the optimal performance of our CUDA implementation in offload mode. Thus for the following tests carried out via the pseudo event-based method, we intend to fix all their bank sizes to 100.

Scaling Analysis

Concerning the offload mode, the overheads of kernel launch and memory transfer become a major point to affect scaling performance. Unlike host mode where its performance improves with the increase of threads continuously, the performance of offload mode is more like having a peak value at a given number of threads. If we use more threads, the benefits obtained from multithreading will be uncompetitive to the overheads of kernel launch and memory transfer, leading to a significant performance degradation.

In particular, for Nvidia GPUs with multi-stream support, too many host threads where each one targets to a CUDA stream make some overheads and latency of kernels and memory transfers unable to be hidden due to the total multi-stream number limited by compute resources. Therefore, it is significant to perform a scalability analysis that concentrates on finding the optimal thread number delivering the best performance. The tests are carried out with the CUDA version of slabAllNuclides on the machine GridCL in offload mode. Both the history-based method and the pseudo event-based method are taken into account. The bank size of the pseudo event-based method is set to 100 which has already been proved as a good choice to exploit the best performance.

Figure 4.7 shows that the pseudo event-based method generally exceeds a magnitude of performance achieved by the history-based method. With reference to the pseudo event-based method, when the number of threads surpasses a value (15 for a single GPU, 30 for two GPUs), its performance reaches up to a high value and becomes quite stable. With reference to the history-based method, there is a peak value obtained by a specific number of threads (4 for a single GPU, 6 for two GPUs). When the number of threads increases over the peak number, its performance degrades little by little. This scalability analysis indicates that the history-based method is a bad candidate for our heterogeneous offloading strategy while the pseudo event-based method exploits more computational capability of devices and delivers much more improving performance.

Performance Evaluation

A series of tests have been carried out on the platform set introduced in section 4.4.3 for the purpose of evaluating the performance of the different programming models. Numerical results of peak performance are illustrated in Figure 4.9: slabAllNuclides performance speedup on Cobalt-hybrid, the higher the better -Baseline for the performance is obtained in host mode with the use of pseudo event-based method and a fixed number of bank size, 100.

Figure 4.10: slabAllNuclides performance speedup on Cobalt-v100, the higher the better -Baseline for the performance is obtained in host mode with the use of pseudo event-based method and a fixed number of bank size, 100. Except for the SYCL version, other implementations of different programming languages and libraries can be compiled successfully and the performance gaps among them are similar to those obtained on GridCL. However, the performance speedups of the CUDA, OpenACC, OpenMP offload and Kokkos versions in offload mode in compar-ison with the OpenMP thread version in host mode are much higher than those tested on GridCL as well as other platforms. For example, the CUDA version may gain a factor of around 22x speedup comparing to the baseline performance. This is due to the fact that the baseline performance in host mode tested on Gorgon via the pseudo event-based method because all compilers achieves no vectorization for the inner-loop in the SIGMA1 on-the-fly Doppler broadening algorithm on OpenPower architecture while the peak performance tested on x86 architecture is obtained by Intel C++ compiler which manages to vectorize the inner-loop with AVX2 extensions.

Besides, the best performance of the OpenMP offload version obtained on Gorgon is better than its peak performance tested on GridCL while the performances of the CUDA and OpenACC versions gained on Gorgon are similar to the peak values tested on GridCL. For instance, the peak performance of the OpenMP offload version tested on GridCL in offload mode with a single device gains 2.6 ⇥ 10 2 particles/s while its peak performance obtained on Gorgon reaches up to 6.5 ⇥ 10 2 particles/s. This may be caused by the OpenMP offload functionalities achieved by different compilers including LLVM Clang and XLM C++.

Figure 4.12: slabAllNuclides performance speedup on Intel NUC, the higher the better -Baseline for the performance is obtained in host mode with the use of pseudo event-based method and a fixed number of bank size, 100.

At last, Figure 4.12 shows the performance speedups of different programming models executed on Intel NUC equipped with an integrated graphical card (Intel Gen9.5). The SYCL version is the single implementation capable of offloading computational kernel to Intel Graphics. However, the corresponding performance in host and offload modes is much worse than the baseline performance obtained by the OpenMP thread version executed in host mode, with a factor of 0.1-0.2x speedup. Concerning other programming models such as OpenMP thread, Ope-nACC and Kokkos which are successfully processed in host mode, the OpenMP thread and Kokkos versions are

Profiling of HB and PEB

The number of threads is set to 10 and the bank size for the pseudo event-based method is fixed to 100. The total number of particles is reset to 5 ⇥ 10 3 . We use nvprof to perform an analysis of the CUDA version of slabAllNuclides using both the history-based method and the pseudo event-based method in default summary mode. The performance is evaluated by the execution time (s), the lower the better.

From Table 4.9 we find that the execution of the CUDA kernel function sigma1DopplerBroadening() and memcpys (HtoD and DtoH) via PEB saves around 90% run time comparing to those via HB. This is in accordance with the performance results introduced in the previous section. To get more detailed information about the CUDA kernels of both HB and PEB methods, we use nvvp to perform a complete profiling analysis which is illustrated in Table 4.10.

Several metrics used in the table above are described as:

• Global memory load/store throughput.

• Global memory load efficiency: ratio of requested global memory load throughput to required global memory load throughput expressed as percentage.

• Global memory store efficiency: ratio of requested global memory store throughput to required global memory store throughput expressed as percentage.

• Achieved occupancy: ratio of the average active warps per active cycle to the maximum number of warps supported on a multiprocessor.

• Branch efficiency: ratio of branch instruction to sum of branch and divergent branch instruction.

• Warp execution efficiency: ratio of the average active threads per warp to the maximum number of threads per warp supported on a multiprocessor.

• Double-precision floating-point efficiency: ratio of achieved to peak double-precision floating-point operations.

• Execution dependency: An input required by the instruction is not yet available. • Memory dependency: A load/store cannot be made because the required resources are not available or are fully utilized, or too many requests of a given type are outstanding.

With reference to the analysis of the history-based method, the results indicate that the performance of the kernel sigma1DopplerBroadening() is most likely limited by the latency of arithmetic or memory operations which delivers low compute throughput and memory bandwidth comparing to the peak performance of V100, as described in Table 4.10 (compute throughput and global memory load/store throughput). The branch efficiency and warp execution efficiency show that the issues of divergent branches and idle threads per warp for the kernel function sigma1DopplerBroadening() is not the key bottleneck limiting the performance. Besides, the grid size of the kernel is too small to keep the GPU busy so as to hide compute and memory latency. As a result, the achieved occupancy of the kernel via the HB method gains 6.22% while the theoretical occupancy equals to 43.8% as the number of blocks per SM is limited by the usage of registers per thread (69, where 64 is the maximal registers per thread that can be used to allow 32 blocks executing concurrently on each SM). The execution dependency and memory dependency are the most important stall reasons that limit the performance where the values corresponding to these two types of stalls are 48.57% and 33.11%, indicating that the latency of the kernel is mainly introduced by execution and memory dependency stalls and it may be potentially reduced by increasing instruction-level parallelism and optimizing memory alignment as well as access patterns. The memory dependency stall can also be reflected by Chapter 4 Portable Implementation of on-the-fly Doppler Broadening in PATMOS 93 the metrics of global memory load and store efficiency, which are respectively 50.39% and 25.00%, meaning that for each memory transaction call, only a half or a quarter of bus volume (32 bytes in total) can be actually used due to memory misalignment and inadequate access patterns.

With reference to the analysis of the pseudo event-based method, although the compute throughput and memory bandwidth improve quite a lot comparing to the metric values of the HB method (64-bit floating-point efficiency varies from 4.42% to 21.94%, global memory load throughput increases from 33.03GB/s to 262.20GB/s, global memory store throughput augments from 0.36GB/s to 2.87GB/s), its performance is still limited by the latency of arithmetic or memory operations. The branch efficiency and warp execution efficiency remain the same as the HB method. The execution dependency and memory dependency are still the most two important reasons causing the performance degradation, taking up respectively 37.79% and 17.97% stalls during the kernel execution, which contribute to less stalls relative to the HB method. Moreover, since the grid size of the kernel for the PEB method is much larger than that of the HB method, the achieved occupancy via the PEB method gains 31.86% which surpasses the value of the HB method significantly and leads to a considerable performance improvement.

In summary, the PEB method has significantly improved the compute throughput and memory bandwidth in comparison with the HB method. The achieved occupancy of the PEB method gains a factor of 5 enhancement comparing to the HB method so as to hide more compute and memory latency. However, even for the PEB method, the utilization of compute resources of Nvidia V100 is still at the medium level and the kernel performance is bound by instruction and memory latency. From one side, the occupancy of the kernel may be improved by reducing the number of registers used for each CUDA thread or by restricting the maximum number of registers per thread with the flag -maxrregcount. From the other side, the execution dependency and memory dependency stalls may be reduced by adding instruction-level parallelism and optimizing memory alignment and access patterns. In what follows we continue to perform a series of performance analysis of different implementations of slabAllNuclides (OpenACC, OpenMP offload, Kokkos) using the PEB method to know better about the performance bottlenecks of all these implementations.

Profiling of different programming models

The following performance analysis has been carried out with the OpenACC, OpenMP offload, Kokkos versions of slabAllNuclides via the PEB method. The input parameters are set to 5 ⇥ 10 3 particles, 10 threads, and the bank size is fixed to 100.

OpenACC Analysis

From Table 4.11, we can see that the computational kernel of the OpenACC version uses 110 registers for each thread which leads to a smaller theoretical block limit and a lower achieved occupancy comparing to the CUDA The warp execution efficiency of the OpenACC version gains around 40% which is much less than the CUDA version (around 100%), indicating that the OpenACC kernel generates more intra-warp divergence than the CUDA kernel.

OpenMP offload Analysis

The OpenMP offload kernel is similar to the OpenACC kernel. From the aspect of the usage of registers, it uses 200 registers per thread which makes the theoretical block limit decreasing to 10 and results in a lower achieved occupancy in comparison with the OpenACC version. The memory bandwidth of the OpenMP offload kernel degrades significantly due to the lack of CUDA multi-stream support. It should also be noted that the OpenMP offload version makes no use of shared memory as the CUDA version. The warp execution efficiency of the OpenMP offload kernel is 27% which is even less than the OpenACC version, indicating that the OpenMP offload version also has much higher intra-warp divergence than the CUDA version.

Kokkos Analysis

The Kokkos version use 70 registers for each thread which is in accordance with the CUDA version. As a results, its achieved occupancy is almost equal to that of the CUDA version. The memory bandwidth of the Kokkos kernel is similar to that of the OpenMP offload kernel as they both provide no support of CUDA multi-streaming. However, the global memory load efficiency of the Kokkos is 25% which is less than other versions, leading to a larger number of memory transactions for memory load operation. This may be the reason that the performance of the Kokkos version is not competitive with the OpenMP offload version. Besides, the Kokkos kernel also makes no use of shared memory as the CUDA kernel. As for the warp execution efficiency of the Kokkos version, it is equal to 100%, meaning that there is no intra-warp divergence occurring during the execution of the Kokkos kernel. The performance is evaluated by the simulation time (s/cycle), the lower the better.

Note that Kokkos also provides the interoperability of CUDA runtime API functions to enable the support of CUDA asynchronous multistreaming. We intend to implement this technique in our Kokkos version as an optimization and the timeline view of nvvp shows that the optimized Kokkos version indeed makes use of CUDA multi-streams and the performance of the optimized Kokkos version has improved 10% relative to the performance of the basic Kokkos version using the PEB method, as showed in Table 4.12.

Optimizations of the CUDA version

According to the performance analysis carried out in the previous section, we intend to accomplish several optimizations for the purpose of mitigating the bottlenecks of the kernel of slabAllNuclides. In order to make fully use of CUDA's feature set, the optimization effort has been done on the CUDA version of slabAllNuclides. The input parameters are set to 5 ⇥ 10 5 particles, 20 threads. The bank size for the PEB method is still set to 100.

Zero Copy

Since our design of algorithms, especially for the history-based method, requires many small data transfers during the simulation, the zero copy functionality can be used to obtain higher bandwidth between the host and the device by using pinned memory. Instead of allocating data in pageable host memory which need to be firstly copied to a temporary pinned memory and then transferred to device memory from the pinned memory once a data transfer between host and device is invoked, the zero copy functionality manages to directly allocate the data in pinned memory so as to avoid the cost of data transfer between pageable and pinned memory with the utilization of CUDA runtime API cudaM allocHost() or cudaHostAlloc() and cudaF reeHost().

Register Usage

Both the HB and PEB methods own the problem of the overuse of registers per thread which may prevent the kernel from fully utilizing the device. Since in our case, the maximum number of registers that a thread can use to execute 32 blocks concurrently on a SM is 64, we may use the -maxrregcount=64 flag to manually set the maximum number of registers for each thread.

Unified Memory

The CUDA Unified Memory is a single memory address space accessible from both host and device which allows the data to be read and written from CPUs or GPUs without explicit invocation of memory transfer operations. In order to retrieve a pointer accessible from any processor, the CUDA runtime API function cudaM allocM anaged() need to be called. In our case, the input data required for kernel execution may be stored in managed memory.

Multiple Kernels

Because the computational kernel of the function sigma1DopplerBroadening() consists mainly two parts, a procedure of precalculation to retrieve lower-upper bound and increment several cross section values to the final result and a procedure of integral computation which handles the execution of several inner-loops. When a warp deals with the cross section calculation of a nuclide, all CUDA threads within this warp can only execute the same workloads of the precalculation part simultaneously which makes no contribution to parallelize this procedure so as to gain improved performance.

As a solution, we may separate the original kernel into two small kernels which are dedicated to handling the precalculation part and integral computation part respectively. Each CUDA thread in the precalculation kernel is responsible for precalculating the data of a set of nuclides required in the integral computation kernel. Each warp in the integral computation kernel is responsible for summarizing the final cross section values of a set of nuclides from inner-loops.

Fine-grained Kernel

As the performance bottlenecks of the kernel using the PEB method may be alleviated by adding instruction-level parallelism and optimizing memory access patterns, we may redesign our computational kernel into a more finegrained one where a set of CUDA threads within a warp are assigned to calculate the cross section of a nuclide. For example, we can use 4 or 8 threads to calculate the cross section value of a nuclide and thus a warp can handle 8

or 4 tasks at the same time. In this case, it is preferable to perform the sum reduction for the inner-loops with the utilization of shared memory. The newly designed function of sum reduction is shown above where NUCS_PER_WARP threads within a warp into two groups where the two groups are responsible for retrieving respectively the lower and upper bounds of all tasks assigned to this warp.

The value of NUCS_PER_WARP can be set to {1, 2, 4, 8, 16}. If we extend this value to 32, meaning that each thread calculates the cross section of a nuclide on its own and there is no vectorization for the inner-loops, we will have no more need to use shared memory since the intra-warp communication among threads is reduced to zero.

A series of tests of slabAllNuclides concerning all these cases have been executed and the corresponding results are listed in Table 4.13. It is obvious that the performance of both the HB and PEB methods improves firstly with the increase of the number of tasks for each warp. After it reaches up to the peak value obtained with NUCS_PER_WARP equaling to 4, the performance degrades little by little. The performance is evaluated by the simulation time (s/cycle), the lower the better.

We can say that there is always a tradeoff between the increase of instruction-level parallelism and the increase of branch divergence. The performance of the CUDA version using the PEB method under the condition that NUCS_PER_WARP = 4 gains approximately 21% improvement comparing to the one given by NUCS_PER_WARP = 1.

Comparison of performances of the CUDA version via different optimizations

According to all optimizations discussed above, we denote them in a simple way to prepare for a thorough performance evaluation.

• warp: the original implementation of the CUDA version.

• zero-copy : an optimized version with the zero copy functionality.

• maxrregcount-64: an optimized version with the maximum usage of registers per thread set to 64.

• maxrregcount-32: an optimized version with the maximum usage of registers per thread set to 32.

• unified-memory : an optimized version with the use of unified memory.

• multi-kernels: an optimized version dividing the original kernel into two small kernels.

• 4-nucs-per-warp: an optimized version of a fine-grained kernel where each warp handles 4 computational tasks simultaneously.

The corresponding results are illustrated in Table 4.14. We can find that the zero-copy strategy improves significantly the performance via the HB method with a factor of 2x speedup while it degrades around 50% of the performance via the PEB method. The performance is evaluated by the simulation time (s/cycle), the lower the better.

As for the strategy of setting up maximum registers per thread, the maxrregcount-64 enhances a little performance for both the HB and PEB methods comparing to the original warp version while the maxrregcount-32

improves the performance adopting the HB method and weakens around 40% of the performance relative to the warp version using the PEB method. It is recommended to make fully use of available registers for each thread without limiting the theoretical maximum occupancy a kernel may achieve since the local memory of a much lower memory bandwidth comparing to registers will be used to accomplish memory operations for spilled data.

With reference to the unified-memory optimization, the results shows that our programming strategy is not a good candidate for the utilization of CUDA Unified Memory since the performance of both the HB and PEB methods using unified memory obtains approximately 60% performance degradation in comparison with the original warp version.

In addition to the optimizations discussed above, the implementations optimized by multi-kernels and 4-nucsper-warp manage to outperforms the original one with a factor of 1.2-1.3x speedup for the PEB method. However, these two optimizations make no contribution to achieve the performance improvement for the HB method, as depicted in Figure 4.13.

Chapter 4

Portable Implementation of on-the-fly Doppler Broadening in PATMOS 100

Overall, we can draw several conclusions listed as follows:

• With respect to the history-based method, the zero-copy strategy is the most useful way to optimize the kernel and improve the performance.

• With respect to the pseudo event-based method, the multi-kernels and 4-nucs-per-warp strategies become the most favorable ways to mitigate performance bottlenecks and achieve optimizations.

• The maxrregcount optimization should be used carefully so as to fully exploit the compute resources of the device and avoid too much usage of local memory.

• The CUDA Unified Memory is inappropriate to be used in our case with the heterogeneous offloading strategy. In conclusion, this chapter introduces explicitly the portable implementations of Monte Carlo neutron transport Chapter 4

Portable Implementation of on-the-fly Doppler Broadening in PATMOS 101 code with the utilization of a heterogeneous offloading strategy via different programming models. The basic performance evaluation and profiling have been done for the purpose of better understanding the bottlenecks of our portable implementations and achieving several optimizations.

Chapter 5

Metrics Comparison among Programming

Models

Based on the portable implementations of a set of programming languages or libraries on the Monte Carlo neutron transport code PATMOS, we need to find ways to evaluate them so as to choose the most suitable one to our Monte Carlo code and accessible machines.

Till now, there are few researches addressing the evaluation of Monte Carlo neutron transport codes on supercomputers with reference to 3P (performance, portability, productivity), not to mention the evaluation of Monte Carlo neutron transport codes adopting a partial offloading strategy. To our knowledge, the only work related to the implementation of MC transport solver with concern of portability is carried on by Bleile and his group [START_REF] Bleile | Investigation of Portable Event-based Monte Carlo Transport Using the Nvidia Thrust Library[END_REF] where they developed a monoenergetic event-based MC neutron transport solver relying on the Nvidia Thrust library and discovered that the Thrust version can only obtain a maximum of 36% performance comparing to the CUDA version.

In order to obtain the correct feedback and take it as reference for future development of next generation Monte Carlo neutron transport code, it is significant to perform an evaluation of our Monte Carlo neutron transport codes in terms of 3P. Thus, in this chapter, we intend to give an explicit evaluation of portable implementations of programming models on PATMOS in respect of performance, portability, and productivity which are leveraged by metrics proposed by ourselves or other researchers. A part of work has already be presented in the conference SNA+MC and will be published recently [START_REF] Chang | Performance Portability Analysis of Different Programming Models for High Performance Monte Carlo Neutron Transport[END_REF].

Definitions of Metrics

According to the 3P Principle, a programming language or library cannot satisfy performance, portability and productivity at the same time and there are always trade-offs among them. For example, Figure 5.1 intuitively showed the relation among 3P for the programming models implemented in our research. The right arrow indicates that after several explicit optimizations, the implementation of a programming model on an application may gain better performance but become less productive. From the perspective of quantitative evaluation, some researchers have already done impressive work in terms of performance portability and productivity. For example, Pennycook and his group have proposed a definition for performance portability which highlights the objectivity and measurability as well as a metric to quantify the performance portability of an application across a set of computing architectures [97][98]. This metric has been used in many studies of performance portability [100][40][73] [START_REF] Zhao | Delivering Performance-Portable Stencil Computations on CPUs and GPUs Using Bricks[END_REF]. The definition is given as: Definition 5.1.1 Performance portability is a measurement of an application's performance efficiency for a given problem that can be executed correctly on all platforms in a given set.

Its corresponding metric measures the performance efficiency of an application for a case set (a, p, H) where performance efficiency includes application efficiency, architectural efficiency and so on. The metric is defined as follows:

P(a, p, H)= 8 > < > : |H| P i2H 1 ei(a,p) if i is supported 0 otherwise (5.1)
where e i (a, p) signifies the performance efficiency of the application a solving a problem p on a platform i in the platform set H.

The reason of e i (a, p) referring to the performance efficiency is that performance efficiency is objective enough to describe the ability of an application exploiting the computing power of a platform for the purpose of solving a problem. However, the use of P always involves comparisons of estimated values between different applications or problems such as P(a 1 ,p,H) against P(a 2 ,p,H). It is meaningless that we only compute an aggregated value by P without doing any comparison. In our opinion, the metrics of absolute performance like execution time or throughput may also be used by P.

Furthermore, P performs the harmonic mean to ensure that the estimated result is directly proportional to the sum of scores across H. Nonetheless, this argument remains to be discussed. Although Smith has proved in his paper [START_REF] Smith | Characterizing Computer Performance with a Single Number[END_REF] that harmonic mean should be used for summarizing performance expressed as a rate (details will be expressed in section 5.1.2.1), which is corresponding to the performance efficiency, such correct aggregation is used to capture an intrinsic property of a particular platform. When it comes to the aggregation of values across a platform set, it is not precise to say that harmonic mean summarizes correctly the measured values to reflect an intrinsic property of a platform set. In fact, whether such intrinsic property of a platform set exists is itself a question to ask. Thus, P may also adopt other aggregating methods such as arithmetic and geometric means to estimate the data set.

According to the two arguments addressed above, we have come up with an original idea to establish a generic model which defines a variety of metrics with regard to 3P, allowing the utilization of different aggregating methods performing from multiple dimensions such as an application set, a problem set, a platform set.

Generic Model

Before the establishment of the generic model, we first give several crucial notions which are repeatedly used for the following definitions. where C(a, p, h) is a metric parameterized by an application, a problem and a platform.

Given a characteristic c, an application set A = {a}, a problem set P = {p} and a platform set H = {h},w e have three functions parameterized by three fixed variables along with a variable set:

• A(c, A, p, h): A generic function that evaluates the characteristic c of a set of applications A solving the problem p running on the platform h.

• P(c, a, P, h): A generic function that evaluates the characteristic c of the application a solving a set of problems P running on the platform h.

• H(c, a, p, H): A generic function that evaluates the characteristic c of the application a solving the problem p across a set of platforms H.

A(c, A, p, h) is generally used to describe an intrinsic property of a given problem and platform. For example, we can set OpenMP thread version of slabAllNuclides via the history-based and pseudo event-based methods as two elements in A and fix a particular problem p. Then we compare the results of A(c, A, p, h 1 ) and A(c, A, p, h 2 ) where h 1 and h 2 differentiate between each other on the type of compiler. In this way, the general influence of compilers for history-based and pseudo-event based applications on an expected property can be evaluated numerically.

P(c, a, P, h) indicates a general property of a given application and platform solving a set of problems. For instance, if we set P = {p i } i21,2,...,10 as ten problems with different number of threads, P(c, a, P, h) will become a value reflecting the scalability of the application a running on the platform h. The comparison between P(c, a 1 ,P,h) and P(c, a 2 ,P,h) shows the corresponding scalabilities of the applications a 1 and a 2 on a particular platform.

Similarly, H(c, a, p, H) is calculated to highlight a property of a specific application and problem across a platform set. Let H be a group of architectures, then H(c, a, p, H) leverages a property relevant to portability for a particular application and problem. Now we transfer A, P, H to three operators F A , F P , F H as described in Equation 5.3 which signify the computations of a data set C(a, p, h) retrieved from one of three dimensions including application, problem and platform so as to eventually calculate an aggregated metric. 5.5:

A(c, A, p, h)=F

F am • C(a, p, h)= 8 > > < > > : 1 N N X i=1 C i (a, p, h)( a) P N i=1 ! i C i (a, p, h)( b) (5.5)
where Equation 5.5(a) is the standard version of arithmetic mean metric using equal weightings while Equation 5.5(b) is the weighted version. Each C i (a, p, h) is multiplied by a weight ! i and the sum of {! i } i2N equals to 1. Note that the standard version is a special case of weighted version with ! i = 1 N .

Geometric Mean

Definition 5.1. [START_REF]POWER9 -Microarchitectures -IBM[END_REF] The metric F gm is an operator computing the geometric mean of a data set of N values {C i (a, p, h)} i=1,2,...,N , defined by Equation 5.6:

F gm • C(a, p, h)= 8 > > < > > : N Y i=1 C i (a, p, h) 1 N (a) Q N i=1 C i (a, p, h) ωi (b) (5.6)
where Equation 5.6(a) is the standard version of geometric mean metric using equal weightings while Equation 5.6(b) is the weighted version. Each C i (a, p, h) is raised to the power of a weight ! i and the sum of {! i } i2N equals to 1. The standard version is a special case of weighted version with ! i = 1 N .

Harmonic Mean

Definition 5.1. [START_REF] Sycl Overview | [END_REF] The metric F hm is an operator computing the harmonic mean of a data set of N values {C i (a, p, h)} i=1,2,...,N , defined by Equation 5.7:

F hm • C(a, p, h)= 8 > > > > > < > > > > > : N P N i=1 1 C i (a, p, h) (a) 1 P N i=1 ! i C i (a, p, h) (b) (5.7)
where Equation 5.7(a)(b) are respectively the standard and weighted versions of harmonic mean metric using equal weightings. Each 1 Ci(a,p,h) is multiplied by a weight ! i and the sum of {! i } i2N equals to 1. The standard version is also a special case of weighted version with

! i = 1 N .
Generally, the arithmetic mean metric F am calculates the sum of a collection of numbers divided by the count of numbers in the collection [START_REF] Jacobs | Mathematics: A Human Endeavor[END_REF]. The geometric mean metric F gm is often used to aggregate the values of a data set containing different numeric ranges [START_REF]TPC-D -Frequently Asked Questions[END_REF] so as to indicate the central tendency by using the products of these values.

The harmonic mean metric F hm is appropriate for situations when the average of rates is desired by calculating the inverse of the sum of the inverses of a collection of values. In order to apply these metrics of Pythagorean means correctly for the purpose of accurate capture of a property, it is significant to exploit conditions that these metrics are applicable to.

Adaptability Analysis

We say a metric is adaptive to a data set C(a, p, h) if it is capable of evaluating a set of elements {C i (a, p, h)} i=1,2,...,N with correctness. Otherwise, the observed result may be overestimated or underestimated. Based on several studies [112][70] addressing this issue, we arrive at the following definitions.

Definition 5.1.9 A single-platform property is characterized by F summarizing a data set C(a, p, h) which varies on the dimension of application or problem. By contrast, a property captured by F with its data set C(a, p, h) varying on the dimension of platform is called across-platform property.

A single-platform property indicates an intrinsic characteristic of the platform whereas an across-platform property expresses an extrinsic property among platforms. Thus, the adaptabilities of arithmetic, geometric and harmonic mean metrics dependent on the property to be characterized. Before further discussions, we firstly define the form of property used in our model.

Definition 5.1.10 A property characterized by C(a, p, h), no matter it is single-platform or across-platform, is generally represented by a base-unit quantity or its inverse as well as a ratio of one base-unit quantity to another base-unit quantity.

F am is not suitable to handle this case accurately. But if we apply F hm to the data set, the result resembles in Equation 5.9 indicating that F hm is adaptive to C(a, p, h) retrieved by second strategy of measurement.

F hm • C(a, p, h)= N P N i=1 ti P = PN P N i=1 t i / 1 P N i=1 t i
(5.12)

Note that the conclusions drawn above expands the ones drawn in Smith's previous work [START_REF] Smith | Characterizing Computer Performance with a Single Number[END_REF] where he concludes that arithmetic mean is a bad candidate for summarizing floating point operations per second and this quantity should be estimated by harmonic mean instead. The reason he arrived at such conclusions is that in his paper, each case is assumed to perform the same number of floating point operations during the execution time. This assumption is actually our second strategy of measurement proposed above.

In fact, the standard versions of F am and F hm adopting two strategies of measurement can be considered 5.10 and Equation 5.12, which are expressed as

as
F • C(a, p, h)= 8 > > > > < > > > > : N X i=1 ! i p i t i = N X i=1 t i P N i=1 t i ! p i t i = 1 N P N i=1 Npi P N i=1 ti = 1 NT P N i=1 p 0 i 1 P N i=1 ! i ti pi = 1 P N i=1 ⇣ p i P N i=1 p i ⌘ t i p i = N P N i=1 pi P N i=1 Nti = PN P N i=1 t 0 i 9 > > > > = > > > > ; = P N i=1 p i P N i=1 t i (5.13)
where:

• ! i = t i P N i=1 t i or p i P N i=1 p i • p 0 i = Np i • t 0 i = Nt i • T = P N i=1 t i • P = P N i=1 p i
According to Equation 5. [START_REF]OpenACC Application Programming Interface[END_REF]

Adaptability to Across-platform Properties

When it comes to address an across-platform property, it is not easy to decide which metric is the most suitable one for a given data set. From the one side, since each measured value indicates the expected property of the corresponding platform, we don't need paying attention to correctly summarize these values like what we have to convert the product of values to the sum of logarithms of values, one may need to pay more attention to the use of zero value.

F gm also has a property that F am and F hm do not have, that is, the geometric mean of the ratio of two data sets (the geometric mean of a normalized data set) equals to the ratio of the geometric mean of one data set to the geometric mean of the other data set.

F gm • ✓ C 0 (a, p, h) C(a, p, h) ◆ = F gm • C 0 (a, p, h) F gm • C(a, p, h) (5.18)
This feature allows F gm being capable of giving consistent results for a normalized data set. However, if we want to aggregate performance values such as time, floating point operations, it is recommended to avoid the use of the geometric mean since multiplications of times or floating point operations are meaningless [START_REF] Smith | Characterizing Computer Performance with a Single Number[END_REF].

Metrics of Portability, Performance, Productivity

Metrics of Portability

After the establishment of generic model F , now we take a look at the metrics of portability. We give a classical definition of portability in HPC domain. According to Equation 5.19, we have:

F gm • C(a, p, h)=F hm • C(a, p, h)= 8 > > > < > > > : 0 if 9h 2 H, h doesn't support (a, p) 1 if 8h 2 H, h supports (a, p) F am • C(a, p, h)= |S| |H| (5.20)
where S is the sub-set of H containing all supported platforms in H for (a, p). |S| is the total number of supported platforms and |H| is the total number of platforms.

Chapter 5

Metrics Comparison among Programming Models 112

Note that the harmonic mean metric F hm only applies to positive number. In this scenario, we use limits 1 0 !1 and 1

1 ! 0 to allow calculations of zero values for F hm . We can say that both F gm and F hm accurately describe the portability of a couple (a, p) on a set of platforms H. On the contrary, F am fails to address the feature of portability across multiple platforms.

Metrics of Performance

Generally, performance is a key property of interest to HPC community. Here we a generic definition of performance for a triple (a, p, h). Definition 5.1. [START_REF]OpenACC Application Programming Interface[END_REF] Performance refers to any measurable property reflecting the running efficiency of an application a to solve a problem p correctly on a platform h.

There are a variety of metrics capable of characterizing the performance, commonly including time-based metrics (e.g. execution time) and energy-based metrics (e.g. energy consumption). We list several time-based metrics as follows:

• Time/Duration: A basic metric describing the performance. It is the easiest metric to be measured, with sec as its unit of measurement.

• Throughput: A metric denoting the total amount of floating point operations which can be processed per unit of time. It can be measured by some profiling tools such as Intel Advisor and nvprof. GFLOPS and MFLOPS are common units of measurement.

• Bandwidth: A metric signifying the total amount of data which can be handled per unit of time. It can also be measured by profiling tools like Intel Advisor and nvprof. gigabytes/sec, megabytes/sec are widely used as units of measurement.

• Particle Tracking Rate: A particular metric targeting to describe the performance of any particles codes such as MC codes. Its unit of measurement is particles/sec.

• Speedup: A metric representing a ratio of a measured value to a baseline value on the platform. This metric has no unit.

• Application Efficiency: A metric describing the achieved performance as a fraction of the best observed performance on the platform [START_REF] Satish | Can Traditional Programming Bridge the Ninja Performance Gap for Parallel Computing Applications?[END_REF].

• Architectural Efficiency: A metric representing the ratio of achieved performance to the peak theoretical performance on the platform.

Metrics of Productivity

The term "productivity" is a broad concept involving many aspects in HPC community. One may be of interest to assess the productivity of a HPC center by measuring the ratio of scientific output to total costs of ownership [START_REF] Wienke | Productivity and Software Development Effort Estimation in HPC[END_REF] [START_REF] Wienke | Productivity and Software Development Effort Estimation in High-Performance Computing[END_REF] or to measure the productivity of an individual by computing the performance of an application that he has been devoted to developing as a fraction of his development effort [START_REF] Zelkowitz | Measuring Productivity on High Performance Computers[END_REF]. When it comes to the productivity of an application, it may be assessed by the LOC (lines of code) required by this application which adopts a particular programming model or algorithm to achieve an expected behavior.

Due to the confounding factors of productivity, the correctness of flattening it into a metric has been challenged [START_REF] Jaspan | No Single Metric Captures Productivity[END_REF] and it seems that to make a metric of productivity more convincing, one should narrow the concept of productivity down as much as possible, only focusing on a single factor to study.

The first level of classification is based on the subject. Productivity of HPC's sites, developers, applications differentiates among each other. For example, we use LOC as a measurement of productivity. For a developer, it's a subjective metric since LOC highly depends on his skillset. But for an application, LOC is quite objective because the basic mechanisms required by this application of a specific programming model or algorithm to perform an expected behavior remain almost unchangeable.

The second level of classification divides productivity into explicit productivity and implicit productivity where explicit productivity bridges the relat on between output and input with a formula like output input while implicit productivity only focuses on the judgement of input. The input may be development effort and the output is usually considered as performance.

In our case, we firstly study the implicit productivity of an application a solving a problem p on a platform h, (a, p, h). The definition is given by: Definition 5.1.14 Implicit productivity of a case (a, p, h) refers to the development effort devoted to a for solving a problem p on a platform h.

Note that the time dependence of development effort is ignored because this quantity is heavily influenced by developer's skillset. Thus, from the perspective of applications, time should not be considered as a parameter in the metrics of implicit productivity.

The typical metrics of implicit productivity are LOC and NW (number of words), which has already been used in many productivity studies [117][83]. Since the code size of an application is relevant to the development style of developers, comparisons of LOC or NW among a set of applications should be done under the condition that all applications are developed in a same coding style, adopting same algorithms.

We propose another metrics of implicit productivity, PS-LOC (platform-specific lines of code) and PS-NW (platform-specific number of words) which are used to measure the lines of code or number of words of a platformspecific code path for an application which allows it executing correctly on a particular platform.

The Halstead complexity metric [START_REF] Halstead | Elements of Software Science[END_REF] is also a choice to characterize implicit productivity of a case (a, p, h). Besides, Harrell et al. [START_REF] Harrell | Effective Performance Portability[END_REF] proposed a metric D(A) called "code divergence" to quantify the difference of LOC between distinct applications targeting different platforms in A. The formula is given by:

Instead of counting lines

D(A)= 0 @ |A| 2 1 A 1 P {ai,aj }⇢A d (a i ,a j ) d (a i ,a j )= |LOC(a i ) LOC(a j )| min(LOC(a i ,LOC(a j )) (5.21)
where d (a i ,a j ) represents the pairwise distance which is quantified by normalized change in LOC between the application a i targeting to a particular platform and the application a j targeting to the same or different platform.

In fact, the pairwise distance d (a i ,a j ) can be used to calculate the corresponding distance of an application between its different platform-specific code paths. We have:

d i,j (a, p, h)=d [(a, p, h i ), (a, p, h j )] = |LOC(a, p, h i ) LOC(a, p, h j )| min(LOC(a, p, h i ),LOC(a, p, h j )) (5.22)
Furthermore, for applications implemented by high level programming model, their code paths of different platforms may differ little between each other despite some addition and removals of lines. For example, an application implemented by OpenACC may use a single line to set device type, such as acc_device_nvidia and acc_device_host. In order to change its device type from Nvidia GPU to host, it is required to delete the line of acc_device_nvidia and add a new line targeting to acc_device_host. The total number of LOC of two code paths remain the same. Therefore, for the purpose of characterizing this feature, a measure called churn is defined [START_REF] Harrell | Effective Performance Portability[END_REF]:

chn i,j (a, p, h)= # lines of add/del l l = 8 < : 1, ∆LOC =0 ∆LOC ∆LOC 6 =0 ∆LOC = |LOC(a, p, h i ) LOC(a, p, h j )| (5.23)
Several relative metrics of measurements mentioned above may also be specified with respect to baseline values as reference, such as relative LOC, relative program effort and so on. Overall, we have a set of implicit productivity metrics described in the following equation: As for explicit productivity, it is in fact a hybrid definition of "performance-implicit productivity" which is described as follow:

C(a, p, h)= 8 > > > > > > > > > > > > > < > > > > > > > > > > > > > : d i,j (
Definition 5.1.15 Explicit productivity of a case (a, p, h) refers to the ratio of the performance of an application a solving a problem p on a platform h to the corresponding development effort devoted to a.

A classical measure of explicit productivity is called RDTP (relative development time productivity), which is shown in Equation 5.25. This metric is originally developed to evaluate parallel code development where the relative speedup and effort are ratios of performance/development effort of parallel code to that of serial code [START_REF] Funk | Application of a Development Time Productivity Metric to Parallel Software Development[END_REF].

Ψ relative = relative speedup relative effort (5.25)

As a matter of fact, the reference version required to compute relative speedup and relative effort in Equation 5.25 may not limited to a serial application. In our scenario, we may set a parallel code as reference, such as OpenMP version, and calculate relative metrics of a set of parallel applications implementing different programming models.

Furthermore, the relative effort in Ψ relative is mainly characterized by LOC, although other attributes can also be considered such as development time, function points [138][67]. We prefer some metrics of implicit productivity such as PS-LOC to account for development effort in Ψ relative and we eventually have:

C(a, p, h)= relative speedup relative metric of implicit productivity (5.26) It should be noted that if the value of relative speedup or metric of implicit productivity can't be zero, otherwise, C(a, p, h) becomes meaningless.

Metrics of Performance Portability

The metric P(a, p, H) proposed by Pennycook which is described in Equation 5.1 can actually be represented by the standard version of portability metric F hm applying to a data set C(a, p, h) retrieved by a performance efficiency The first type mainly applies to several performance metrics indicating the absolute performance, including time, throughput, bandwidth as well as particle tracking rate. Because the absolute performance is unable to reflect the scale of performance on a particular platform, we can use the percentage of change to implicitly express the scaled contribution on this platform.

On the contrary, the second type is for metrics measuring the relative performance such as application efficiency and architectural efficiency. Since these metrics represent itself a scaled value on a particular platform, we can directly use the amount of change to present the scaled contribution on this platform. Overall, with respect to Property 5.1.5 F hm is the least adaptive metric in F to evaluate performance portability since the percentage of change for F hm is not platform-independent in all cases and in turn fails to meet the requirement described in Definition 5.1.17.

Strict and Relaxed Performance Portability

One of key features of P is that it highlights the property of portability which can be expressed as a boolean value.

In this way, if an application is not supported by all platforms in H, then the metric value will equal to zero. Similarly, F gm and F hm applying to a data set C(a, p, h) retrieved by a performance metric also keep this feature. We call the performance portability maintaining this feature a strict performance portability, which can be defined as:

Definition 5.1.18 Strict performance portability (SPP) inherits the entire feature of portability saying that for an application a solving a problem p, if this application is not supported by all platforms in the platform set H, the estimated value will be zero.

However, strict performance portability addresses significantly the state of portability at the cost of lack of expressiveness for the metric to indicate which application can hit the highest score value in terms of performance portability when a large number of applications may not be supported by all platforms in H.

For example, in Table 5.2, if we calculate the strict performance portability of a set of applications {a 1 ,a 2 ,a 3 ,a 4 } solving a problem across a platform set H = {h 1 ,h 2 ,h 3 ,h 4 ,h 5 ,h 6 }, the estimated values of all applications will be zero since each application is not supported by all platforms in H. In fact, we can't retrieve non-zero values of {a 1 ,a 2 ,a 3 } and {a 4 } at the same time. Moreover, in order to obtain non-zero values of {a 1 ,a 2 ,a 3 }, the maximal platform set is limited to H = {h 3 ,h 5 }. Such behavior is impractical and requires users to cautiously choose the platform set so as to attain comparable results from metrics of performance portability. After all, it is meaningless that most of applications gain the same score, zero. In pursuance of exploiting complete information, it is obligated to apply the metric of strict performance portability to a variety of sub-sets of platforms which makes the procedure of evaluation redundant and complicated. With respect to Definition 5.1.18 and Definition 5.1.19, we know that F am is a natural metric calculating relaxed performance portability while F gm and F hm are metrics aggregating strict performance portability. We need to redesign F gm and F hm in order to make them capable of characterizing relaxed performance portability. But for convenience, we make modifications of F gm and F hm as well as F am .

Variants of F

Generally, the criterion for design of variants of F can be summarized as follows:

1. Quantify both the performance and portability.

2. The metric equals to zero when all platforms in the set H do not support (a, p).

3. With the increase of non-supported platforms for an application, the metric value degrades. Here is our proposition of a transfer operator T 1 to get the first variant of F :

(T 1 F ) • C(a, p, h)= |S| |H| F 0 (5.28)
where |S| is the total number of supported platforms and |H| is the total number of platforms. Based on Equation 5.20, we have T 1 F am = F am .

T 1 F directly use the estimated value of F applying to the data set retrieved from all supported platforms in the platform set H and add a multiplier factor, the ratio of total number of supported platforms S to the total number of platforms H, to weight it. Alternatively, we can also treat F 0 as a part of terms for all elements in C(a, p, h) which are retrieved from non-supported platforms and we have:

T 2 C(a, p, h)={C 0 i (a, p, h)} 8i=1,2,...,|H| C 0 i (a, p, h)= 8 < : f (F 0 ) if C i (a, p, h)=0 C i (a, p, h) if C i (a, p, h) > 0 (5.29)
where T 2 is the second transfer operator applying to the data set C(a, p, h). f is a function parameterized by F 0 with a positive value as output. F am , F gm and F hm have its corresponding function f to meet the third criteria proposed above, saying that with the increase of non-supported platforms for an application, the metric value decreases.

Thus, we can calculate their critical factors k am , k gm , k hm to learn about the maximal and minimal values of f .

Let S and N the sub-sets of a platform set H containing respectively the supported and non-supported platforms in H for a couple (a, p), K 1 =

P |S| i=0 Ci(a,p,h) |S| , K 2 = Q |S| i=1 C i (a, p, h) 1 |S| , K 3 = |S| P |S| i=0 1 C i (a,p,h)
. The criticality equations of F am , F gm and F hm are described in Equation 5.30:

P |S| i=0 C i (a, p, h)+|N |k am |S| + |N | = K 1 |S| Y i=1 C i (a, p, h) 1 |S|+|N | |N | Y i=1 k 1 |S|+|N | gm = K 2 |S| + |N | P |S| i=0 1 Ci(a,p,h) + |N |k hm = K 3 (5.30)
And we have:

k am = F 0 am k gm = F 0 gm k hm = 1 F 0 hm (5.31) • if f (F 0 am ) <k am : F am • (T 2 C(a, p, h))
decreases with the increase of non-supported platforms for the couple (a, p).

• if f (F 0 gm ) <k gm : as above.

• if f (F 0 hm ) >k hm : as above.

With the intention of highlighting the effect of the increase of non-supported platforms on F , we choose f (F

0 am )= kam |N | , f (F 0 gm )= kgm |N | , f (F 0 hm )= |N |
k hm for numerical evaluations afterward. At last, we propose the third transfer operator along with the third variant of F described as follow:

(T 3 F ) • C(a, p, h)=Integer[|S|].Decimal[ f F 0 ] (5.32)
where Integer[|S|] means that we use |S| as the integer part of metric and Decimal[ f F 0 ] as the decimal part of metric.

f F 0 is the normalized value of F 0 with f F 0 = F 0 F 0 max . F 0 max
is the maximal F 0 obtained among a set of applications. Note that if f F 0 equals to 100%, the decimal part will be noted as 0. 9.

T 3 F is a more intuitive formula which addresses more portability comparing to the other ones. It allows numerical comparison of performance portability among a set of supported/non-supported platforms where the performances of different applications are only compared when these applications have the same number of supported platforms.

Overall, the variants of F that we adopt enable the numerical comparison of performance portability among applications across different architectures when not all applications are supported by all architectures. This feature is more adaptive under the circumstances that the supported platforms for each application in an application set differentiates significantly among each other, resulting in small application and platform sets.

Metrics of Productivity Portability

As the definition of performance portability, we give the definition of productivity portability:

Definition 5.1.20 Productivity portability is a measurement adopting various aggregating methods which estimates the average productivity of an application a solving a given problem p across a platform set H with the intention of comparison.

In this way, F becomes the generic metric of productivity portability when it applies to a data set retrieved by a productivity metric.

Metrics of Implicit Productivity Portability

When F applies to implicit productivity metrics, it characterizes the average development effort for a couple (a, p) across a set of platforms. The metrics that indicate platform-specific code paths of an application should be adopted, such as PS-LOC. Since there is no meaning to aggregate the PS-LOC of a non-supported platform, the platform set should be limited to a sub-set of platforms supporting this application.

The adaptability of F towards metrics of implicit productivity portability (IPP) is unlike F towards performance portability metrics. Since metrics such as LOC, PS-LOC is not scaled on each platform, we are unable to choose the aggregating method (arithmetic, geometric, harmonic means) with reference to the platform-independent percentage of change. By contrast, we give all aggregating methods a chance to be examined.

Moreover, the metric of pairwise distances d i,j (a, p, h) may also be applied to F and we obtain the code divergence of the application of different platform-specific code paths as described in Equation 5.33 based on Equation 5.22. The number of combinations for taking r things out of n things is denoted as Comb(n, r).

F • C(a, p, h)= 8 > > > > > > > > > > > > < > > > > > > > > > > > > : F am • C(a, p, h)= P {i,j}⇢H d i,j (a, p, h) Comb(|H|, 2) F gm • C(a, p, h)= Q {i,j}⇢H d i,j (a, p, h) 1 Comb(|H|, 2) F hm • C(a, p, h)= Comb(|H|, 2) P {i,j}⇢H 1 d i,j (a, p, h) (5.33)
Similarly, according to Equation 5.23, we may apply F towards the metric of churn chn i,j (a, p, h), as described in Equation 5.34 where chn i0,i refers to the churn distance between the reference code path i 0 and another code path i for an implementation. This value expresses the average addition and removals of lines required for an implementation of a given coda path to be ported to another platform in a given platform set. Particle tracking rate (particles/s), the higher the better; Application Efficiency (%), the higher the better.

and application should be reset to smaller ones where all applications in the application subset are supported by all platforms in the platform subset. Although we may draw graphs of performance portability varying over a sequence of subsets of platforms, as what Deakin has done in his paper [START_REF] Deakin | Performance Portability across Diverse Computer Architectures[END_REF], this requires a large amount of work for repeated calculations of performance portability over a range of architectures. As an alternative, we have used the variants of F to retrieve the values of performance portability.

All results retrieved by different aggregating methods are calculated as reference. According to Property 5.1. In conclusion, the results listed in Table 5.5 and Table 5.6 prove that F hm 6 F gm 6 F am . For example, the values retrieved by the variant T 1 F hm , T 1 F gm , T 1 F am for the CUDA version in Table 5.5 increase from 17.0 ⇥ 10 2 particles/s to 18.8 ⇥ 10 2 particles/s. Excluding the variant of F transferred by the operator T 3 which maintains its own scoring mechanism, the CUDA version aggregates the highest value of particle tracking rate and the OpenACC version gains the highest application efficiency when they are evaluated by the metric F as well as its variants with the utilization of two transfer operators T 1 and T 2 . As for the metric T 3 F , the OpenACC version obtains the best performance in terms of particle tracking rate and application efficiency, since this variant addresses much more portability than performance, which makes the OpenACC and Kokkos versions the most competitive implementations of slabAllNuclides running across the given ten platforms in terms of performance portability. where it is obliged to define the Kokkos execution space during the compilation. The differences of code paths for each implementation are summarized as follows:

Productivity Portability Evaluation

Implicit Productivity Portability

• OpenMP thread: none.

• CUDA: none.

• OpenACC: set acc_device_host or acc_device_nvidia.

• OpenMP offload: nb_devices = omp_get_num_devices(), the program automatically decide whether it is executed on host or device according to nb_devices. Considering the evaluation of explicit productivity portability (EPP or productive performance portability), we first need to choose a subset of the original platform set H so as to make at least one application in the application set capable of running across the platforms in the subset of H. Fortunately, we may exclude the platform of the machine Intel NUC associated with an Intel GPU and reset the platform set to H 0 with |H 0 | =9. In this way, both the OpenACC and Kokkos versions are good candidates for reference and the relative speedup and implicit productivity (PS-LOC) of different programming models are illustrated in Table 5.9 as we use the results of the OpenACC version as baseline values.

Explicit Productivity Portability

Based on the values of Ψ relative listed above, we can use the variants of F to evaluate the productive performance portability of different implementations of slabAllNuclides across the platform set H 0 . Note that the blank cells in Table 5.9 signify the meaninglessness and cannot be replaced by zero values. Thus, the original metric of F am is also unable to aggregate meaningful results.

From Table 5.10, we find that F am , F gm and F hm retrieve the same priority order. According to Property 5.1.4, F am is the most adaptive metric for the evaluation of relative value, so we give an explication of results obtained by F am which are highlighted with green color. Overall, excluding the metric T 3 F • C(a, p, h) which addresses more portability than productive performance, the results shows that the OpenMP thread version obtains the best productive performance across a set of platforms Chapter 5

Metrics Comparison among Programming Models 132 which makes it the best candidate for implementation in terms of performance, portability and productivity with the given platform set. The OpenACC and Kokkos versions are more competitive to the CUDA version in terms of 3P which are mainly due to their much wider range of portable platforms and much higher abstraction of code.

The OpenMP offload version gains bad scores since it achieves bad performance and covers a narrow range of platforms in the platform set. Lastly, the SYCL obtains bad results however these values may be biased since this programming model is preliminary and its compiler support is underdeveloped. In order to retrieve unbiased estimates, we need carry out more tests across more architectures for the SYCL version with the use of the next version of SYCL compilers.

It should also be addressed that the conclusions drawn above depend on the applications, problems and platforms (environments included) that we use in the thesis [31][32]. For example, it is safe to say that the OpenMP thread version of our Monte Carlo neutron transport codes obtains the best productive performance across a set of platforms to solve a fixed source Monte Carlo problem using SIGMA1 on-the-fly Doppler broadening approach.

However, if we change the applications, the problems, or the platforms, other conclusions may be reached instead.

Hence rather than giving a guide for users who want to make a choice of programming model and environment for their applications with the use of our generic metric, our interest lies in offering users a practical tool to analyze their special cases and help them making their own choice.

Chapter 6

Conclusion and Future Work

Conclusion

This thesis project depicts the proposal of a generic metric as well as its variants to evaluate Monte Carlo neutron transport code in terms of performance, portability and productivity through the implementation of a benchmark based on the Monte Carlo neutron transport prototype PATMOS via different programming languages or libraries.

The development work is based on the SIGMA1 on-the-fly Doppler broadening method which calculates cross section data on-the-fly during the simulation. Because this method has added more FLOP workloads and thus manages to mitigate the memory-bound issue caused by random memory accesses which are mainly incited by the binary search of the typical pretabulated method. Besides, we adopt a heterogeneous offloading strategy for our portable implementations of Monte Carlo codes, which performs all the particle tracking and scoring on host, and offloads the calculations of total microscopic cross sections of nuclides to accelerators since the cross section computation accounts for up to 95% of total run time in PATMOS. To achieve this heterogeneous offloading strategy, the main development effort is to rewrite codes which are dedicated to porting nuclide data and microscopic cross section calculations to devices via implementations of several programming models.

The first part of work in this thesis project focuses on the portable implementations of microscopic cross section computational kernel in PATMOS. Our general programming model can be abstracted as OpenMP thread + {X}, where {X} can be any languages or libraries which are capable of parallel programming on modern accelerators including CUDA, OpenACC, OpenMP offload, Kokkos, SYCL or none at all. For the purpose of reducing the large amount of overheads introduced by kernel launch and memory transfer operations with the utilization of conventional particle tracking approach, the history-based method, we have proposed another approach, called the pseudo eventbased method, which reorganizes the procedure of history tracking into three events which are respectively assigned to handle data banking, microscopic cross section calculation and the rest of processes (macroscopic cross section summarization, distance sampling, particle movement and interaction).

The tests show that the pseudo event-based method generally exceeds a magnitude of performance achieved by the history-based method. It corresponds to the fact that the pseudo event-based method manages to batch small transfers into one larger transfer and merge small kernels into one larger kernel which eliminates significantly the data transfer and kernel launch overheads and exploits higher computational capability of accelerators. The numerical quantities such as the achieved occupancy have also proved that the pseudo event-based method achieves to improve the compute throughput and memory bandwidth of Nvidia GPUs in comparison with the history-based method. However, from the perspective of the utilization of compute resources of Nvidia GPUs, both the pseudo event-based and history-based methods are far from making full use the compute resources. The profiling analysis indicates that our heterogeneous offloading code may be mitigated by reducing the number of registers used for each CUDA thread, adding instruction-level parallelism and optimizing memory alignment and access patterns.

As for the comparison of performance achieved by different portable implementations of slabAllNuclides, the results show that in host mode, the OpenMP thread and Kokkos versions which can be compiled by Intel compiler with the vectorization options enabled gain better performance than other implementations. In offload mode, the CUDA and versions always renders better performance than other implementations via both the pseudo event-based and history-based methods.

The second part of work in this thesis project intends to evaluate our portable implementations of Monte Carlo neutron transport codes in terms of performance, portability and productivity. A generic metric F has been established based on the metric of performance portability P(a, p, H) proposed by Pennycook and his group. Comparing to Pennycook's metric, our generic metric allows for the aggregation of values retrieved from an application set or a problem set. Besides, it integrates three aggregating methods which are respectively the arithmetic mean, the geometric mean and the harmonic mean. The adaptabilities of three aggregating methods applying to different types of data set have been analyzed theoretically. According to the definitions of performance, we have proved that F gm is the most adaptive metric applying to a data set retrieved by a metric indicating the absolute performance and F am is the most adaptive metric applying to a data set retrieved by a metric indicating the relative performance.

Furthermore, we have proposed several variants of F to make it capable of evaluating relaxed performance portability where the estimated value will be zero only if the given application is blocked by all platforms in the platform set. Such variants make users easier to obtain meaningful values without the obligation of finding the appropriate application or platform set where each application in the application set is supported by all platforms in the platform set.

We have used the generic metric F • C(a, p, h) and its variants to evaluate the performance portability and productive performance portability of different implementations of slabAllNuclides (OpenMP thread, CUDA, OpenACC, OpenMP offload, Kokkos and SYCL) across a set of platforms.

The main contributions of this research work can be concluded into several points:

• Proposition of a new particle tracking algorithm, pseudo event-based method for the heterogeneous offloading strategy of SIGMA1 on-the-fly Doppler broadening method.

• Implementations of SIGMA1 on-the-fly Doppler broadening method via multiple programming models.

• Proposition and analysis of the generic metric and its variants F • C(a, p, h) for the quantitative evaluation of portable implementations of Monte Carlo neutron transport codes in terms of portability, performance and productivity.

Future Work

There is some work that we intend to do for future development. From evaluation metric side, the operator F need to be testified by much more different cases which are widely used in the domain of scientific computing. More aggregating methods may also be integrated into our generic metric to enrich the capability of F handling complex problems. Besides, more tests need to be done so as to cover a wider range of architectures (AMD, Intel GPUs, etc.)

for the evaluation of our portable implementations of Monte Carlo neutron transport code in terms of performance, portability and productivity. We also intend to find an optimal version using portable programming model (on Nvidia GPUs) such as the OpenACC version and compare the cost the portability for moving it to other architectures.

From programming model side, we have interest to optimize the OpenMP offload, Kokkos and SYCL versions.

For example, Kokkos may enable the support of MTMA strategy in future which makes it possible to allow multiple threads targeting to multiple devices at the same time. Several implementations of SYCL such as DPC++ may completely support Nvidia GPUs. All of these possible changes need to be tested and integrated in our portable implementations of Monte Carlo neutron transport code.

From Monte Carlo neutron transport simulation side, since we have demonstrated that porting total microscopic cross section calculation with Doppler broadening techniques to accelerators may contribute to a significant performance improvement, we can offload partial cross section calculation as well in order to make our implementation more adaptive to other complex cases. Moreover, we may also offload other methods of cross section calculation such as pre-tabulated and multipole methods to accelerators. In this thesis we propose to evaluate the different programming models available for addressing manycore and / or heterogeneous architectures within the framework of the Monte Carlo transport codes. A simple but representative application test case will be considered in order to cover a fairly wide range of so-lutions and compare them in terms of performance, portability of performance, ease of implementation and maintainability. The target architectures are 'classic' CPU, Intel Xeon Phi and GPU. The most relevant programming models will then be set up in a Monte Carlo transport code.

vii 5 . 4

 54 Metrics of performance of different implementations of slabAllNuclides across a set of platforms. . . . . . . . . . .

5. 7

 7 Metrics of implicit productivity for implementations of slabAllNuclides across a set of platforms. . . . . . . . . . . . 5.8 Implicit productivity portability evaluated by F for implementations of slabAllNuclides across a set of supported platforms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Metrics of explicit productivity for implementations of slabAllNuclides across a set of platforms. . . . . . . . . . . . 5.10 Evaluation of productive performance portability for different implementations of slabAllNuclides across a subset of platforms via the generic model F , Ψ relative , the higher the better. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 1 :

 21 Figure 2.1: Von Neumann Architecture.

Figure 2 .

 2 Figure 2.2 show the workflow of instruction and data stream during program execution for the computer. Instruction stream refers to the flow of instructions fetched from main memory and operated by CPU. Data stream is the flow of operands loaded or stored between main memory and CPU. The multiplicity of stream has two states, single or multiple, indicating that if there is only one stream or multiple streams during any one clock cycle.

Figure 2 . 2 :

 22 Figure 2.2: Instruction and Data Stream.

2. 1 . 2 . 1

 121 SISD organization SISD computing machines are serial computers which handles only one stream of instructions and data at any point in the execution. Some old type of computers such as CRAY1 and IBM360 adopt this organization. The diagram of SISD organization is illustrated in Figure 2.3.

Figure 2 . 3 :

 23 Figure 2.3: SISD Organization.

  .

  4.

Figure 2 . 4 :

 24 Figure 2.4: SIMD Organization.

Figure 2 . 5 :

 25 Figure 2.5: MISD Organization.

  .6.

Figure 2 . 6 :

 26 Figure 2.6: MIMD Organization.

Figure 2 . 7 :

 27 Figure 2.7: Shared Memory Architecture.

Figure 2 . 8 :

 28 Figure 2.8: Distributed Memory Architecture.

Figure 2 .

 2 Figure 2.11 denotes the die diagram of a Broadwell microarchitecture with 10 cores. As we can see, the diagram can be divided into two main areas: Uncore and Core. Uncore is comprised of multiple modules designed specifically for the server and workstation market space [62], such as LLC (Last Level Cache), Cbo (Caching Agent, including Cache Box and Core Box), QPI (Intel QuickPath Interconnect), HA (Home Agent), iMC (Integrated Memory Controller), IIO (Integrated I/O module), PCU (Power Control Unit) and UBox (Config Agent). Core contains mainly the cores executing instructions and operations as well as two levels of caching hierarchy (L1 and L2).

Figure 2 .

 2 Figure 2.12: POWER9 SMT Core Diagram [7].

Figure 2 .

 2 Figure 2.13: POWER9 SMT4 Chip Diagram where the cache hierarchy is illustrated [108].

Figure 2 .

 2 Figure 2.12 shows the building blocks of POWER9 core. As we can see, a SMT core is composed of 1 or 2 super-slices (SMT4 or SMT8), along with an instruction fetch unit (IFU) and an instruction sequencing unit (ISU).

Figure 2 .

 2 Figure2.15 shows the main components of Intel Gen9.5. Overall, Intel Gen9.5 is comprised of three parts: unslice, display, as well as slice. Among them, slice is responsible for calculations. It is a cluster of subslices. Each subslice consists of 6 or 8 EUs which are dedicated to execute 3D shaders and computational kernels. Based on different scaling models, Intel Gen9.5 is allowed to have maximum 72 EUs and minimum 12 EUs. Each EU is

Figure 2 . 15 :

 215 Figure 2.15: Intel Gen9.5 Diagram [5].

Figure 2 . 16 :

 216 Figure 2.16: Nvidia Pascal GP100 Streaming Multiprocessor Unit [87].

  , each SM is composed of a L1 instruction cache, a 128KB combined L1 data cache and shared memory subsystem where shared memory size can be set up to 96KB, 4 texture units as well as 4 processing blocks. Each processing block is comprised of 16 FP32 cores, 8 FP64 cores, 16 INT32 cores, 2 mixed-precision tensor cores (capable of handling FP64, FP32, FP16 operations), a L0 instruction cache, a warp scheduler, a dispatch unit and a 64KB register file. In total, a Tesla V100 contains 5120 FP32 cores, 2560 FP64 cores, 5120 INT32 cores, 640 tensor cores, 320 texture units and 20MB register file size.

Figure 2 .

 2 Figure 2.18 show the CUDA thread hierarchy. Theoretically, CUDA provides the organization of threads in three dimension. But in practice, a grid is usually organized as a 2D arrays of blocks and a block is composed of a 3D array of threads.

Figure 2 .

 2 Figure 2.18: CUDA Thread Hierarchy [33].

Figure 2 . 20 :

 220 Figure 2.20: SYCL Implementations [8].

Figure 3 . 2 :

 32 Figure 3.2: Effect of temperature on cross sections [36].

Figure 3 . 3 :

 33 Figure 3.3: Neutron histories tracked from their birth to death through a set of volumes where green cycles refer to active neutrons and red cycles refer to desactived neutrons.

Figure 3 . 4 .

 34 It is mainly comprised of three steps including initialization, particle tracking and finalization. The initialization step is responsible for setting up materials and generating particles from source. The particle tracking step is devoted to tracking particle histories one by one in which the cross section calculations account for the most consuming part of Monte Carlo neutron transport simulation. The finalization step is in charge of tally computation in order to retrieve estimates of physical quantities.

Figure 3 . 4 :

 34 Figure 3.4: Serial Monte Carlo neutron transport simulation workflow.

1 foreach particle distributed to the processing thread do 2 while particle is alive do 3 calculation 4 • 5 • 6 • 8 if

 234568 of macroscopic cross section, with several vectorizable parts; outer loop for calculating macroscopic cross section; inner loop of energy lookup to find E low and E up ; inner loop of integral injected by on-the-fly Doppler broadening method; 7 sample distance to collision in material; new position is still inside material then 9 move particle to new position, and do collision;

Algorithm 3 : 2 foreach active particle do 3 4 end 5 foreach vector of event type do 6 perform

 323456 Event-based Monte Carlo transport algorithm.1 while any particles are alive on the processing thread doIdentify next event type and move it to the vector of corresponding event type; event for all particles in vector;

3 foreach 13 foreach 17 foreach

 31317 particle in the vector of migration do 4 calculation of macroscopic cross section; 5 sample distance to collision in material; 6 if new position is still inside material then 7 move particle to new position; 8 shuffle particle to the vector of corresponding collision; particle in the vector of scattering do 14 do scattering; 15 shuffle particle to the vector of migration; 16 end particle in the vector of absorption do 18 do absorption; 19 add produced particles (due to fission) to the vector of migration;

5 )⇡ Z 1 a e z 2 dz ( 4 . 6 ) 8 : 1 =0 7 += 4 Portable

 51468174 where erf c and exp are respectively the complementary error function (see Equation 4.6) and the base e exponential function. erf c (a)= 2 p Algorithm Function of compute_integral(). Input: a given energy E, a coefficient ↵, a lower-upper bound lb, ub and an energy grid eg[N ], a cross section grid csg[N ] Result: accumulated cross section after integral computation of the first element */ 2 index = lb, ea = eg[index],eb = eg[index + 1]; 3 x = p ↵ • ea, a = x y; 4 calculate F 0 (a), F 1 (a), F 2 (a), F 3 (a), F 4 (a); /* According to Equation 4.5 */ 5 calculate A and B based on F n (a); /* According to Equation 4.4 */6 siab = csg[index + 1] csg[index] ↵(eb ea) ; A(csg[index] siab • ↵ • ea)+B • siab;/* Accumulation from the second element to the second last element */ 8 for index lb +1to ub 1 by 1 do9 ea = eg[index 1],eb = eg[index],ec = eg[index + 1]; x = p ↵ • eb, a = x y; calculate F 0 (a), F 1 (a), F 2 (a), F 3 (a),F 4 (a); /* According to Equation 4.5 */ calculate A and B based on F n (a); /* According to Equation 4.4 */ siab = csg[index] csg[index 1] ↵(eb ea) ; sibc = csg[index + 1] csg[index] ↵(ec eb) ; += (sibc siab)(B A • ↵ • eb); end /* Accumulation of the last element */ index = ub, ea = eg[index 1],eb = eg[index]; x = p ↵ • eb, a = x y; calculate F 0 (a), F 1 (a), F 2 (a), F 3 (a), F 4 (a); /* According to Equation 4.5 */ calculate A and B based on F n (a); /* According to Equation 4.4 */ siab= csg[index] csg[index 1] ↵(eb ea) ; = A(csg[index] siab • ↵ • eb)+B • siab;Based on Equation4.3, we have the algorithm of the function compute_integral() with a given energy, lower and upper bounds of the interval of integration, a coefficient ↵ = M k B (T T0) and cross section grids as inputs (seeAlgorithm 8).Chapter Implementation of on-the-fly Doppler Broadening inPATMOS 57 The calculations of F n (a) fucntions call erf c and exp functions multiple times which accounts for a large part of the computing time. A CPU-based profiling of a typical Monte Carlo neutron transport benchmark, slabAllNuclides, implemented in PATMOS in terms of run time percentage is illustrated in Table4.1 in which all results were retrieved by C++ native clock and perf (more detail about slabAllNuclides is elaborated in section 4.4.1).

  . The derived classes of Nuclide are implemented via different programming languages or libraries to accomplish the nuclide data transfer from host to devices. The UML (unified modeling language) diagram of class Nuclide along with its derived classes is depicted in Figure 4.1.

Figure 4 . 1 :

 41 Figure 4.1: Inheritance relationship between base class Nuclide and its derived classes.

Figure 4 . 2 :

 42 Figure 4.2: Inheritance relationship between base class NeutronMedia and its derived classes.

Figure 4 . 3 :

 43 Figure 4.3: Inheritance relationship between base class NeutronMediaNavigator and its derived classes.

Figure 4 .Algorithm 9 : 2 total_sigma = 0; 3 for 5 } 6 for 12 E1= 25 }

 4923561225 3 shows the diagram of the class NeutronMediaNavigator. Pseudo-code of OpenMP thread implementation.1 void NeutronMediaNavigator::computeSigma(...){ (int t=0; t<N; ++t){ /* calculate microscopic cross sections of N nuclides */ 4 sigma1DopplerBroadening(t, h_macroscopic_cs, ...); (int t=0; t<N; ++t){ /* calculate macroscopic cross section of N nuclides */ 7 total_sigma += h_macroscopic_cs[t] * h_concentrations[t]; ...;E 2=...; 13 lb = getLowerBound(E1, ...); ub = getLowerBound(E2, ...)+1 ; 14 h_macroscopic_cs[i] += compute_integral(lb, ub, ...);In general, the parallelization effort is made with the function sigma1DopplerBroadening(), as a part of the function computeSigma() in the class NeutronMediaNavigator as well as its derived classes to vectorize the SIGMA1 on-the-fly Doppler broadening algorithm on host or accelerators. Furthermore, since the function getLowerBound() takes up a rather negligible run time percentage comparing to that of the function compute_integral(), we basically focus on the parallelization of compute_integral(), more specifically, the parallelization of the for loop in the integral computation.

  coded in the functions allocateN ucOnDevices() and deleteN ucOnDevices() of derived classes inherited from the base class Nuclide. Furthermore, the functions allocateN ucDataOnDevices() and deleteN ucDataOnDevices() of child classes inherited from parent class NeutronMedia perform a for loop to recursively allocate, copy and remove data of all nuclides. The microscopic cross section calculations are offloaded by parallelizing the for loop in the function computeSigma() of derived classes inherited from the base class NeutronMediaNavigator (Algorithm 9 at lines 3-5).The parallelization is achieved by distributing each call of a function sigma1DopplerBroadening() in the outerloop to a work-group of device where each work-group is composed of a set of work-items which are responsible for vectorizing the inner-loop in the function compute_integral(). For example, the default strategy for Nvidia GPUs is to assign the workload of microscopic cross section calculation of each nuclide to a CUDA warp so as to achieve a tasklevel parallelism. And all CUDA threads within a warp are used to vectorize the inner-loop of integral computation as a completion of data-level parallelism. Our reference implementation on GPU achieved by CUDA is shown in Algorithm 10.The functions allocateN ucOnDevices() and deleteN ucOnDevices() in the class CUDANuclide calls the functions of the CUDA runtime API cudaM alloc(), cudaM emcpy() as well as cudaF ree() to allocate, copy and delete the data of a given nuclide on multiple devices (at lines 1-9). Note that the function cudaSetDevice() is used to set the current device associated to the corresponding host thread for kernel execution.

Algorithm 10 : 4 Portable

 104 Pseudo-code of CUDA implementation. void CUDANuclide::allocateNucOnDevices(...){ for(int idev=0; idev<nb_devices; ++idev){ /* allocate and copy the nuclide data to devices */ cudaSetDevice(idev); cudaMalloc(...); cudaMemcpy(..., cudaMemcpyHostToDevice); } } void CUDANuclide::deleteNucOnDevices(...){ for(int idev=0; idev<nb_devices; ++idev) /* remove the nuclide data from devices */ cudaSetDevice(idev); cudaFree(...); } void CUDANeutronMediaNavigator::computeSigma(...){ dim3 blocksize(32, 1); int gridsize = N; /* set CUDA kernel dimension */ cudaMalloc(...); cudaMemcpyAsync(..., cudaMemcpyHostToDevice, stream); kernel_launchngridsize, blocksize, 0, streamo(...); /* launch CUDA kernel */ cudaMemcpyAsync(..., cudaMemcpyDeviceToHost, stream); cudaFree(...); total_sigma = 0; ... /* macroscopic cross section calculation */ } __global__ void kernel_launch(...){ int t = blockIdx.x; /* ensure t < N */ sigma1DopplerBroadening(t, ...); } __device__ inline double compute_integral(...){ double sigma = 0; double sigma_tmp = 0; int index = lb; sigma += ...; for(int i=lb+1+threadIdx.x; i<ub; i+=32) sigma_tmp += ...; for(int mask=16; mask>0; mask/=2) sigma_tmp += __shfl_xor_sync(...); sigma += sigma_tmp; index = ub; sigma -= ...; return sigma; } As for the outer-loop in the function computeSigma() of class CUDANeutronMediaNavigator calling a function sigma1DopplerBroadening() during each iteration, a CUDA kernel kernel_launch() is replaced to parallelize this loop along with two asynchronous memory copies cudaM emcpyAsync() where all three functions are enabled with Chapter Implementation of on-the-fly Doppler Broadening in PATMOS 63 CUDA multi-stream support owning a single stream associated to the corresponding host thread (at lines 12-14). This CUDA asynchronous multistreaming technique allows for concurrent CUDA calls initiated by different host threads. By contrast, operations (data transfer and kernel calculation) with the same stream called by a single host thread are executed sequentially and cannot overlap. With reference to the computational kernel kernel_launch() (at lines 17-20), its default dimension is set to N threadblocks, where each threadblock contains one CUDA warp. Note that N is the gridsize of kernel (at line 11). Each warp is assigned to calculate the microscopic cross section of a single nuclide recorded by the task index t=b l o c k I d x . x . At last, the inner-loop in the function compute_integral() is vectorized by CUDA threads within a warp. Each CUDA thread handles a part of the workload and store the partial result in sigma_tmp. A sum reduction is then operated by resorting to a warp shuffle instruction __shf l_xor_sync() to aggregate the total value of sigma_tmp among CUDA threads within the warp.

19 }

 19 20 KOKKOS_INLINE_FUNCTION double compute_integral(...,K o k k o s : : T e a m P o l i c y < > : : m e m b e r _ t y p e &v e c t o r ) { 21 double sigma = 0; double sigma_tmp = 0; 22 int index = lb; sigma += ...;

  the function allocateN ucOnDevices() of the class SYCLNuclide. All nuclide data are stored in the cl::sycl::buffer class of one dimension as an input parameter of a command group which is submitted to a queue by invoking its member function submit(). The access of SYCL buffer is achieved by the class cl::sycl::accessor which provides an entry to data either on host or within a command group enqueued to SYCL queues. The constructor of SYCL queues adopted in our implementation accepts an instance of the class cl::sycl::default_selector which may select a host device or accelerators. The memory copy operation is undertaken by the member function copy() of the cl::sycl::handler class which copies the host data into the memory object accessed by SYCL accessors.The parallelization of outer-loop in the function computeSigma() is accomplished by the invocation of an hierarchical kernel which calls the function parallel_for_work_group() of the class cl::sycl::handler with two values of the class cl::sycl::range representing respectively the number of work-groups and the number of work-items per work-group as well as a user-defined functor encapsulating the function sigma1DopplerBroadening() in its operator() function. In our case, the number of work-groups is set to N where each work-group is assigned to a task of microscopic cross section calculation of a nuclide. The number of work-items per work-group denoted as work_group_size is hardware-dependent which needs developers to manually set its value. Besides, since by default all variables declared inside the parallel_for_work_group scope are allocated in local memory which are

Algorithm 14 :}

 14 Pseudo-code of SYCL implementation. void SYCLNuclide::allocateNucOnDevices(...){ cl::sycl::queue sq(cl::sycl::default_selector{}); /* create a SYCL queue of a particular device type */ cl::sycl::buffer<...> buffer_eg(h_eg, ...); /* wrap the host data in a SYCL buffer */ sq.submit([&] (cl::sycl::handler& cgh){ /* Enqueue a kernel */ auto d_eg_accessor = buffer_eg.get_access<...>(cgh); /* request access to the buffer */ cgh.copy(h_eg, d_eg_accessor); /* copy h_eg to the memory object accessed by d_eg_accessor */ sycl::private_memory<double> compute_integral(...,c l : : s y c l : : g r o u p < 1 > &w o r k _ i t e m s ) { size_t num_items = work_items.get_local_range(0); double[num_items] item_sigmas = {0}; cl::sycl::private_memory<double> sigma(work_items); cl::sycl::private_memory<int> index(work_items); work_items.parallel_for_work_item([&](cl::sycl::h_item<1> item){ sigma(item) = 0; sigma_tmp = 0; index(item) = lb(item); sigma(item) += ...; for(int i=lb(item)+1; i<ub(item); i+=item.get_get_local_range(0)) sigma_tmp += ...; item_sigmas[item.get_local_id(0)] = sigma_tmp; ) = ub(item); sigma(item) -= ...; Since it is possible to call the function parallel_work_for_item() of the class cl::sycl::group multiple times within the parallel_work_for_group scope, the inner-loop in the compute_integral() function can be vectorized by this function (at lines 25-26). However, due to the fact that SYCL of version 1.2.1 does not provide a pre-defined function performing the sum reduction, we are obligated to manually implement a reduction loop to summarize the final result with the utilization of shared memory (at lines 28-31). Although the latest version of SYCL released in 2020 offers the function parallel_reduce() of the class cl::sycl::handler along with the cl::sycl::reduction interface and cl::sycl::reducer class to perform reduction over all work items, this function has not been expanded

1 foreach particle distributed to the processing thread do 2 while particle is alive do 3 calculation of macroscopic cross section, with several vectorizable parts; 4 • 5 ⇤ 6 ⇤ integral computation; 7 •

 234567 do microscopic cross section lookups =) offloaded; binary search to find lower-upper bounds; sum up macroscopic cross section; 8 sample distance, move particle, do interaction; 9 end

Algorithm 16 : 2 while particles remain in bank do 3 foreach 5 end 6 • 7 ⇤ 8 ⇤ integral computation; 9 foreach 10 •

 16235678910 Heterogeneous offloading pseudo event-based Monte Carlo neutron transport algorithm. 1 foreach bank of particles distributed to the processing thread do remaining particle in bank do /* First event */ 4 bank required data for microscopic cross section lookups; do microscopic cross section lookups =) offloaded; /* Second event */ binary search to find lower-upper bounds; remaining particle in bank do /* Third event */ sum up macroscopic cross section; 11 sample distance, move particle, do interaction;

Figure 4 .

 4 Figure 4.4 depicts the run time per cycle (the lower the better) and strong scalability (the higher the better)

  #15305: vectorization support: vector length 2 remark #15301: OpenMP SIMD LOOP WAS VECTORIZED remark #15475: --begin vector cost summary -remark #15476: scalar cost: 450 remark #15477: vector cost: 195.000 remark #15478: estimated potential speedup: 2.280 remark #15488: --end vector cost summary --LOOP END Auto-vectorization enabled LOOP BEGIN at sigma1xse.cxx(88,3) remark #15305: vectorization support: vector length 4 remark #15475: --begin vector cost summary -remark #15476: scalar cost: 419 remark #15477: vector cost: 91.750 remark #15478: estimated potential speedup: 4.320 remark #15488: --end vector cost summary --

  version of slabAllNuclides on GridCL. If the bank size is set to 1, the program will run in a history-based way. If it is superior to 1, a pseudo event-based way. The tests are carried out with certain thread numbers and a single Nvidia V100 card, as depicted in Figure4.6.

Figure 4 . 6 :

 46 Figure 4.6: Bank size analysis of the CUDA version of slabAllNuclides in offload mode, the higher the better.

Figure 4 . 7 :

 47 Figure 4.7: Scaling performance of the CUDA version of slabAllNuclides in offload mode, the higher the better.

Figure 4 . 8 :

 48 Figure 4.8: slabAllNuclides performance speedup on GridCL, the higher the better -Baseline for the performance is obtained in host mode with the use of pseudo event-based method and a fixed number of bank size, 100.

Figure 4 .

 4 Figure 4.8 provides a straightforward view for comparison of performance speedup among different programming models on GridCL. From the figure we find that in general the pseudo event-based method (PEB) outperforms the history-based method (HB) in offload mode. Comparing to the baseline performance obtained by the pseudo eventbased method fixing the bank size to 100, the CUDA version may reach up to approximately 6.0x speedup and the

Figure 4 .

 4 Figure 4.9 shows the performances of different programming models on Cobalt-hybrid. In comparison with the baseline performance in host mode, the CUDA version may obtain at maximal a factor of 5.5x speedup and the OpenACC version reaches up to a factor of 4.5x speedup with the utilization of two Nvidia P100s. If we only target to a single accelerator, the CUDA version gains around 2.6x performance speedup and the OpenACC version gets a factor of 2.3x speedup. Except for these two implementations, the Kokkos version in offload mode obtains worse performance than its performance in host mode.

Figure 4 .

 4 Figure 4.10 depicts the performance speedups of different programming models on Cobalt-v100. By using 4 Nvidia GPUs, the CUDA and OpenACC versions may respectively reach up to at maximal 10.2x speedup and 5.2x speedup compared with the baseline performance in host mode. At the same time, they may obtain approximately 6x and 3x speedups with the utilization of 2 Nvidia GPUs. When targeting to a single device, the CUDA version gains a factor of 4.1x speedup and the OpenACC version a factor of 1.7x speedup. The OpenMP offload and SYCLversions are unable to be compiled and executed on Cobalt-v100. As for the Kokkos version, its performance in offload mode is similar to that on GridCL and Cobalt-hybrid, which is much worse than the performances of the CUDA and OpenACC versions.

Figure 4 . 11 :

 411 Figure 4.11: slabAllNuclides performance speedup on Gorgon, the higher the better -Baseline for the performance is obtained in host mode with the use of pseudo event-based method and a fixed number of bank size, 100.

Figure 4 . 13 :

 413 Figure 4.13: Performance speedup of the CUDA version of slabAllNuclides adopting different optimization strategies tested on GridCL with a single Nvidia V100 and 20 CPU threads, the higher the better.

Figure 5 . 1 :

 51 Figure 5.1: An example of 3P principle of the implemented programming models including OpenMP thread, OpenMP offload, OpenACC, Kokkos, SYCL.

  of code, we separate each line into a group of operators and operands and then obtain the program length N by summarizing N 1 and N 2 where N 1 and N 2 are the total number of occurrences of operators and operands. Similarly, the vocabulary n is aggregated by n 1 and n 2 where they represent respectively the total number of distinct operators and operands. Based on N and n, the program volume V and program difficulty D are defined as V = N log 2 n and D = n1 2 • N2 n2 . The program effort E is calculated by V and D with the formula E = V • D. Parameters mentioned above mainly including N , V and E can be used as metrics to evaluate implicit productivity. Moreover, in order to address the program length of a platform-specific code path for an application, a new metric PS-N (platform-specific program length) is proposed the same as PS-LOC and PS-NW.

  a, p, h) and chn i,j (a, p, h) program volume/effort N (program length) and PS-N LOC and PS-LOC, NW and PS-NW relative metrics derived from items above

Property 5 . 1 . 1 andProperty 5 . 1 . 3 FProperty 5 . 1 . 4 F

 511513514 Property 5.1.2, we have a few properties listed as follows: gm is the most adaptive metric applying to a data set C(a, p, h) retrieved by a metric indicating the absolute performance since in this scenario the percentage of change for F gm is platform-independent. am is the most adaptive metric applying to a data set C(a, p, h) retrieved by a metric indicating the relative performance since in this scenario the percentage of change for F am is platform-independent.

Similarly, we explainT 1 F

 1 in detail the results retrieved by the variants of F am . With reference to the metric T 1 F gm • C(a, p, h), the results are equal to those obtained by the original metric F am , which corresponds to the equation am = F am . With reference to the metric F gm • (T 2 C(a, p, h)), the results shows that the OpenMP thread and OpenACC versions gain the highest value, 59% while the OpenMP offload version obtains the lowest one, 9%. The priorities of different programming models for better application efficiency become OpenMP thread = OpenACC > CUDA > Kokkos > SYCL > OpenMP offload. With reference to the metric T 3 F • C(a, p, h), the priority order resembles in the one obtained in Table5.5, which is OpenACC > Kokkos > OpenMP thread > CUDA > OpenMP offload > SYCL.

Table 5 .

 5 7 shows a variety of metrics of implicit productivity for implementations of slabAllNuclides across a set of platforms. The measured source files include the partial offloading part to transfer nuclide data to device and calculate microscopic cross sections on device. With respect to LOC, the results are measured by using the command line wc -l. The CUDA version requires the largest number of LOC whereas the OpenMP thread version the smallest. The PS-LOC of each implementation for Intel x86-based, IBM Power-based, Intel GPU-based and Nvidia GPU-based platforms are same and thus the corresponding code distances d i,j between two of four types of platforms for all implementations are equal to zero. It indicates that the OpenACC, OpenMP offload, Kokkos and SYCL versions are of high abstraction which leads to little difference between distinct platform-specific code paths. The churns of different code paths for each implementation reflect that the OpenMP offload version share a single code path for both Intel x86-based, IBM Power-based and Nvidia GPU-based platforms and the largest difference between code paths is introduced by the Kokkos version

Title:

  Evaluation of programming models for manycore and / or heterogeneous architectures for Monte Carlo neutron transport codes Keywords: Manycore architectures, Heterogeneous architectures, Particles transport Abstract:

2.4.1.1 Message Passing Interface

	Message Passing Interface (MPI) is a message-passing library interface specification designed for functioning on
	a wide variety of parallel computing architectures, allowing convenient C, C++, Fortran bindings [47]. It addresses
	principally the message-passing parallel programming model, in which multiple processes communicate by calling
	library routines to send and receive messages to other processes.
	MPI is a standard, not an implementation. There are many implementations of MPI, such as OpenMPI (an
	open source MPI-2 implementation) [51] or MPICH2 (a freely available, portable implementation of MPI developed
	originally by Argonne National Laboratory) [77]. In most MPI implementations, a fixed set of processes is created at

initialization of program. Moreover, these processes may execute different programs, which is a big distinction from SPMD model where every process executes the same program. MPI has been designed by a community of parallel computing vendors, computer scientists, and application developers. Nowadays MPI is still the most widely used model in HPC domain.

Table 2 .

 2 1 shows the comparison of several programming models implemented in our thesis work, including CUDA,

	Kokkos, OpenMP (thread + offload), OpenACC, and SYCL. From the perspective of thread hierarchy, all program-
	ming models provide a number of abstract levels of parallelism, enabling themselves to map to hardware with
	appropriate features except for CUDA, which only provides three levels of parallelism targeting specifically to the
	Nvidia GPUs.

Table 2

 2

		.1: Comparison of implemented programming models.
	Programming Models Thread Hierarchy	Memory Hierarchy	Targeting Architectures
	CUDA	block, warp, thread		Nvidia GPUs
	Kokkos	team, thread, vector	Programmable shared memory	CPUs, Nvidia GPUs
	SYCL	work-group, work-item		CPUs, Intel GPUs
	OpenMP OpenACC	team, thread, simd gang, worker, vector	Lack of interface to on-chip memory	CPUs, Nvidia GPUs CPUs, Nvidia GPUs
	From the perspective of memory hierarchy, CUDA, Kokkos, and SYCL provide an interface to explicitly man-
	age shared memory while OpenMP, OpenACC both manage on-chip memory implicitly with the utilization of com-
	piler directives or high-level functions, giving users no interface to programme shared memory explicitly. Note
	that CUDA uses __shared__ to declare variables on shared memory. Kokkos makes use of the member function
	team_shmem().get_shmem() to explicitly get an allocation on the scratch pad memory. SYCL provides access to
	allocated shared memory via local memory if a SYCL accessor has the access target target::local.
	From the perspective of targeting architectures, CUDA is dedicated to targeting to Nvidia GPUs. Kokkos sup-
	ports multicore/manycore processors (x86, OpenPower, ARM, etc) and Nvidia GPUs. SYCL is capable of targeting
	to any CPU (x86, OpenPower, ARM), Nvidia/AMD/Intel GPUs as well as Intel/Xilinx FPGAs based on different
	implementations. OpenMP and OpenACC may port codes to multicore CPUs and Nvidia GPUs.

  In this case, an energy lookup with input energy E is indeed to find the lower-upper bound E k and E k+1 with the relation E k 6 E 6 E k+1 . The corresponding microscopic cross section is then retrieved by one of linear-linear interpolations as described in Equation 3.8.

Input: a group of M couples of microscopic cross sections, atomic densities, { i ,N i } i2M , corresponding to M nuclides which compose a given material Result: calculated macroscopic cross section Σ of the material 1 Σ=0; 2 foreach composing nuclide of the material, indexed with i, i 2 M do 3 Σ+=N i i ; 4 end

3.4.4.2 Pre-tabulated Microscopic Cross Section Calculation Algorithm

Pre-tabulated approach mainly consists two steps, an energy lookup and an interpolation. Typically, the continuousenergy cross section data is stored as a long table of cross sections (cross section grid) which is bijective to another long table of energy points (energy grid).

  H, 3 H) whereas other nuclides have more than a million energy points (e.g. 238 U). This feature leads to considerable variations of workloads between energy lookup processes through different nuclides and makes efficiency of energy lookup algorithms nuclide-dependent via different search schemes.

Table 4 .

 4 1: Typical Monte Carlo neutron transport run time percentage in PATMOS.

	Processing Step

  shown in Algorithm 9 (at line 19). Besides, this loop can be vectorized by auto-vectorization of Intel compiler since version 16. It should be noted that the OpenMP thread version is the reference implementation on CPU.

  .2, we find that OpenMP thread can be compiled by most of compilers except NVCC where NVCC installed. hipSYCL[START_REF] Alpay | SYCL beyond OpenCL: The Architecture, Current State and Future Direction of hipSYCL[END_REF] is a LLVM-based compiler building on CUDA, HIP and OpenMP backends instead of OpenCL which allows it to target multicore CPUs, Nvidia GPUs and AMD GPUs. It is wrapped with a clang plugin to provide support for SYCL code and thus interoperable with CUDA and HIP codebases. ComputeCpp is a LLVM-

	based device compiler capable of compiling SYCL implementations on OpenCL platforms which mainly include
	Intel CPUs, Intel GPUs, Intel FPGAs, AMD GPUs with SPIR(-V) support and Nvidia GPUs with experimental PTX
	support

is Nvidia's CUDA Compiler based on LLVM compiler infrastructure which is dedicated to compiling CUDA code.

Concerning OpenACC, PGI starts to support OpenACC 1.0 specification from its 2012 release and the full features of directives-based parallel programming specified in OpenACC 2.6 as well as some features of OpenACC 10.2

Table 4 .

 4

	Machines

3: List of machines to run tests of slabAllNuclides benchmark.

Table 4 .

 4 4: Configurations of processors composing accessible machines.

	Processors	Cores Threads Frequency GHz Memory Bandwidth GB/s FP64 TFLOPS
				Skylake
	Intel Xeon Gold 6138	40	80	2.00	207.4	2.6
	Intel Xeon Gold 6148	40	80	2.40	207.4	3.1
				Broadwell
	Intel Xeon E5-2680 v4	28	56	2.40	76.8	1.1
				Kaby Lake G
	Intel Core i7-8705G	4	8	3.10	37.5	0.2
				Power9
	Power9 SMT4	40	160	3.20	340	1.0
				Nvidia Pascal
	P100	56	1792	1.33	720	5.3
				Nvidia Volta
	V100	80	2560	1.30	900	7.8
				Intel Gen9.5
	Intel HD Graphics 630	24	168	0.35	37.5	0.03

Table 4 .

 4 5: Performance of the OpenMP thread version of slabAllNuclides using the history-based method via different thread affinity policy.

	Number of threads	cores/close	cores/spread
	10	1269.89	1113.39
	20	732.26	639.34
	30	472.89	462.95
	40	374.05	373.86
	50	352.12	351.85
	60	331.73	329.19
	70	310.20	307.66
	80	294.50	291.38

Table 4 .

 4 6: Comparison of the performances of slabAllNuclides in offload mode adopting different thread affinity policies tested on GridCL with a single Nvidia V100.

	Number of threads	HB	cores/close	PEB	HB	cores/spread	PEB
	2	1168.85	393.20	1278.35	399.46
	4	725.82	215.61	955.73	207.73
	6	1383.32	153.42	1794.62	149.67
	8	1731.11	124.94	2004.87	122.95
	10	1845.44	112.84	1592.36	107.69
	20	1905.26	96.90	1868.24	96.87
	40	2475.58	96.69	2556.41	96.83
	60	3047.89	97.16	4006.39	96.96
	80	3980.08	141.11	4062.55	143.19

Table 4 . 7 :

 47 Comparison of the performances of slabAllNuclides in offload mode adopting different thread affinity policies tested on GridCL with two Nvidia V100.

	Number of threads	HB	cores/close	PEB	HB	cores/spread	PEB
	2	1087.49	388.09	1071.43	377.84
	4	587.34	202.77	632.48	195.39
	6	442.70	141.48	483.95	137.02
	8	397.23	110.35	864.31	107.81
	10	594.39	99.27	969.56	90.39
	20	621.76	62.15	933.36	57.77
	40	927.63	49.20	943.89	49.15
	60	1669.81	49.39	1699.01	49.85
	80	1971.01	49.83	1966.18	49.99

Table 4 .

 4 8.

Table 4 .

 4 8: Particle tracking rate via different programming models on a set of platforms, the higher the better. OMPth refers to OpenMP thread, OMPth+offload means OpenMP host and offload functionalities.We find that among the implemented programming models, the CUDA and OpenACC versions allow to obtain better performance. For example, at maximal 101.7 ⇥ 10 2 and 59.8 ⇥ 10 2 particles/s on GridCL, 147.0 ⇥ 10 2 and 80.6 ⇥ 10 2 particles/s on Gorgon). The OpenMP offload version is not competitive with the OpenACC version, for instance, merely leading to a 5 10% tracking rate of the CUDA version on GridCL and Gorgon. The Kokkos version obtains even much worse performance than the OpenMP offload version, with the peak performances of 0.7 ⇥ 10 2 particles/s on GridCL and 1.2 ⇥ 10 2 particles/s on Gorgon. As for the SYCL version, it is the single version capable of offloading computational kernel to Intel Graphics while there is a lack of compiler support to successfully compile and execute the SYCL version on other machines equipped with Nvidia GPUs. The peak performance of the SYCL version in offload mode is 0.6 ⇥ 10 2 particles/s on Intel NUC.

				slabAllNuclides
	Machine		Programming Model	(⇥10 2 particles/s)
				HB	PEB
	GridCL	CPU (40 cores, HT2)	OMPth	17.2	16.9
			OMPth+ACC	7.7	6.9
			OMPth+offload	6.5	5.5
			OMPth+Kokkos	12.8	8.1
		1V100	OMPth+CUDA	5.2	51.8
			OMPth+ACC	3.5	34.7
			OMPth+offload	0.4	2.6
			OMPth+Kokkos	0.4	0.7
		2V100	OMPth+CUDA	10.5	101.7
			OMPth+ACC	7.0	59.8
			OMPth+offload	0.6	4.5
	Cobalt-hybrid	CPU (28 cores, HT2)	OMPth	10.2	9.5
			OMPth+ACC	4.5	4.4
			OMPth+Kokkos	7.1	5.9

where

Table 4 .

 4 9: Profiling of the CUDA version of slabAllNuclides in offload mode on GridCL using nvprof.

	GPU activities	Time Percentage(%)	Execution Time (s)
		HB	PEB	HB	PEB
	sigma1DopplerBroadening	85.81	72.41	9.17	1.62
	CUDA memcpy HtoD	9.36	22.33	1.00	0.50
	CUDA memcpy DtoH	4.84	5.27	0.52	0.12

Table 4 .

 4 10: Metric values of the computational kernel sigma1DopplerBroadening() of the slabAllNuclides CUDA version using nvvp.

	Metrics

Table 4 .

 4 11: Metric values of the computational kernel sigma1DopplerBroadening() of the OpenACC, OpenMP offload, Kokkos versions of slabAllNuclides using the PEB method.

	Metrics	OpenACC	OpenMP offload	Kokkos
	Grid size	[25996,1,1]	[25996,1,1]	[25996,1,1]
	Block size	[32,1,1]	[32,1,1]	[1,32,1]
	Registers/Thread	110	200	70
	Shared Memory/Block	276B	0B	0B
	Theoretical Block Limit	16	10	28
	Global Memory Load Throughput	195.51GB/s	57.31GB/s	62.24GB/s
	Global Memory Store Throughput	2.19GB/s	0.63GB/s	0.31GB/s
	Global Memory Load Efficiency	64.94%	63.82%	25.00%
	Global Memory Store Efficiency	25.00%	25.00%	25.00%
	Achieved Occupancy	18.02%	12.59%	20.16%
	Double-precision Floating-point Efficiency	14.13%	5.61%	5.22%
	Branch Efficiency	93.16%	97.50%	100.00%
	Warp Execution Efficiency	39.96%	27.00%	100.00%

version. The memory bandwidth of the OpenACC version is also worse than the CUDA version. Furthermore, the OpenACC version makes use of shared memory automatically and each block contains 276 bytes shared memory.

Table 4 .

 4 12: Comparison of the performances of the Kokkos versions of slabAllNuclides tested on GridCL with a single Nvidia V100.

	Particle Tracking Method	Simulation Time (s/cycle) Kokkos basic Kokkos optimized
	HB	426.95	393.26
	PEB	97.63	87.60

Table 4 .

 4 [START_REF]OpenACC Application Programming Interface[END_REF]: Granularity analysis of the CUDA version of slabAllNuclides performed on GridCL with a single Nvidia V100 and 20 CPU threads.

	NUCS_PER_WARP	Simulation Time (s/cycle) HB PEB
	1	2257.97	93.12
	2	2086.65	74.95
	4	1927.87	73.39
	8	2228.79	85.98
	16	2740.79	109.98
	32	3080.32	159.34

Table 4 .

 4 14: Comparison of the performances of the CUDA implementations of slabAllNuclides tested on GridCL with a single Nvidia V100 and 20 CPU threads.

	Optimizations	Simulation Time (s/cycle)
		HB	PEB
	warp	1868.24	96.87
	zero-copy	933.20	198.48
	maxrregcount-64	1600.47	91.69
	maxrregcount-32	1664.74	222.70
	unified-memory	4775.12	237.65
	multi-kernels	2121.47	80.99
	4-nucs-per-warp	1927.87	73.39

Definition 5.1.2

  An application is a program which adopts specific algorithms and programming models to handle one or more problems and produce an output.

	c = C(a, p, h)	(5.2)
	Definition 5.1.3 A problem is actually a suite of input parameters that describes a real problem.
	Note that in our case, slabAllNuclides fed with different number of threads, number of particles, number of volumes,
	bank size are considered as different problems.
	Definition 5.1.4 A platform is a hardware and software infrastructure that supports the execution of an application
	(architecture, operating system, compiler, etc.).
	Definition 5.1.5 A characteristic is a value describing a property such as performance, productivity, portability,
	maintainability, energy consumption and so on. This value is obtained by other metrics with the following formula:

  Since our objective is to establish a generic model which applies to all three dimensions, the mathematical expressions should be independent of each dimension. Thus we define F as a simplified notation of {F A , F P , F H } Note that the distribution of all elements in the data set {C i (a, p, h)} i=1,2,...,N is better one-dimensional because it is quite meaningless to aggregate values diversifying on multiple dimensions (e.g. summarization of three values C(a 2 ,p 1 ,h 1 ), C(a 1 ,p 2 ,h 1 ) and C(a 1 ,p 1 ,h 2 ) does not help us to identify the exact influence of each dimension on the property). Therefore, we exclude possibilities of complicate interactions between different dimensions in our generic model.There are a variety of mathematical models which are extensively used to perform statistical computations on data sets. In what follows we define our generic operator F by Pythagorean means (arithmetic, geometric, harmonic)[99][44] which are respectively denoted as F am , F gm and F hm . The metric F am is an operator computing the arithmetic mean of a data set of N values {C i (a, p, h)} i=1,2,...,N , defined by Equation

	F • C(a, p, h)=F •{C i (a, p, h)} i=1,2,...,N =	8 > > > < > > > :	C(a, p, {h i })=C(a, p, H) if 8h i 2 H C({a i },p,h)=C(A, p, h) if 8a i 2 A C(a, {p i },h)=C(a, P, h) if 8p i 2 P	(5.4)
	5.1.1.1 Arithmetic Mean
	Definition 5.1.6
			A • C(a, p, h)
	P(c, a, P, h)=F P • C(a, p, h)	(5.3)
	H(c, a, p, H)=F H • C(a, p, h)
	and we obtain:

  weighted versions of F am and F hm applying to C(a, p, h) with C(a, p, h)={

	resemble in Equation	p 1 t 1	,	p 2 t 2	,...,	p N t N	} and the formulas

  ti and we can always find a weighted version of F am or F hm to obtain it. However, for the standard versions of F am or F hm , it is not safe to directly apply them to any

	, we know that for a data set C(a, p, h)={ single-platform property, its optimal estimate is equaling to P N i pi P N i data set in a form like C(a, p, h)={ p 1 t 1 , p 2 t 2 ,..., p N t N }. We better fix one quantity and measure the other. p 1 t 1 , p 2 t 2 ,..., p N } characterizing an intrinsic t N

  Definition 5.1.12 Portability refers to the ability of an application a to solve a problem p correctly on a given set of platforms H.Thus, the generic metric of portability is F applying to a specific data set C(a, p, h) and we have three different sub-types, F am , F gm and F hm . Since the portability of an application solving a problem on a platform can be described by a boolean value, 0 or 1, referring to non-portable or portable, the values of elements {C i (a, p, h)} i=1,2,...,N

	can be expressed as:
	8i =1, 2,...,N, C i (a, p, h)=	8 < :	0 if (a, p) fails to run correctly on h 1 if (a, p) runs correctly on h	(5.19)

Table 5 .

 5 1: Execution time of applications {a 1 ,a 2 ,a 3 } solving a given problem across a set of platforms {h 1 ,h 2 }, the lower the better.

	Applications	Platforms	Performance Improvement	Expected F • C(a, p, h)
		h 1	h 2
	a 1	10	100	-	-
	a 2	7.5	100	25%	equal
	a 3	10	75	25%	equal

Table 5 .

 5 2: Portability of applications {a 1 ,a 2 ,a 3 ,a 4 } solving a given problem across a set of platforms {h 1 ,h 2 ,h 3 ,h 4 ,h 5 ,h 6 }.

	Applications			Platforms
		h 1	h 2	h 3	h 4	h 5	h 6
	a 1	1		111
	a 2		1	111
	a 3	111		1
	a 4						1

As an alternative, we propose a new type of performance portability which is called as relaxed performance portability. The definition is given as follow: Definition 5.1.

[START_REF] Binder | Applications of the Monte Carlo Methods in Statistical Physics[END_REF]

Relaxed performance portability (RPP) limits the boolean feature of portability and addresses more about performance, which means that for an application a solving a problem p, the estimated value will be zero only if this application is blocked by all platforms in the platform set H.

  4. Not to change the nature of F am and F gm as described in Property 5.1.1 and Property 5.1.2.Let S and N the sub-sets of a platform set H with S [ N = H. All platforms in S support an application a solving a problem p while all platforms in N don't support the couple (a, p). According to the criterion, we may firstly calculate F • C(a, p, h) h2S , denoted as F 0 , the estimated value of F applying to a data set C(a, p, h) where each h is in the platform sub-set S and then make use of this value in variants of F .

Table 5 .

 5 4: Metrics of performance of different implementations of slabAllNuclides across a set of platforms.

	Machines	#GPUs	Programming Models OpenMP thread CUDA OpenACC OpenMP offload Kokkos SYCL
			Particle Tracking Rate (⇥10 2 particles/s)
	GridCL	0	16.9	6.9	5.5	8.1
		1	51.8	34.7	2.6	0.7
	Cobalt-hybrid	0	9.5	4.4		5.9
		1	25.7	21.8		2.4
	Cobalt-v100	0	13.6	6.4		7.2
		1	55.7	21.8		2.8
	Gorgon	0	6.8	7.4	1.5	6.0
		1	55.0	36.6	6.5	1.2
	Intel NUC	0	3.1	1.1		2.2	0.1
		1					0.6
			Application Efficiency (%)
	GridCL	0	100	41	33	48
		1	100	67	5	1
	Cobalt-hybrid	0	100	46		62
		1	100	85		9
	Cobalt-v100	0	100	47		53
		1	100	39		5
	Gorgon	0	92	100	20	81
		1	100	67	12	2
	Intel NUC	0	100	35		71	3
		1					100

3 ,

 3 F gm is the most adaptive metric applying to a data set C(a, p, h) of particle tracking rate since in this case the percentage of change for F gm is platform-independent to meet the requirement of adaptability of performance portability expressed in Definition 5.1.17. Therefore, we mainly give a detailed analysis of results obtained by F gm which are highlighted with green color.With respect to the metric T 1 F gm • C(a, p, h), the highest value is obtained by the CUDA version, 18.0 ⇥ 10 2

particles/s, while the lowest value is retrieved by the SYCL version, 0.0 particles/s by keeping one decimal place. The priorities of different programming models for better performance become CUDA > OpenACC > OpenMP thread > Kokkos > OpenMP offload > SYCL.

Table 5 .

 5 5: Evaluation of performance portability for different implementations of slabAllNuclides across a set of platforms via the generic model F and the absolute performance metric, particle tracking rate, the higher the better. , T 2 and T 3 are three transfer operators defined in section 5.1.4.3 to obtain variants of F . Note that the value of T 3 F • C(a, p, h) is a score which has no unit.With respect to the metric F gm • (T 2 C(a, p, h)), the priority order is exactly the same as the metricT 1 F gm • C(a, p, h).The CUDA version gains the best result which surpasses the OpenACC version with a factor of 1.5x speedup. With respect to the metric T 3 F • C(a, p, h), the highest score is retrieved by the OpenACC version while the lowest score is obtained by the SYCL version. The OpenACC and Kokkos versions benefit from the largest integer part which represents the number of supported platforms for a given programming model. As for the comparison between the OpenACC and Kokkos version, the decimal part of the OpenACC version is larger than that of the Kokkos version, which means that the normalized aggregating result of the OpenACC version surpasses the one of the Kokkos version. Overall, this variant of F gm renders a priority order as OpenACC > Kokkos > OpenMP thread > CUDA > OpenMP offload > SYCL. Now we consider the evaluation of performance portability for different implementations of slabAllNuclides via the generic model F and the relative performance metric, application efficiency, as illustrated in Table 5.6. Similarly, the values highlighted with green color are gained by the utilization of F am which is the most adaptive metric for the evaluation of application efficiency according to Property 5.1.4. Since in this scenario the percentage of change , T 2 and T 3 are three transfer operators defined in section 5.1.4.3 to obtain variants of F . Note that the value of T 3 F • C(a, p, h) is a score which has no unit.

	Metrics		Particle Tracking Rate (⇥10 2 particles/s) OpenMP thread CUDA OpenACC OpenMP offload Kokkos SYCL
		F am	5.0	18.8	14.1	1.6	3.7	0.1
	F • C(a, p, h)	F gm	0.0	0.0	0.0	0.0	0.0	0.0
		F hm	0.0	0.0	0.0	0.0	0.0	0.0
	F 0 • C(a, p, h)	F 0 am F 0 gm F 0 hm	10.0 8.6 7.1	47.1 44.9 42.4	15.7 9.9 5.2	4.0 3.4 2.9	4.1 3.1 2.2	0.4 0.2 0.2
		F am	5.0	18.8	14.1	1.6	3.7	0.1
	T 1 F • C(a, p, h)	F gm	4.3	18.0	9.0	1.4	2.8	0.0
		F hm	3.5	17.0	4.7	1.2	2.0	0.0
		F am	6.0	23.5	15.7	2.0	4.1	0.1
	F • (T 2 C(a, p, h))	F gm	3.8	15.3	9.9	1.2	3.1	0.0
		F hm	2.4	10.6	5.2	0.7	2.2	0.0
		F am	5.21	4.99	9.33	4.09	9.09	2.01
	T 3 F • C(a, p, h)	F gm	5.19	4.99	9.22	4.08	9.07	2.01
		F hm	5.17	4.99	9.12	4.07	9.05	2.00
	where T 1

Table 5 .

 5 9: Metrics of explicit productivity for implementations of slabAllNuclides across a set of platforms.

	Machines	#GPUs	Programming Models OpenMP thread CUDA OpenACC OpenMP offload Kokkos SYCL
				Relative Speedup (x)
	GridCL		2.45		1.00	0.80	1.17
				1.49	1.00	0.07	0.02
	Cobalt-hybrid		2.16		1.00		1.34
				1.18	1.00		0.11
	Cobalt-v100		2.13		1.00		1.13
				2.56	1.00		0.13
	Gorgon		0.92		1.00	0.20	0.81
				1.50	1.00	0.18	0.03
	Intel NUC		2.82		1.00		2.00	0.09
				Relative PS-LOC (x)
	GridCL		0.80		1.00	0.98	0.98	1.17
				2.14	1.00	0.98	0.98
	Cobalt-hybrid		0.80		1.00	0.98	0.98	1.17
				2.14	1.00	0.98	0.98
	Cobalt-v100		0.80		1.00	0.98	0.98	1.17
				2.14	1.00	0.98	0.98
	Gorgon		0.80		1.00	0.98	0.98
				2.14	1.00	0.98	0.98
	Intel NUC		0.80		1.00	0.98	0.98	1.17
	GridCL		3.06	Ψ relative	1.00	0.82	1.19
				0.70	1.00	0.07	0.02
	Cobalt-hybrid		2.70		1.00		1.37
				0.55	1.00		0.11
	Cobalt-v100		2.66		1.00		1.15
				1.20	1.00		0.13
	Gorgon		1.15		1.00	0.20	0.83
				0.70	1.00	0.18	0.03
	Intel NUC		3.52		1.00		2.04	0.08

Table 5 .

 5 10: Evaluation of productive performance portability for different implementations of slabAllNuclides across a subset of platforms via the generic model F , Ψ relative , the higher the better. The variants of F am transferred by T 1 and T 2 obtain the similar results. The priorities of different programming models for better productive performance portability are both OpenMP thread > OpenACC > Kokkos > CUDA > OpenMP offload > SYCL where the highest values gained by the OpenMP thread version of two variants are respectively 1.454 and 1.745 and the lowest values retrieved by the SYCL version of two variants are respectively 0.009 and 0.018. As for the metric T 3 F am • C(a, p, h), the OpenACC version renders the highest score, 5.99, while the SYCL version obtains the lowest score, 1.03. The priority order retrieved by this variant resembles in the one attained by the same variant addressing performance portability, which is OpenACC > Kokkos > OpenMP thread > CUDA > OpenMP offload > SYCL.

	Metrics

  Évaluation de modèles de programmation pour les architectures manycore et / ou hétérogènes pour les codes de transport neutronique Monte Carlo Mots clés: Architectures manycore, Architectures hétérogènes, Transport des particules Résumé: Dans cette thèse nous nous proposons d'évaluer les différents modèles de programmation disponibles pour adresser les architectures de type manycore et / ou hétérogènes dans le cadre des codes de transport Monte Carlo. On considèrera dans un premier temps un cas test d'application simple mais représentatif pour couvrir un éven-tail assez large de solutions et les comparer en terme de performance, de portabilité de la performance, de facilité de mise en oeuvre et de maintenabilité.Les architectures cibles sont les CPU 'classique', Intel Xeon Phi et GPU. Les modèles de programmation les plus pertinents seront ensuite mis en place dans un code de transport Monte Carlo.

Titre:

Carlo développé au CEA, appelé PATMOS. Une version de déchargement hétérogène du code a été développée via une méthode basée sur l'historique et une méthode basée sur des pseudo événements. Plusieurs modèles de programmation (thread OpenMP, OpenMP offload, OpenACC, CUDA, Kokkos et SYCL) ont été implémentés pour l'évaluation et les tests ont été réalisés sur différentes architectures (x86 ou OpenPower + GPU). Un ensemble de métriques a été adopté pour évaluer notre simulation Monte-Carlo de déchargement hétérogène en termes de portabilité, de performances et de productivité. L'analyse des métriques a été réalisée avec la proposition d'une métrique générique et son qui atténue les défauts des métriques d'origine, les rendant plus adaptables pour gérer la comparaison numérique des applications en termes de portabilité des performances et de portabilité de la productivité à travers une variété d'informatique architectures.En un mot, l'objectif de la thèse est d'évaluer les métriques existantes en termes de portabilité, de performance et de productivité et de proposer des améliorations de ces métriques sur un exemple d'application réelle réalisant une simulation de transport de neutrons de Monte-Carlo. Les principaux objectifs sont les suivants:1. proposer un nouvel algorithme de suivi des particules, la méthode pseudo event-based.2. analyser et régler les performances de la version CUDA.3. mettre en oeuvre des codes de transport de neutrons Monte-Carlo portables via des modèles de programmation de haut niveau.4. proposer et analyser une métrique générique et ses variantes.5. évaluer les implémentations portables des codes de transport Monte-Carlo par la métrique générique.

Acknowledgements

First and foremost, I would like to express my deep and sincere gratitude to my thesis supervisor, Christophe CALVIN, for the continuous support of my Ph.D study throughout the past three years, for his motivation, enthusiasm and vision which have deeply inspired me every time I got stuck on research. I am extremely grateful for having the opportunity to study and work under his guidance. I would like to thank my advisor, Emeric BRUN, for his patience, insightful comments and hardware support. I am grateful for his aids on helping me carrying out tests on the Cobalt-hybrid, Cobalt-v100 platforms. Moreover, I thank him and Philippe Thierry for providing me an Intel NUC machine to accomplish my research work. I thank the rest of my thesis committee: Michael HEROUX and Jean-François MEHAUT, for being my reviewers and writing me the literature of review and Raymond NAMYST, Marc VERDERI, Emmanuelle SAILLARD, Francieli ZANON BOITO, for being the examiners at my thesis defense. My sincere thanks also goes to Fausto MALVAGI, the director of my laboratory SERMA/LTSD in CEA Saclay and Edouard Audit, the director of Maison de la Simulation, for allowing me studying and working in two places throughout the research. I would also like to express my gratitude to my fellow labmates in these two laboratories, for their administrative, academical and hardware supports as well as insightful talks with me. I would like to thank David Chamont for offering me access to the GridCL platform and being my academical advisor. I thank Pierre Kestener for giving his useful comments on the implementation of Kokkos in the Monte Carlo neutron transport code. I thank Serge G. Petiton, for guiding my master internship in Maison de la Simulation and enlightening me the first glance of research in the field of HPC.

I am very much thankful to my parents and twin brother for their love and sacrifices for supporting me continuously to chase my dream.

The clause present(...) ensures that the data required by the function sigma1DopplerBroadening() is already copied to the device. The clause async(stream) where stream may be a nonnegative scalar integer expression makes the host multi-threads possible to process kernel computations or data operations in parallel, which equals to the CUDA asynchronous multi-stream mechanism. Besides, the reduction operation can be achieved by using 3.0 specification such as deep copy directives have been implemented in recent releases. Besides, GCC provides implementation of OpenACC specification from GCC 5 and its latest release GCC 10 starts to support some features of OpenACC 2.6 specification.

With respect to OpenMP offload, GCC offers support for offload of OpenMP 4. [START_REF] Reyes | Codeplay Contribution to DPC++ brings SYCL Support for NVIDIA GPUs[END_REF]. This implementation should work on Nvidia GPUs of capability 5.0 or higher with CUDA compiled by Intel compiler with vectorization support enabled. This is the reason that their performances outperform the performance of the OpenACC version which is compiled by PGI compiler.

Overall, based on the peak performances expressed above, we end up with several conclusions listed as follows:

1. Pseudo event-based method (PEB) surpasses significantly history-based method (HB) with a magnitude of performance speedup in offload mode while its performance is a little less than the performance of historybased method in host mode.

2. In host mode, the implementation that is able to be compiled by Intel compiler with the vectorization options enabled such as the OpenMP thread and Kokkos versions gain better performance than those implementations requiring to be compiled by other compilers such as the OpenACC, OpenMP offload versions.

3. In offload mode (Nvidia P100, V100, Intel Gen9.5):

• The CUDA version renders the highest performance via PEB comparing to other programming models.

• The OpenACC version generally gains approximately 50%-60% of the CUDA performance via both HB and PEB.

• The OpenMP offload version provokes a huge performance degradation via either HB or PEB which is caused by underdeveloped support of CUDA multi-streams for our OpenMP offload implementation achieved by LLVM-Clang and IBM XLC compilers comparing to the CUDA and OpenACC versions.

• Our Kokkos version is limited to offload computational kernels to a single device and does not support CUDA multistreaming. Its performance in offload mode is less competitive to the baseline performance in host mode which is completely unsatisfactory.

• Our SYCL version can only be compiled on Intel NUC associated with Intel Graphics. However, its performance in offload mode is much worse than the baseline performance with a factor of 0.1-0.2x speedup. More optimizations and tests are required with the use of more architectures and the next version of SYCL compilers since SYCL is a preliminary programming model and its implementations achieved by different compilers are immature.

Profiling and Optimization

In this section, we intend to perform a series of profiling to understand better about the conclusions drew in the previous section. The analysis is carried out on GridCL focusing on the offload mode by using the profiling tools such as nvprof and nvvp to ensure that most of implementations of different programming models can be analyzed.

Unfortunately, the SYCL version is not able to be included in the following analysis. We start the analysis from the comparison between the history-based method and the pseudo event-based method.

refers to the number of nuclides assigned to a warp and sigma_tmp is a table of size 32 allocated in shared memory for data storage.

Original sum reduction

for (int i = 16; i > 0; i »= 1) sigma += __shfl_xor_sync(0xffffffff, sigma, i); } Furthermore, since the precalculation part includes two invocations of the function getLowerBound(), we may also execute these two functions concurrently with the utilization of shared memory. The idea is to divide all 32

Adaptability to Single-platform Properties

Since a single-platform property reflects an intrinsic characteristic of the platform, in order to estimate it correctly, it is important to summarize the measured values for each quantity. Hence, we define the adaptability of F as follow:

Alternatively, if C(a, p, h) is denoted as the inverse of time, the standard version of harmonic mean metric F hm will be adaptive towards the data set C(a, p, h)={ 1 t i } i=1,2,...,N because the final result is inversely proportional to the sum of measured time, as expressed below:

As for the third case, we consider C(a, p, h) a ratio such as floating point operations per second (GLOPS) with

Since this estimate is proportional to P N i=1 p i , we can say that F am is an accurate measurement of floating point operations per second when all values in C(a, p, h) are measured with the same execution time. For example, an application counting the total single-precision floating point operations within the same time step. However, if we adopt another strategy of measurement by which each C(a, p, h) is evaluated with fixed number of floating point operations, denoted as P , then we will have:

The result is proportional to the sum of the inverses of t i instead of the inverse of the sum of t i . Therefore, discussed in previous section. From the other side, how to choose a metric and properly use it highly depends on the objective of users and the property to be aggregated. Therefore, we need to know more about the nature of three metrics, F am , F gm , F hm .

In general, F am , F gm and F hm have a relation described below, the equality among them is achieved under the condition that all elements in the data set C(a, p, h) are equal.

ci and we have:

)

(5.17

In addition, F am and F gm accept both positive numbers and zero while F hm only applies to positive number.

Note that F gm =0if any element in C(a, p, h) equals to zero. However, since F gm usually cooperates with logarithm metric with C i (a, p, h)=e i (a, p), and we have:

With reference to discussions above, the metric P can be expanded from adopting the metrics of relative performance to other metrics of absolute performance with the utilization of multiple aggregating methods. We arrive at a new definition of performance portability along with a new metric of performance portability.

Definition 5.1.16 Performance portability is a measurement adopting various aggregating methods which estimates the performance of an application a running correctly on a platform set H to solve a given problem p with the intention of comparison.

In this way, F becomes the generic metric of performance portability when it applies to a data set C(a, p, h) retrieved by a metric of performance.

Adaptability of F towards performance metrics

Now that we have several metrics of performance portability and a set of performance metrics including execution time, application efficiency and so on, a study of adaptability of F towards these performance metrics is going to be of great significance.

Because the objective of summarization of values is no longer to reflect an intrinsic property, but to simply aggregate measured data for comparison, the common rule that we should follow is that each platform in the platform set should be treated independently and equally. In other words, since each platform has its own scale of performance for an application to solve a problem, the performance improvement with reference to the maximal performance on platforms h 1 and h 2 should make the same contribution to the estimated value of F . Therefore, we have the following definition: Definition 5.1.17 Adaptability of a performance portability metric F towards a data set C(a, p, h) retrieved by a performance metric is maintained only if the contribution made by any platform in a platform set H is treated independently and equally.

For instance, consider an application set A = {a 1 ,a 2 ,a 3 } solving a given problem p across a platform set H = {h 1 ,h 2 } as illustrated in Table 5.1. We assume that comparing to the application a 1 , a 2 and a 3 respectively gains a 25% performance improvement on platforms h 1 and h 2 . According to Definition 5.1.17, they should make same contribution to the estimated value of F , thus F • C(a 2 ,p,h)=F • C(a 3 ,p,h).

Based on Definition 5.1.17, we need to find out in which way the contribution of performance on platforms is expressed by different performance metric. In general, there are two ways which are individually described as a percentage of change and an amount of change.

Metrics of Explicit Productivity Portability

When F applies to explicit productivity metrics, it actually reflects the productive performance portability which represents the average performance improvement per program utility unit across a platform set. In order to obtain meaningful relative speedup or effort, it is required that we have to at least find an application supported by all platforms in the platform set and treat it as the reference implementation. Otherwise we may use different applications as reference on distinct platforms and the aggregation of results will be quite meaningless. Based on this limitation, the variants of F addressing relaxed performance portability cannot be directly used to capture productive performance portability. We have to find the application supported by a maximal number of platforms in the platform set and perform the evaluation of productive performance portability simply across the new sub-set of platforms.

The adaptability of F towards metrics of explicit productivity portability (EPP), also called as productive performance portability, is similar to that of performance portability. We should maintain the feature of platformindependent percentage of change. Since our productive performance portability metrics adopt speedup as performance metric, according to Property 5.1.1, F am is the most adaptive metric to aggregate productive performance portability of an application solving a problem across a set of platforms.

Evaluation of Benchmark

The benchmark slabAllNuclides is used to perform further evaluations. SIGMA1 on-the-fly Doppler broadening approach and pseudo event-based method are chosen along with the input parameters fixed to 5 ⇥ 10 5 particles and 100 bank size. Ten platforms are composed from five machines introduced in the previous chapter in the consideration of single CPU or CPU + 1GPU.

Portability of six versions of slabAllNuclides implemented via programming models of OpenMP thread, CUDA, OpenACC, OpenMP offload, Kokkos as well as SYCL across ten platforms is shown in Table 5.3.

Note that green cells refer to applications portable to the corresponding platform while gray cells represent the non-portability of applications towards the related platform. We find that the OpenMP version is only targeted to CPU platforms and it is non-portable to Nvidia GPUs and Intel GPU. The CUDA version is limited to be executed on heterogeneous architectures of CPU + Nvidia GPU. Both the OpenACC and Kokkos implementations are portable Chapter 5 Metrics Comparison among Programming Models 123 to nine of ten platforms composed of Intel CPUs and Nvidia GPUs whereas the OpenMP offload implementation is only portable to the heterogeneous architecture of GridCL and Gorgon due to lack of compiler support on other platforms. The SYCL version is the single implementation portable to the architecture of Intel CPU+GPU, however, it is not supported by other platforms owing to lack of compilers to support the offload functionality.

Table 5.3: Portability of slabAllNuclides implemented via different programming models across a set of platforms.

Platforms Programming Models

OpenMP thread CUDA OpenACC OpenMP offload Kokkos SYCL

Performance Portability Evaluation

We use a metric of relative performance, application efficiency and a metric of absolute performance, particle tracking rate to perform the evaluation of performance portability with our generic model F . The values of particle tracking rate are exactly those listed in Table 4.8. The data set of application efficiency are measured by dividing the performance (particle tracking rate) of an application on a given platform by the best performance achieved on this platform.

Table 5.4 shows all values of performance metrics that will be applied to F and its variants as the data set C(a, p, h). Based on these values, we can calculate the results of F •C(a, p, h) evaluating the performance portability of different implementations of slabAllNuclides which are illustrated in Table 5.5 and Table 5.6. Note that all blank cells represent that the implementation of a given programming model fails to run on a given platform. These blank cells are corresponding to the gray cells illustrated in Table 5.3 and will be considered as zero values for computation.

From Table 5.5 we can find that for the given platform set and application set, F am is the single metric which renders meaningful results for comparison. By contrast, F gm and F h m both obtain zero values for all implementations of benchmark. In order to make these two metrics capable of retrieving non-zero values, the groups of platform for F am is platform-independent to meet the requirement of adaptability of performance portability expressed in Definition 5.1.17.

F am manages to get meaningful values and the results shows that the priority order is OpenACC > OpenMP thread > CUDA > Kokkos > SYCL > OpenMP offload where the OpenACC version gains the best application efficiency 53% and the OpenMP offload version renders the worst, 7%. F gm and F hm are still unable to evaluate the strict performance portability for the given application set and platform set in our case. All values expressing the relaxed performance portability are calculated by the variants of F as reference.

• Kokkos: set Kokkos::CudaSpace or Kokkos::HostSpace during the compilation.

• SYCL: set sycl::cpu_selector{} or sycl::gpu_selector{}. Now we may apply the metric F towards the data sets of implicit productivity listed above with the intention of evaluating the implicit productivity portability for implementations of slabAllNuclides across a set of supported platforms. From Table 5.8 we find that the estimated values of PS-LOC for each implementation retrieved by F am , F gm , F hm are equal since the platform-specific code distances between different code paths of each implementation are equaling to zero. We can only observe several differences from results of the metric churn. For instance, with reference to the Kokkos version, the estimated value of platform-specific churn aggregated by F am is 8, meaning that given a code path as reference, it requires on average 8 lines of modification to make the Kokkos version successfully running on all supported platforms. It should be noted that F gm and F hm always return zeros if any value in its data set is equal to zero. Thus, both of them are not qualified to capture the feature of churns for each implementation of different code paths.

In general, for implementations which are portable to multiple architectures including the OpenACC, OpenMP offload, Kokkos, SYCL versions, their code divergences and churns are all at a low level, indicating that their maintenance effort is rather little. However, because the code paths of each implementation targeting to different platforms are similar to each other, the retrieved data sets are not suitable to examine different sub-types of F . We can only say that in our scenario, F am is the single metric capable of evaluating the churns of several implementations.