面向天河2A系统的基于蒙特卡罗方法的粒子输运异构协同计算 *

李 彪，刘 杰，2

(1.国防科技大学并行与分布处理国家重点实验室，湖南 长沙 410073,2.复杂系统软件工程湖南省重点实验室，湖南 长沙 410073)

1 引言

1.1 粒子输运简介

输运理论是研究微观粒子在介质中迁移统计规律的数学理论。这里“微观粒子”指的是中子、光子、分子、电子和离子等。输运理论研究大量粒子在空间或者某种介质中运行时，由于各粒子位置、动量和其他特征量的变化而引起的各种有关物理量随时空变化所表现出来的非平衡统计运动规律[1]。

19世纪中后叶的分子运动论导致粒子输运理论发展的开端。1872 年，Boltzmann 推导出了Maxwell分子速度分布函数随时空变化的非线性积分-微分方程，该方程也常被称为粒子输运方程或者Boltzman方程[1]。粒子输运方程描述的分布函数具有7个独立变量，包括3个空间变量、2个方向角度变量、1个能量变量和1个时间变量。由于实际情况的复杂性，通常不能采用解析方式进行直接求解，只能通过数值模拟求解。

目前常用的求解粒子输运问题的数值方法分成2类[2]：一类是“确定性方法”，这类方法通过对变量进行离散，将粒子输运方程转化为一组线性代数方程，然后通过求解此代数方程组来获得精确解或近似解；另一类方法是“蒙特卡罗MC(Monte Carlo)方法”，也称非确定性方法或概率论方法。它是基于统计理论的数值方法，对所要研究的问题构造一个随机概率模型来模拟大量粒子的运动，以获得感兴趣的物理量分布。

本文研究的是关于粒子输运中的非确定性方法。

1.2 蒙特卡罗方法在粒子输运上的应用

蒙特卡罗MC方法亦称随机模拟法，是20世纪40年代中期电子计算机诞生后发展起来的一门计算科学[3]。它是通过计算机来对中子行为进行随机模拟的数值方法，是一种以概率统计理论为指导的非常重要的数值计算方法。

通常蒙特卡罗方法可以粗略地分成2类[3]：一类是所求解的问题本身具有内在的随机性，利用计算机直接模拟这种随机的过程。例如，在中子输运理论中，要求建立单个中子在给定几何系统中的真实运动历史，通过对大量中子历史的跟踪，得到充分的随机试验值(或称抽样值)，然后用统计的方法得出随机变量某个数值特征的估计量，用此估计量作为问题的解。另一类是所求解问题可以转化为某种随机分布的特征数，比如随机事件出现的概率，或者随机变量的期望值。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的多维积分问题。

1.3 基于MC方法的异构协同计算

由于功耗和散热问题日益突出，单块芯片的性能已经不能单纯地从提高频率来获取提升，摩尔定律将会消亡[4]，于是工业界逐渐转向多核处理器。但是，随着技术的发展，多核处理器中每个内核的设计越来越复杂，并且频率较高，功耗较大，严重阻碍了多核处理器中核数的扩展。为了进一步提升性能，同时满足功耗和散热的要求，出现了由众多相对简单的内核构建的众核协处理器，如 NVIDIA公司推出的 GPU[5]、Intel公司推出的MIC协处理器[6]和本文所探讨的国产加速器Matrix2000等。与多核处理器相比，这些众核加速器的内核相对简单，且频率较低，但数量众多，能够以合理的功耗实现更高的性能。

近十年来，异构加速器快速发展，TOP500排行中排名靠前的超算系统的体系结构均采用了异构协同模式，许多研究人员对已有的求解粒子输运问题的方法在各种加速平台上做了许多相关的研究工作。

杨博[7]提出了一种针对CPU/GPU混合异构系统的深穿透粒子输运MC模拟异构协同算法，针对GPU的计算和访存特点，给出了一种基于粒子数的任务划分方法和高效并行数据结构以及线程间的并行规约方法。相比运行在Intel Xeon X5670上的MC程序，该算法获得了3.4的加速比，并在TianHe-1A的64个结点进行了测试，结果表明该算法具有良好的性能和可扩展性。

崔显涛[8]基于MC方法，提出了一种面向CPU/MIC混合异构系统的粒子输运并行算法，针对MIC访存特点提出了适于程序并行的高速数据结构和基于静态分配的任务划分方式。相比运行在CPU端的程序串行程序，改进的MC程序获得了8.6倍的加速比。

