Enhancing Monte Carlo Particle Transport For Modern More-core Architecture--背景和相关工作-CSDN博客

本文链接：https://blog.csdn.net/weixin_52190646/article/details/136748931

background and relative work

2.1 介绍

超级计算机与蒙特卡洛粒子传输的介绍。最近二十多年来，超级计算机的架构发生剧烈变化，超级计算机的架构由更快更复杂的多核CPU向着更多数量、更简单更缓慢的简单处理器发展，也就是CPU架构向着专用处理器发展，例如GPU，TPU等。这种变化可以使计算机能够高度并行计算。本节分为四个部分：

1.蒙特卡洛粒子传输的背景

2.蒙特卡洛粒子传输计算的最新研究的相关工作的讨论

3.专注于蒙特卡洛在GPU上的计算

4.多种架构上的可移植性

2.2 什么是蒙特卡洛粒子传输？

Eckhardt (1987)提供了斯坦·乌拉姆和约翰·冯·诺伊曼关于纸牌游戏中的未发表对话，这成为启动蒙特卡罗输运方法的基础。"关于蒙特卡罗方法，我最初的想法和尝试是在1946年我病愈休养时玩纸牌时产生的。当时我在思考这样一个问题：52张纸牌的康菲尔德纸牌布局能成功完成的几率是多少？在花费了大量时间尝试用纯组合计算来估计之后，我想知道是否有比“抽象思维”更实用的方法，比如将它布局一百次，然后简单观察和计数成功的游戏次数。随着快速计算机新时代的开始，这已经可以想象到，于是我立即想到了中子扩散等数学物理问题，更一般地是如何将某些微分方程描述的过程转化为等效的可解释为一系列随机操作的形式。后来... [在1946年，我] 将这个想法描述给了约翰·冯·诺伊曼，我们开始计划实际计算。" - 斯坦·乌拉姆 1983年约翰·冯·诺伊曼对斯坦·乌拉姆的想法产生了兴趣，并概述了如何解决裂变设备中的中子扩散和增殖问题。自那时以来，Eckhardt（1987）告诉我们，蒙特卡罗方法一直是解决许多中子输运问题的主要方式。

2.2.1定义：

在《中子输运的计算方法》中，Lewis和Miller（1993）将蒙特卡罗输运描述为利用随机数发生器模拟一定数量的粒子历史。对于每个计算的粒子历史，会生成随机数并用于采样描述粒子可能经历的不同物理事件的概率分布，比如散射角度或者碰撞之间的距离等。Lux和Koblinger（1991）在他们的书《蒙特卡罗粒子输运方法：中子和光子计算》中进一步扩展了之前的定义：“在蒙特卡罗方法的所有应用中，都会构建一个随机模型，其中某个随机变量（或多个变量的组合）的期望值等于待测物理量的值。然后通过对上述引入的随机变量进行多次独立采样的平均值来估计这个期望值。为了构建一系列独立样本，会使用遵循待估计变量分布的随机数。”（第5页） Lewis和Miller（1993）解释了估计一个数量的数学形式如下：

$\hat{x}=\frac{1}{N}*\sum_{1}^{n}$ $X_{n}$

，其中xn表示该数量的第n次历史的贡献。对于蒙特卡罗方法，我们对每个粒子历史的Xn进行统计，以计算期望值 $\hat{x}$ 。一个非常重要的问题是估计值ˆx与真实值 x之间的比较。事实证明， $\hat{x}$ 的不确定性随着粒子历史数量的增加而减少，并且通常按照 $N^{-1/2}$ 的渐近比例下降，其中N是粒子的数量。如图所示，横坐标是粒子数量，而纵坐标是是 $N^{-1/2}$ ，代表 $\hat{x}$ 的不确定性。

2.2.2对应方程

Monte Carlo中子传输方程要解决的的方程是线性化波尔兹曼传输方程。这个方程的每个部分使用了大量的粒子来创建准确的估计。由于波尔兹曼方程可以用不同的方式写下，我们关注由Lewis和Miller（1993年）描述的变化。方程和组件的描述如下：

