Enhancing Monte Carlo particle transport for Modern Many-core Architectures--第三章，跟踪算法

整体框架图：

3.1介绍

本章内容主要包括两个部分，第一个部分是介绍了两种经典的跟踪算法：基于历史的跟踪算法和基于事件的跟踪算法。第二部分是作者对于两种方法的一些优化。

3.2part 1

3.2.1基于历史的方法

History-based Monte Carlo algorithm
foreach particle history do
    generate particle from boundary condition or source
    while particle not escaped or absorbed do//在粒子未逃逸或被吸收的情况下执行以下循环
        sample distance to collision in material
        sample distance to material interface
        compute distance to cell boundary
        select minimum distance, move particle, and perform event//选择最小距离，移动粒子并执行相应的事件，例如碰撞、散射等
        if particle escaped spatial domain then
            update leakage tally
            end particle history
        if particle absorbed then
            update absorption tally
            end particle history

3.2.2基于事件的算法

Event-based Monte Carlo algorithm
    foreach batch of particle histories (fits in memory constraint) do
        generate all particles in batch from boundary condition or source
        determine next event for all particles (collision, material interface
        crossing, cell boundary crossing)
    while particles remaining in batch do
        foreach event E in (collision, material interface crossing, cell boundary crossing) do
            identify all particles whose next event is E
            perform event E for identified particles and determine next event
               for these particles
        if particle escaped spatial domain then
            update leakage tally
        if particle absorbed then
            update absorption tally
        delete particles absorbed or leaked

3.2.3Thrust（不知道怎么翻译最好）

thrust使用类似STL的模板接口的C++头文件库，是NVIDIA的一个并行算法和数据结构库，为GPU编程提供了灵活的高级接口，无需单独的编写CUDA代码。通过使用Thrust，C++开发人员只需要编写几行代码就可以使用GPU来进行排序、扫描、转换和归约。它实现了CPU和GPU之间的性能可移植性，相关的官方资料如下：