在多能群MC粒子输运领域，近些年来有许多研究和探索,Bergmann 等人[9]基于GPU开发了一个MC粒子输运框架WARP(Weaving All the Random Particles)。WARP的目的是使以前的基于事件的输运算法适应新的基于GPU硬件而开发的连续能量Monte Carlo中子输运代码。Hamilton等人[10]基于Shift程序包开发了一个在GPU上运行的连续能量Monte Carlo中子输运求解器，并把程序移植到了超算系统Summit上，在耗尽燃料的SMR(Small Modular Reactor)核心配置上，进行了1 024个结点的弱扩展测试，结果显示其并行效率接近93%，这说明GPU对性能提升显著。

本文的工作与文献[8,9]的接近，都是只考虑单群的MC程序，不过本文的异构算法设计是基于不同于GPU和MIC体系结构的国产加速器Matrix2000，且研制的程序可以移植到大规模结点上，另外，还针对MC程序的数据通信进行了优化，取得了显著的性能提升。

2 粒子输运MC模拟

MC模拟是通过对大量中子历史的跟踪来获取大量的随机试验结果，并对这些结果进行统计方法处理，最终以得到的统计量作为问题的解。 所谓一个中子的历史，是指该中子从源出发，在介质中随机游走，经过各种核反应作用，直到中子历史结束或称中子“死亡”。所谓“死亡”是指中子被吸收、穿出系统、被热化，或达到能量权下限或时间上限。其中时间、能量的载断是无条件的，而权截断是有条件的，由俄罗斯轮盘赌决定[11]。

图1是一个中子历史的循环过程。初始时，从粒子源出发的一个粒子的位置、方向、能量等参数是确定的，随后粒子的运动方向和与原子核的碰撞类型服从特定的概率分布，在经历了一系列碰撞过程之后，粒子历史趋向于结束。有几种不同类型的结束，比如达到时间与能量界限、逃脱出系统或达到权下限。在这一个过程中，下一次碰撞仅与当前碰撞后粒子的位置、能量以及方向相关，完全独立于先前的碰撞。 因此，这一过程是一个典型的马尔可夫过程。所以，只要知道粒子与原子核碰撞的规律，那么粒子的轨迹就可以用蒙特卡罗方法正确地模拟出来，从而得到所关心的解[12]。

Figure 1 Simulation of particle transport process图1 模拟粒子输运过程

3 通信模式优化

3.1 串行阻塞通信模式

通过测试MC程序的强扩展性发现，当进程数大于256时，效率急速下降，达到1 024进程时基本上已经没有加速效果了，这对移植MC程序到大规模系统上产生了巨大障碍。经过仔细分析程序的通信模式发现，由于计算结果数据收集的通信模式是主从模式，导致等待时间过长。

首先0号进程初始化任务，接着把数据发送给其他所有从进程，从进程接收到数据后开始计算，并在从进程的计算结束后设置同步栅栏，直到所有从进程的计算都结束后，0号进程才开始收集从进程发送过来的结果数据。数据的收集过程中所有从进程都向主进程阻塞发送数据，主进程进行阻塞接收，先到的信息先被处理。对于主进程来说这是一个串行处理过程，此通信方式会造成通信瓶颈，核数越多，对程序性能的影响越大。图2展示了P个进程的通信过程。

Figure 2 Serial data collection mode图2 串行数据收集模式

3.2 优化后的二叉树通信模式

考虑到MC程序0号进程对收集到的数据只进行简单的统计处理，所以可以通过逐层进行两两进程通信，实现局部统计处理，最终把统计结果汇总到0号进程。本文采用二叉树通信模式来实现这种通信过程，假设总进程数为2K+res，其中,K,res都是大于或等于0的整数，且满足条件0≤res<2K，为了满足二叉树的通信模式，排在前面的2*res个进程先就近奇偶进程号两两通信，并把通信过程中已发送过信息的进程剔除，不参与下次通信，从而下一阶段通信的进程数变为2K。接着对2K个进程号重排，前2*res个进程号除以2，其余进程依次减去res，得到新的2K个进程的标识号，下一阶段就是进行二叉树通信，算法1描述了具体的通信过程。为了方便描述，图3展示了7个进程的通信过程，括号里的数字表示重排后的进程标识号。

Figure 3 Binary tree data collection mode图3 二叉树数据收集模式

改进后的数据收集模式，通信复杂度由2K+res减少为log (2K+res)，极大地减少了通信时间，也避免了大规模通信时从进程同时向主进程发送数据导致的程序阻塞。

4 基于CPU/Matrix2000的异构协同算法

4.1 Matrix2000加速器

天河2A系统中，每个结点由2颗Intel Xeon 微处理器和2颗Matrix2000加速器组成，如图4所示。每个Intel Xeon微处理器包含12核，工作频率为2.2 GHz，采用英特尔Ivy Bridge微架构，峰值性能0.211 2 TFLPOS。每个Matrix2000加速器包含128核，由4个超结点组成，每个超结点包含32个计算核，其中超结点支持64核超线程技术，峰值性能2.457 6 TFLOPS@1.2 GHz，有8个DDR4内存通道，支持×16 PCIE 3.0 EP工作模式。