$\Psi (\vec{r},\hat{\Omega },E,t)$ 表示角通量，指单位事件内通过单位立角体的粒子流量或能量流量，是中子在不同能量和角度下的流量

$\sum ^{}_{t}(\vec{r},E)$ 表示所有粒子反应的宏观总截面，表示粒子与介质中的所有原子核或电子发生相互作用的总概率
$\sum_{s}(\vec{r},E{'}\rightarrow E,\hat{\Omega} {'}\cdot {\hat{\Omega}}))$ 表示粒子散射的宏观截面，表示粒子从一个能量状态散射到另一个能量状态的概率

$\sum_{f}(\hat{r},E)$ 表示来自裂变反应的粒子产生的宏观截面，表示粒子参与裂变反应的概率

$S_{ext}(\vec{r},\hat{\Omega},E,t)$ 表示外部源，表示外部对系统的输入

$\chi (E)$ 表示裂变的次级粒子谱

$v(E)$ 表示每个裂变产生的平均粒子数

$v$ 是粒子的速度； $\vec{r}$ 表示位置矢量，表示粒子在空间中的位置； $\hat{\Omega }$ 是指立体角的单位矢量，代表粒子传播方向；E是粒子的能量；t是时间。

对于这个公式有想法的可以参考：https://zh.wikipedia.org/wiki/%E7%8E%BB%E5%B0%94%E5%85%B9%E6%9B%BC%E6%96%B9%E7%A8%8B

2.2.3 算法

主要介绍了基于历史模拟算法(history-based approach)，它跟踪单独粒子的历史，直到模拟了一定的数量的粒子。为了模拟粒子，我们必须计算在粒子发生反应前每个粒子所运动的距离，并且与可能发生反应的距离进行对比。最短的反应距离被选择，随之更新基于运动和反应发生的距离粒子和计数数据。

如下展示了基于历史的算法：

parse input//分解输入

for each Cycle do{

Cycle_init ( ) {//周期初始化阶段
    source in particles//初始化粒子，生成初始粒子数量
    population control//粒子数量控制，人口控制策略，保证粒子不会过多或者过少
}
Cycle_tracking ( ) {//周期跟踪阶段
    for all particles {//对每个粒子进行循环-----并行化发生在此
        do {
              compute distance to census//计算到下一普查点的距离，就是粒子不发生反应的距离
              compute distance to facet//计算到面交叉的距离，到网格的距离
              compute distance to reaction//到达下一个反应（碰撞或者吸收）的距离
              do segment with shortest distance increment tallies//执行最短距离的事件，增量计数
             } until census, absorbed, escaped//直到它们被吸收、逃逸或者到达普查
      } 
}
Cycle_finalize ( ) {//周期结束阶段
    reduce all tallies //该函数处理时间步结束时的一些后勤工作，如计算所有统计数据的全局归约，得到最终的结果
}}

gather tallies//收集数据

初基于历史的算法还有基于事件的算法，将在2.4提及。

2.3蒙特卡洛的研究现状

蒙特卡洛传输问题的研究历史、改进以及最近的研究方向。

理解这一发展历程以及为之设计的计算机可以帮助指导对更近期努力的分析。最近的蒙特卡洛传输研究主要涉及以下三个主题之一：

（1）GPGPU计算，

（2）并行算法改进（不涉及GPGPUs）

（3）物理学改进。

关于与GPGPU相关的研究的回顾在第2.4节（最新技术：GPU研究）中展示。本节主要关注非GPU蒙特卡洛传输研究，即在CPU架构上的并行性能、并行负载平衡、核数据查找的优化以及方差缩减技术。

2.3.1 并行性能