Thrust. (2014). (https://developer.nvidia.com/Thrust)

Thrust - parallel algorithms library. (2018). Retrieved from http://thrust.github.io/

Thrust算法用于在一个批次中实现所有粒子的过程。这些算法执行诸如数据并行映射、归约、收集、分散或扫描等操作，这些操作由Blelloch（1990）定义。Thrust还提供了用于管理GPU设备内存的数据类型，如thrust::device vector和thrust::host vector数据结构，它们类似于C++的std::vector，但在必要时可以在主机和设备之间自动进行内存复制。这些数据类型允许简单的内存管理方案，适用于基于GPU和CPU的架构。

3.2.4 算法细节

这个是使用Thrust对基于事件算法的分组操作的改进，详细方案如下：

对于基于事件的算法，其重点是在所有执行相同事件的粒子之间进行数据并行操作。为了执行这些操作，需要额外的开销来对将要操作的粒子分组，并确定粒子的访问模式。而对于粒子的重新分组可能代价很大，而且与解决传输问题没有直接关系。

Thrust提供了置换迭代器（置换，排列组合的意思），通过索引映射对数据元素进行非对齐访问。使用这种迭代器方案，数据元素在每次操作中不需要被复制到新的位置。然而，这种方法的代价是执行非连续的内存访问来读取和写入信息。

为了使用这种方案在粒子上执行事件操作，需要使用一系列数据并行操作来建立置换迭代器的正确索引映射。该方案的步骤如下：

1. 使用thrust::transform填充一个所有执行事件E的粒子的1和0的模板映射（每个下一个事件为E的粒子将在其索引位置的模板映射中得到一个1）。
2. 使用thrust::reduce计算模板中标记为1的元素的数量（确定将执行事件E的粒子数）。
3. 检查元素数量是否大于0（检查是否有任何粒子执行事件E）。
4. 使用thrust::exclusive scan从模板映射生成索引映射的索引（为每个执行事件E的粒子生成索引）。
5. 分配一个适当大小的新映射（用于保存执行事件E的所有粒子的索引）。
6. 将扫描中的索引散布到新的索引映射中（将独占扫描生成的索引减少到仅包含足够执行事件E的粒子的映射中）。
7. 在置换迭代器循环中使用新的索引映射（将索引映射与置换迭代器结合使用，允许循环只对索引映射中选择的粒子执行操作）。

3.2.5实现

在这个段落中，作者使用 Nvidia CUDA 编程模型实现了 基于事件的ALPSMC版本，并使用 Nvidia C++ Thrust 库。ALPSMC 的 Thrust 实现利用了数据并行操作和 Thrust 数据类型来管理内存。相同的 Thrust 基于事件的实现可以使用 CUDA 在 GPU 上编译，也可以使用 OpenMP 在 CPU 上编译，从而实现了在不同平台上的可移植性。

在 ALPSMC 的本机 CUDA 实现中，我们发现在构建各种映射时继续使用 Thrust 算法非常有用，并使用 CUDA 直接启动事件内核，而不是调用 Thrust::foreach。 Thrust 和 CUDA 实现的 ALPSMC 提供了与原始基于历史记录的实现完全相同的物理结果。

本研究的 CUDA 实现与 Thrust 实现中的算法相匹配。性能差异来自于本机 CUDA 编程提供的能力，这些能力无法通过 Thrust 实现。直接使用 CUDA 可以在内核级别实现更精细的控制，并且可以重要地访问不同的内存空间，如 GPU 共享内存。CUDA 实现包括调度算法，以优化每个内核调用的活动线程数量。此外，CUDA 实现包括使用不同可用内存空间，例如常量和共享内存。例如，蒙特卡罗粒子最初在 GPU 全局内存中分配，然后在内核中的所有操作中复制到共享内存。所有问题常量，如截面和平均弦长值，都放置在 GPU 常量内存中。在某些条件下，这些优化对 GPU 内核的性能有显著影响。

3.2.6初步结果

实验目的：通过改变问题规模（MC particlehistories）和实现方式（Thrust or CUDA），来探究改进后的基于事件的算法域最初的串行基于历史算法的加速情况，结果以加速比表示。

实验案例：Case 1a， 空间域为10cm

实验平台：

计算机名称	host core（主机核心）	Accelerators（加速器）
Rzgpu计算机	Intel Xeon Westmere 2.8GHz	70个NvidiaTesla M2070 GPU
Max计算机	Intel Sandy Bridge 2.6GHz	Nvidia Tesla K20X GPU
Rzhasgpu	Intel Xeon Haswell 3.2GHz	Nvidia Tesla K80 GPU(由于没有实现代码的MPI变体，大概只能利用其一半的计算能力）

实验结果如下：

图中结果为对比之前串行的基于历史算法的基于事件的算法的加速比。

纵轴表示两个版本CUDA和Thrust，都在三台实验平台进行了测试；

横轴显示了不同的粒子历史规模。

结论：由上表可以看出Thrust的加速比大约从1.4到2.6，这个加速比与预期的GPU硬件的加速比不符。而使用CUDA明显大于Thrust，加速效果大概为Thrust三四倍。对于基于事件的MC算法，直接使用CUDA会更加有效。

实验二：

实验背景：基于事件的Thrust事件代码实现还使用了OpenMP编程，用于在主机CPU上，也就是它具有可移植性。

实验目的：研究基于事件的Thrust版本的缩放情况

实验变量：

Thrust事件版本、CUDA事件版本、串行历史版本

粒子历史数量

实验平台：Rzhasgpu的Intel Xeon Haswell CPU，使用16个线程/核心

实验结果：横轴-粒子历史数量、纵轴-挂钟时间

实验结论：从图中可以看出，大约从10^5开始，基于事件版本挂钟时间少于串行历史版本。粒子数量越高，CUDA版本>Thrust版本>串行历史版本，见图后半段。

3.3part2

这一部分主要讲述了提升ALPSMC再GPU上的速度而做出的两个重要更改。

1.一个是将粒子的数据结构从结构数组改成的数组结构；

2.一个是基于时间算法粒子移除方案；

最后讲述了其他的一些小修改。

3.3.1结构数组与数组结构

结构数组：也就是每个结构或者说对象作为一个整体存储在数组里面。也就是说这个本质上是一个数组，而数组中的每个元素是一个结构或者说对象。

假设由一个粒子结构，其中包含位置、速度、方向等属性，每个粒子的所有属性都会连续存储在一个结构里，然后将一个一个这样的粒子结构存储到一个数组中。

数组结构：将同一类型的所有数据元素组织在一起为一个数组，不同的相同类型的数据元素数组组合成为一个结构。也就是说这个本质上是一个结构，结构中的一些元素是数组。

继续使用粒子的例子，你会有一个数组专门存储所有粒子的位置，另一个数组专门存储所有的粒子的速度，以此类推。

剩下部分请看本文资源