算法1二叉树通信伪代码

//假设总进程数为2K+res，0≤res<2K

1mask=1;

2if(myid< 2*res)then/*前2*res个进程先通信，使下一阶段通信进程个数为2K*/

3if(myid%2!=0)then//奇数号进程

4dest=myid-1;

5new_id=-1/*进程号置为-1，不参与下次通信*/

6myid进程SENDmsgtodest;

7else//偶数号进程

8src=myid+1;

9new_id=myid/2;//重排后新的进程号

10myid进程RECVmsgfromsrc

11endif

12else

13new_id=myid-res;/*其余的进程重排后的进程号*/

14endif

15if(new_id!=-1)then/*重排后的2K个进程进行通信*/

16while(mask<2K)do//二叉树通信模式

17if(mask&new_id)then/*获取重排号之后的偶数进程号*/

18newsrc=new_id|mask;

19if(newsrc

20src=newsrc*2

21else

22src=newsrc+res

23endif

24 RECVmsgfromsrc

25else

26newdest=new_id& ～mask

27if(newdest

28dest=newdest* 2

29else

30dest=newdest+res

31endif

32 SENDmsgtodest

33 exit//发送过消息后，此进程就退出通信

34endif

35mask= 2*mask

36endwhile

37endif

Figure 4 Structure of Tianhe-2A system node图4 天河2A系统结点结构

4.2 异构编程模型

天河2A系统支持3种异构通信模式OpenMP4.5、ACL和BCL，其中BCL是一种简单高效的对称传输接口，CPU和协处理器之间数据通信利用PCIE总线进行传输，底层通过SCIF来实现。BCL相比OpenMP4.5更底层，程序移植也更复杂，但是传输速率更快，移植后的程序灵活性更好。基于BCL接口的异构程序需要编译2套程序，CPU和加速器分别编译一套程序，2套程序分别在CPU端和加速器端同时运行，其中运行时加速器端的程序需要通过ACL传送到加速器端并启动。

图5给出了基于CPU/Matrix2000的MC程序的异构模式流程。首先，CPU端0号进程启动MPI进行进程初始化，0号进程负责读取文件数据，并把总计算任务按照粒子数均等的方式分配给从进程对应的子任务，再将子任务的初始化数据传输给对应的从进程。每个从进程控制一个Matrix2000超结点，通过ACL把加速器端的程序传输到Matrix2000端并启动程序，利用BCL使CPU端与 Matrix2000超结点建立连接，CPU端把初始化数据传输到Matrix2000端，加速器端完成CPU端数据接收后就利用多线程技术并行跟踪每个粒子历史，直到所有粒子历史全部完成计算，加速器统计结果数据并把数据通过BCL传输给对应的CPU进程，再由MPI实现进程间的传输。整个异构控制的过程中，CPU主要负责传输数据，Matrix2000主要负责计算。

Figure 5 Flowchart of heterogeneous logic 图5 异构模式逻辑流程图

4.3 OpenMP线程级并行

Matrix2000中一个超结点包含32核，支持64核超线程技术，为了充分发挥Matrix2000的性能，通过OpenMP指令实现线程级细粒度并行，算法2描述了具体算法伪代码。加速器端接收CPU端发送过来的数据，接着做初始的工作，针对粒子输运模块进行多线程并行。每个线程跟踪一个粒子，完成当前粒子跟踪后，进入time_to_stop判断是否所有粒子数都已完成，如果未完成，当前线程临界更新粒子数变量nps++，更新后的第nps个粒子分配给当前线程进行跟踪模拟。这种调度方式类似迭代块大小等于1的dynamic调度。

算法2Matrix2000加速器上的OpenMP线程级并行

1 假设线程数为T，加速器分配到M个粒子，粒子编号区间为(nps,nps+M-1);

2 #pragma omp parallel;/*开启多线程，每个线程跟踪一个粒子，完成后继续跟踪下一个粒子*/

3forall threadsdo

4 各个线程初始化自己的私有变量;

5omp_set_lock(&lock);/*设置锁，当前只有获得锁的线程才能更新共享数据*/

6 更新共享变量;

7while(!time_to_stop())do/*判断当前是否完成M个粒子数追踪*/

8nps=nps+1;/*未完成，更新nps，追踪第nps个粒子*/

9 更新共享变量;

10omp_unset_lock(&lock);/*释放锁，使得多线程可以调用hstory并行追踪粒子*/

11 callhstory();/*线程各自调用hstory()函数，互不干涉，完成输运过程*/

12omp_set_lock(&lock);/*设置锁，进入time_to_stop()中访问共享数据*/

13endwhile

14 更新共享变量;/*完成M个粒子数追踪，进行数据更新*/

15omp_unset_lock(&lock);

16endfor

相比Intel MIC协处理器Offload模式多线程技术以及OpenMP4.5提供的异构模式，本文的异构模式中的CPU与Matrix2000加速器并不是主从关系，而是对等的关系。Matrix2000加速器端启动的是单独的一套程序，所有的数据传输和初始化工作在启动多线程之前已经完成，所以Matrix2000上的OpenMP不需要通过编译指导语句来传输数据和设置各种变量属性，使得程序结构设计更简单、更灵活。

5 实验及结果分析

5.1 测试环境及参数

测试平台为天河2A系统，由于ACL和BCL指令只提供C/C++的接口，本文需要对MC程序插入C语言的通信控制接口，CPU端的程序采用Intel的编译器和基于高速网的mpich3.2进行编译；加速器端采用自定义交叉编译器进行编译，支持OpenMP指令。

单结点硬件配置为2颗Intel E5-2692 v2 @2.20 GHz(12 核/颗)+8个超结点Matrix2000，每个超结点通过超线程技术实现45核多线程并行。 MC程序的输入文件参数除问题规模随具体实验设置外，均保持默认值。

5.2 通信优化测试

考虑到数据的通信仅在CPU端进程之间进行，而并不涉及到加速器，所以本节仅在CPU上进行对比测试，表1给出了通信优化前后的测试结果。

Table 1 Test results before and after communication optimization

分析图6可以发现,原始程序测试结果中，512核是拐点，小于512核时程序还有加速效果，大于512核时测试时间反而急速上升，已经没有加速效果，表明此时计算能力不是瓶颈，数据通信占主要耗时。相反，优化了通信模式的程序测试时，运行时间随进程数递增逐渐下降，没有出现进程数越多测试时间反而上升的现象，当大于2 000核时耗时趋于稳定，这是由于粒子数计算规模受限，增加进程数的优势已经不明显导致的。可以设想如果计算规模增加进程数越多，优化通信后的程序会比原始程序效果更好。

Figure 6 Comparison of two communication modes图6 通信模式对比图

5.3 异构算法测试

本节对基于CPU/Matrix2000的异构MC程序进行弱扩展测试，每个结点启动10个CPU进程，其中8个进程分别和结点内的8个Matrix2000超结点进行异构计算，剩余2个进程协同计算，每个超结点使用48核多线程，表2给出了弱扩展测试结果，图7给出了相应测试结果的柱状图。

Table 2 Weakly scalable test

Figure 7 Bar and curve chart of weakly scalabletest result on Tianhe-2A图7 天河2A弱扩展测试结果柱状图

表2说明弱扩展到45万核时，相对5万核并行效率保持在 22.54%。 分析图7可以得到，虽然研制的程序可以运行到45万核，但是可以看到，随着核数增大并行效率的下降趋势很明显。结合程序与加速器体系结构分析，导致并行效率下降迅速的原因主要有2点：首先，MC程序的粒子输运部分的计算不是数据密集型，在中子的历史过程中充满了大量粒子条件状态的判断，程序的逻辑控制语句居多，这与Matrix2000加速器适合程序逻辑控制简单、数据计算密集的条件相冲突；其次，在程序的OpenMP多线程部分，虽然开启了多线程技术，但是粒子历史的模拟追踪过程中各个线程都涉及更新共享数据，需要设置临界区，在临界区内同时只能有一个线程更新共享数据，其他线程只能等待，这会影响多线程的并行效率。

6 结束语

本文在现有粒子输运蒙特卡罗模拟程序的基础上，提出了一种面向CPU/Matrix2000异构系统的粒子输运异构协同算法，基于天河2A系统的异构通信模式BCL和ACL，针对CPU和加速器各自设计了一套不同的代码，在CPU端通过ACL实现加速器端代码的启动，利用BCL进行CPU与Matrix2000之间的通信，进而提出了一种CPU与加速器Matrix2000之间的简单高效的对称通信模式。优化了原MC程序的串行数据收集通信模式，提出了新的二叉树通信模式，极大地减少了通信时间，加速比可达17.7。通过优化通信模式，以及基于MPI-ACL-OpenMP编程框架，本文实现的基于CPU-Matrix2000异构协同计算的并行程序，可以弱扩展到45万核，相对5万核并行效率保持在22.54%。作为未来工作，将探索如何继续增强该异构程序的扩展性，以更好地发挥出HPC平台的算力优势。