自蒙特卡洛粒子传输诞生六十多年来，随着蒙特卡洛应用在每一代新架构上的性能最大化，其增长势头十分迅猛。1947年开发的第一个模型需要五个小时才能计算100次碰撞，而今天这一任务可以在毫秒内完成。在1940年代和1950年代，蒙特卡洛代码主要使用最早期的计算机上非常低级别的语言编写。1960年代至1980年代，随着计算能力的提升和代码功能的增强，蒙特卡洛代码的能力显著增强。1980年代，蒙特卡洛代码开始采用并行/向量机器。1990年代，蒙特卡洛代码变得更加常见，并行性增加到通过PVM（2011年）或消息传递界面（Clarke, Glendinning和Hempel（1994））实现的数百到数千处理器。2000年代，Brown（2011年）解释说多核处理器意味着线程化变得更加常见，将本地和全局并行性混合使用可达到数万处理器。这种计算能力的增长也可以从内存访问风格的角度进行分类。早期系统几乎完全是共享内存环境。然后，分布式内存系统变得流行起来。最后，混合并行系统（即同时使用分布式和共享内存并行性的系统）变得流行起来。以下讨论将围绕这三种并行性展开。

共享内存系统的相关研究。共享内存系统指的是所有处理器都可以访问相同内存空间的机器或模型。

向量机器上的相关研究。 描述了向量集上的一些研究，尤其是基于事件算法的提出。

多线程架构系统 除向量机之外的共享内存系统的研究，引出现代系统中运用到的共享内存的思想，也就是使用open MP线程模型共享内存处理。

分布式内存系统 超级计算机由向量计算转变为分布式内存计算，讲述了一些相关研究

分布式+共享内存 主要讨论共享内存和分布式内存两种模型的结合，一个主要例子便是open MP和MPI的结合

2.3.2 负载平衡和域分解

域是问题空间或者几何图形的一部分，取决于应用的场景。相当于把大问题化为小问题。域的两种模型分别是域复制域域分解。域复制中的负载平衡简单的将粒子分裂给对应的处理器，负载平衡经常和域分解一块讨论

什么时候域负载平衡？ 把执行负载平衡的计算代价和负载不平衡的代价结合考虑。

第一步是通过比较当前并行效率（εC）和如果处理器重新分配负载后的并行效率（εLB）来计算加速比。第二步是通过使用前一周期执行时间（τP hys）、加速比（S）以及计算负载平衡本身的时间（τLB）来预测运行时间。最后一步是比较有负载平衡和无负载平衡时的预测运行时间，以确定该操作是否值得进行。

扩展域分解 涉及到网格域分解 大规模并行计算机的负载平衡

2.3.3核数据

核数据主要是材料在各种条件下域粒子发生相互作用的信息。包括所有可能的反应以及原子碰撞的围观横截面数据，同时也包括宏观界面信息。

核数据的查找。核数据存储在表格当中，通常在转盟的库中解释和处理。存储和查找和核截面数据有两个方法，连续能量方法和多组方法。连续能量方法花费时间久，但查找准确。多组方法，查找时间短，但不是很准确。而对于核数据查找的研究重点便是如何加速在给定截面的搜索？线性搜索、二分搜索、哈希搜索等。现有的研究由哈希化、联合搜索、分数级联等方法。

2.3.4方差缩减技术

在蒙特卡洛传输的解是以统计数据的形式输出的，所以减少方差有利于提高准确性和计算最后的解。没有缩减方差会导致花费大量的时间去寻找解。

以下是一些常用的方差缩减技术：

共同随机数（Common Random Numbers）：通过比较两个或多个不同的配置，而不仅仅是单个配置，来实现方差缩减。这也被称为样本的相关性。通过确保问题的所有配置都使用相同的随机数来找到解决方案，可以引入一种正相关性的元素，从而实现方差缩减。
对偶变量法（Antithetic Variates）：对于每个采样路径，采用其对偶路径，例如对于给定的路径 {ε1, ..., εM }，还会采用路径 {−ε1, ..., −εM }。这种方法减少了所需的样本数量，并减小了采样路径的方差。
控制变量法（Control Variates）：通过使用已知量的信息来创建相关系数，以减小未知量中的误差。这种方法等价于解决一个最小二乘系统，因此通常被称为回归采样。
重要性采样（Importance Sampling）：这种方法涉及从不同于所需分布的分布中生成样本，并估计特定分布的属性。通过更频繁地对重要值进行采样，以及对不重要的值进行较少的采样，强调了重要性值。通常通过分裂或俄罗斯轮盘方法实现这一点。
分层抽样（Stratified Sampling）：通过将总体成员分成同质群体，然后进行采样，来减小采样误差，并且产生的加权平均值比总体的算术平均值具有更小的变异性。何况

2.4 GPU的研究的最新技术

如图为Cpu和GPu的示意图，CPU拥有多个核，核少但功能强大，用于顺序处理。

而GPU具有大规模并行架构，包括数千个更小更高效的核心，能够同时处理多个任务，处理多个任务时，吞吐量大。

2.4.1 蒙特卡洛传输在GPU上的应用

本小节从几个关键领域（准确性、性能和算法选择）来讨论了相关的研究，随着时间的流逝，硬件也在不断进步。

准确性： 虽然随着时间的发展，现在GPU上的准确性已经好的多，但还是需要考虑硬件的准确性。尤其是单精度浮点数的准确性、以及其在CPU和GPU结果的不同、以及IEEE-754规范的合规性。

性能：主要探讨了从CPU硬件迁移到GPU硬件的性能的变化，并根据不同的研究进行了不同的性能比较。（例如：光子传输、中子传输、伽玛射线传输、耦合的电子光子传输等）

算法：基于前面的相关性能的算法，详细的讨论了一些重要的算法：粒子块算法、基于事件的算法。

2.4.2蒙特卡洛和医学

在医学领域，Monte Carlo传输经常被忽略。放射输运用于对患者进行剂量估计，讨论了三个医学蒙卡运输应用的描述：GMC中的电磁Monte Carlo运输、质子疗法中的gPMC、DPM中的电子-光子输运。

2.4.3 蒙特卡罗与光线追踪

蒙特卡罗传输涉及一个关键且经常需要大量计算的步骤，即确定粒子是否会与任何背景几何体发生相互作用或者转移到不同的材料区域。这个过程类似于光线追踪，在计算机图形学中用于通过追踪光线路径穿过像素并模拟其与虚拟对象的相互作用来生成图像。

2.4.4基于事件的蒙特卡洛算法

Algorithm 3: The basic iteration event
for event n = 0, 1, 2, ... do//在每次迭代中，从事件n开始，依次进行操作
    Fetch Γ^n//获取事件n对应的数据集Γn
    Perform free flight analysis://获取自由飞行分析：
        gather the cross section data and geometry data tabulated by particle,
//收集粒子表格化的截面数据和几何数据
        Σ ← S,//截面数据赋给Σ 
        ρ ← R;//几何数据赋给ρ 
        using Σ, sample a vector of distances to collision, dc//使用求和，生成一组到碰撞的距离组成的向量dc
        using ρ, determine vector of minimum distances to boundary, db//使用ρ，确定到边界的最小距离向量db
        determine the minimum distances to the end of event,//决定到事件结束的最小边界
            d_min = min[dc, db];//两者取最小值
        update the particle coordinates,//更新粒子坐标，
            r^(n+1) = r^n+ Ω^n · r_min
        Perform collision analysis://执行碰撞分析
            gather particle attributes,//收集粒子属性
                Ω ← Γ^n, E ← Γ^n;//将Γn赋给Ω 和E 
            evaluate collision physics for new direction cosines and energies,//评估新方向矢量和能量对应的物理碰撞
                Ω′ ← Ω, E′ ← E
            scatter new particle attributes back into bank,//散射型例子属性回数据库
                Ω′ ← Γn, E′ ← Γn
        Perform the boundary analysis://执行边界分析
            gather particle zone indices Z,//收集粒子区域索引Z
                Z ← Γ^n
            determine new zone indices,//确定新的区域索引
                Z′ ← Z :
            scatter new zone indices back into bank.//散射新的区域粒子缩影回数据库
                Z′ → Γ^n
        Update the particle bank,//更新粒子数据库
            Γ^n ⇒ Γ^(n+1) (with Ln+1 particles)
            (e.g. compress out terminated particles).
        If L[n+1] != 0, continue