基于RMC的蒙特卡罗程序性能优化*

徐海坤,匡邓晖,刘 杰,2,龚春叶,2

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

1 引言

在核反应堆设计、临界安全分析及屏蔽计算中，中子输运问题是数值计算核心。通常采用确定论方法和蒙特卡罗MC(Monte Carlo)方法(又称为随机模拟法)[1,2]进行数值模拟。其中基于组合几何的蒙特卡罗MC方法具有几何仿真能力强、物理建模完善、模拟近似程度小、对问题维度不敏感等优点。但是，利用蒙特卡罗方法进行数值模拟时也有明显缺点，例如，收敛速度慢，模拟时间长[3 - 5]。随着现代超级计算机技术的飞速发展，蒙特卡罗方法天然的粒子并行性使得这些缺点得以明显改善，但在超大规模的核反应堆数值模拟计算中，模拟时间长的问题依然比较突出。

堆用蒙特卡罗分析程序RMC(Reactor Monte Carlo Code)[6,7]是由清华大学工程物理系核能科学与工程管理研究所反应堆工程计算分析实验室(REAL团队)自主研发的,用于反应堆堆芯计算分析的三维输运蒙特卡罗程序。RMC程序针对反应堆计算分析中的基本需求，同时结合新概念反应堆系统设计时几何结构灵活、中子能谱复杂及材料组分多样、各向异性及泄露强(某些特定情况)等特点进行研发，是多物理多尺度耦合核能系统数值分析平台的物理计算核心。

RMC程序有以下几大特点：

(1)MPI+OpenMP 2层并行。RMC程序实现了MPI+OpenMP 2层并行计算模式，上层使用MPI实现对等模式下的并行计算，下层采用OpenMP实现多线程并行计算，用任务分解方法实现粒子并行。同时，MPI也分成2级并行，第1级采用区域分解方法划分几何区域，以减少内存占用量，第2级实现对粒子的并行计算。

(2) 大量采用vector作为基本数据类型。RMC程序使用C++作为开发语言，程序中大量使用容器vector代替数组作为基本数据类型。vector可变容量特性使得编程极为方便，但容量变化的同时需要不停地从内存中申请和释放空间，尤其是在多线程的环境中，内存申请和释放会引起线程之间的资源竞争，从而导致程序性能下降。

(3) 文件I/O量大。RMC程序实现了断点续算功能，续算功能由主进程控制，每隔一段时间主进程负责收集各个进程中粒子的相关信息，然后由主进程输出到相关文件中。当计算规模比较大时，续算文件的大小往往能达到几十GB甚至上百GB，文件读写均需要较长的时间。

RMC程序在实际进行反应堆数值分析计算时，对于千万网格的计算规模，在数千核的计算机上完成单个燃耗步的计算需要数小时之久，完成一次大规模的数值分析需要几天甚至一周时间。在有限的计算资源上如何提高计算的效率是本文研究的主要目的。

2 采用的优化方法

根据RMC程序的特点，本文先后开展了以下工作。

2.1 基于TCMalloc的动态内存分配优化方法

RMC程序在内存管理上有以下几个特点：

(1) 内存分配量大，在进行千万网格数值模拟时，单个进程所占内存在70 GB以上。

(2) 多线程中使用vector容器，内存申请释放比较频繁，线程之间内存锁操作开销较大，且容易产生内存碎片。

在多线程计算环境中，要保证内存数据的访问效率和稳定性，首先要考虑多线程内存管理的有效性，即内存分配释放的快速性、碎片率、利用率和扩展性问题。RMC程序使用vector容器作为基本数据类型，在多线程的计算过程中需要不断地申请和释放内存资源，需要选择高效的内存管理算法来减少多线程环境中的锁竞争问题。

多线程内存管理方式主要采用多堆的组织方式，内存管理算法用于处理堆、处理器和线程之间的对应关系。目前PTMalloc是Linux的主流内存管理算法，该算法为一种基于线程随机挑选内存堆的算法。该算法每次分配内存首先随机查找未被使用的内存堆，然后在该内存堆中分配内存。 释放内存时需将分配的内存堆收回该内存堆。在多线程的环境中，该内存管理算法可能存在多个线程共享一个内存分配区，使各个线程在申请释放内存时都需要加锁、解锁操作，且PTMalloc采用的是比较耗时的互斥锁，这在一定程度上使得PTMalloc在多线程竞争时的性能受到影响。

TCMalloc是Google开发的面向多线程的内存管理库。与PTMalloc相比，TCMalloc具有更高的内存分配效率和内存利用率，提高了多线程并发环境中内存管理的效率。TCMalloc的内存结构为三级缓存模型，分为线程局部缓存区(ThreadCache)、中央缓存区(CentralCache)和全局的页面堆(PageHeap)。线程局部缓存区负责线程内的内存管理，当局部缓存区能满足线程分配内存需求时，线程内无需锁操作即可进行内存分配；当局部缓存区无法满足线程内存分配需求时，则采用锁机制从中央缓存区和全局页面堆分配内存。频繁地内存操作时，线程局部缓存区和中心缓存区之间可能存在周期性的抖动现象。

本文基于TCMalloc的内存管理算法，结合RMC程序的特点，对TCMalloc大内存的管理算法进行了以下几点改进。

2.1.1 预先分配大量内存块

在内存结构上，本文沿用TCMalloc的三级缓存结构，以保证内存分配效率。由于RMC程序占用内存量比较大，因此在初始化管理时，本文根据RMC程序的内存需求，提前设置好中央缓存区中的内存数量、运行时线程数以及每个线程局部内存区的内存数量，预先分配大量内存块到线程局部缓存区和中央缓存区，以减少线程在分配内存时向中央内存区申请内存和释放内存所带来的锁开销。

Figure 1 Loop iteration task queue of openMP schedule strategy图1 OpenMP调度策略下循环迭代的任务队列

2.1.2 一次申请多个大内存块

TCMalloc算法将大小在32 KB以内的内存块称为small内存，将32～128 KB的内存块称为large内存，将超过128 KB的内存块称为huge内存。对于huge内存TCMalloc直接使用mmap映射内存的方法来进行内存申请操作。当线程局部缓存区的内存无法满足huge内存的申请操作时，TCMalloc采用锁机制一次性从全局页面堆申请多个内存堆用于huge内存的分配操作，这样避免了在大内存分配时的多次锁操作。预计一次性需要申请的内存堆大小按式(1)计算：

(1)

其中,n为线程数，α为内存堆中的空间系数，M为全局内存堆的大小，Li为局部缓存区内存块总大小。

2.2 OpenMP线程调度策略

RMC程序中最主要的计算部分为粒子输运计算，即对所有粒子进行追踪(track)、采样(sample)和碰撞(collide)计算。如在使用32进程、每进程创建10个线程计算共50万粒子的算例中，粒子输运计算部分占据了RMC程序约80%计算时间。

在RMC代码中粒子输运部分体现为对粒子循环的for循环结构，由于各粒子模拟之间不存在数据依赖，RMC使用OpenMP对该循环体进行了线程级并行，这同MPI并行一样，是对粒子并行的粗粒度并行结构。由于输运计算内部有众多判断分支，每个粒子输运时间必然都不尽相同，且无法预测，因此需要对该粗粒度的omp for循环体采用合适的并行调度策略，以尽量减少线程级并行计算时间。

OpenMP对循环体有3种调度策略：静态调度(static)、动态调度(dynamic)和指导调度(guided)，并通过子句schedule(type,chunksize)来控制所采取的调度策略。循环体可以看成是各循环迭代块组成的任务队列，循环块大小由chunksize设置。

图1演示了当循环结构迭代900次、由3个线程并行、chunksize=100时，3种调度策略的循环迭代块划分方式，以及迭代块任务队列分配方式。其中的方块代表迭代块，块中数字表示迭代块大小。

值得注意的是，guided调度策略中，队列第k个迭代块大小由式(2) 确定，依次减小，直至减小到chunksize。

(2)

其中，LCk为第k个迭代块之后剩余的迭代次数，NP为线程数，p为比例因子，编译器中通常设置为1。

对于无法预测每个粒子输运时间的随机循环结构，应该采用dynamic或guided动态调度策略，并让迭代块chunksize尽量小，这将使得各线程的计算负载更加平衡。不过static调度迭代块分配方式固定，而dynamic和guided调度是由运行时系统动态决定的，相较static方法，dynamic、guided方法和较小的chunksize会增加系统调度开销。但是，考虑到大规模问题中粒子抽样、追踪花费的计算时间较长，可以忽略增加的调度开销。总而言之，针对蒙特卡罗类粒子输运程序，采用chunksize小的动态调度策略，理论上能得到更短的整体计算时间。

2.3 vector内存对齐优化方法

许多计算机系统对基本类型的数据在内存中存放的位置有限制，它们会要求这些数据的首地址是某个数的倍数，通常是4或者8的倍数，这就是所谓的内存对齐。在多数情况下，经过内存对齐后，CPU的内存访问速度会大大提升，也更容易对数据做向量化之类的数据并行调优。因此，如果要提升应用程序的性能，所有的程序数据都应尽可能地对齐。现在大多数的编译器都已经实现了自动对齐策略，但手动的数据对齐仍然是常用的应用程序性能优化手段。

RMC程序以vector作为基本数据类型，vector是C++容器中最常用的基本数据类型之一。C++中的vector数据类型原型声明如下所示：

template 〈typename _Tp,typename _alloc=std::allocator〈_Tp〉〉

class vector::protected _Vector_Vase〈_Tp,_Alloc〉;

vector底层默认调用系统的std:allocator作为内存分配器，而std::allocator分配器通过Linux系统默认内存分配函数malloc实现，分配后的内存首地址是否对齐取决于编译器。使用自定义的对齐内存分配器(AlignedAllocator)替换vector默认分配器std::allocator，便可强制实现vector内存首地址对齐。不同于std::allocator利用malloc函数，AlignedAllocator利用C语言底层内存分配函数posix_memalign实现对齐内存分配。

其中AlignedAllocator部分定义如下所示：

template 〈typenameT,AlignmentAlign〉

classAlignedAllocator

{

public:

typedefT*pointer;

typedef constT*const_poster;

pointer allocate(size_tn,typenameAlignedAllocator〈void,Align〉::const_pointer=0)

{

const size_talignment=static_cast〈size_t〉(Align);

void *ptr=allocate_aligned_memory(alignment,n*sizeof(T));

if(ptr== nullptr) {

throw std::bad_alloc();

}

returnreinterpret_cast〈pointer〉(ptr);

}

void * allocate_aligned_memory(size_talign,size_tsize)

{

assert(align>=sizeof(void*));

if(size== 0) {

returnnullptr;

}

void*ptr=nullptr;

intrc=posix_memalign(&ptr,align,size);

if(rc!= 0) {

returnnullptr;

}

returnptr;

}

}

2.4 基于HDF5的并行I/O技术

HDF5(Hierarchical Data Format Version 5)是近年来流行的一种针对科学数据进行高效存储和管理的新型数据存储方案,RMC程序同样采用HDF5作为其中间续算文件数据存储方案。

HDF5使用类似UNIX文件系统的层次结构来描述所存储的数据集。如用/foo/zoo来表示一个名为foo的group之下名为zoo的dataset。dataset是HDF5数据的基本类型，主要包括数据本身和描述数据用的元数据。

如图2a所示,HDF5的串行I/O方式通常是按照file→group→dataset→dataspace→data流程依次创建各对象，最终由单个进程将数据读/写进dataset中。HDF5也能够进行并行I/O，它通过底层调用MPI-2并行I/O接口来实现并行I/O，可以由多个进程同时读写文件中同一个dataset的各部分。如图2b所示，编写并行代码需要在串行流程上增加一些操作，首先需要给file和dataset赋予mpio属性，将dataspace划分成归属各子进程的slabspace，各子进程在slabspace基础上，同时读写划分成子数据的subdata。这样便完成了并行I/O操作。

Figure 2 Programming process of HDF5图2 HDF5编程流程

RMC的续算文件通常需要保存几何、材料、粒子、计数器和燃耗等数据，这些数据主要与网格规模和粒子规模有关。在大规模算例中，续算文件通常达到几十、上百GB。例如在320万粒子的千万级网格算例中，文件大小高达148 GB。由于RMC对粒子进行了数据分解、对网格进行了区域分解，从而续算文件所需数据分散在各进程中。而RMC程序通常采用HDF5串行I/O方式访问续算文件，即主进程在进行数据读写之前，必须先从其他进程中收集完整信息，再将信息数据写入/读出文件。而采用HDF5并行I/O方式，主进程无需收集数据，各进程便可直接将子数据写入/读出文件，从而节省MPI通信时间。尤其对于RMC具有几十、上百GB续算文件数据的大规模算例，并行I/O方式能够节省的通信开销非常可观。

此外，相比串行I/O，并行I/O具有更高的读写速度，不过并行I/O性能受到计算机系统存储硬件配置的限制。一般来说，数据存储服务器结点越多，I/O并行度越高，并行I/O的加速优势越能体现出来。

3 测试与分析

本节对以上优化方法开展单项性能测试，测试在高性能计算系统上进行，该系统每个计算结点包括2个Intel Xeon E5-2660 处理器(主频2.6 GHz、10个核，每核包括32 KB L1数据Cache和32 KB L1指令Cache，256 KB L2 Cache,10核共享25 MB L3 Cache)和256 MB内存，结点间使用国防科技大学自主研发的Express高速互连网络，通信带宽25 Gbps。系统采用Lustre并行文件系统，共有16个存储结点。编译环境使用基于gcc 4.8.5的Intel C/C++Compiler 16.0.3，并行环境使用针对Express-2高速网络优化的mpich3。

3.1 动态内存分配优化测试

将优化后的TCMalloc库编译成静态库，在RMC程序编译阶段链接静态库，程序运行期间会自动调用优化后的内存分配方法。对RMC程序进行实际性能测试，测试规模尽可能重现工业算例所需的运行模式，固定进程数为32，每个进程创建10个线程，改变粒子数，分别测试不同粒子数情况下程序的执行时间，结果如图3所示。

Figure 3 Test results using optimized dynamic memory allocation method图3 动态内存分配优化测试结果

在1 200万网格数量的实验测试中，分别测试了200万粒子数、150万粒子数和100万粒子数情况下程序总体的运行时间。在粒子数为200万时，未使用内存优化管理方法、使用TCMalloc和使用优化的TCMalloc 3种情况下程序的运行时间分别为783 s, 683.2 s和671.4 s，使用本文优化的内存管理方法，较使用系统默认的内存管理方法，程序总体性能提升了16.25%，较使用TCMalloc提升了1.7%。

3.2 OpenMP调度策略测试

通过在omp制导指令中调整schedule子句,便可指定OpenMP调度策略，如使用如下制导指令便可将调度策略修改为迭代块大小为100的动态调度策略：

#pragma omp parallel for schedule (dyna- mic,100)

RMC程序中对粒子追踪和随机抽样过程的循环结构采用的原始调度策略是chunksize=nIter/nThread的guided调度，这种调度结果其实和缺省static调度一样，没有达到负载平衡效果。

为了改善这种不合理的调度策略，本文采用了(guided,1),(guided,100),(dynamic,1),(dyna- mic,100)4种较小迭代块的动态调度策略，并同原始调度(guided,nIter/nThread)和缺省static调度进行了对比。表1展示了对于6百万级网格算例，使用32个进程，每个进程创建10个线程，计算200万个粒子，迭代200次时，采用6种不同OpenMP调度策略的RMC计算时间。原始调度策略耗时最长，约252.4 s，采用(guided，1),(guided，100),(dynamic，1)的计算时间最短，约为234.1～238.5 s，整体性能提升约5.82%～7.81%。这表明，对于由粒子模拟存在众多条件语句分支的蒙特卡罗程序，采用chunksize小的OpenMP动态调度策略，能够在一定程度上提升性能。

Table 1 Computation time of particle transportation module under different scheduling strategies

3.3 vector内存对齐优化测试

在RMC程序中，把自定义的内存对齐分配器(AlignAllocator)应用到vector的数据类型定义中的方式如下所示：

vector〈int,AlignedAllocator〈int,64〉〉 Array;

本节对于1 200万网格算例，使用32个进程，每进程创建10个线程，分别测试了100万粒子数和50万粒子数情况下程序的运行时间，表2为测试结果。

Table 2 Test results using vector memory alignment optimization表2 vector内存对齐优化测试结果

从表2中可以看出，使用内存对齐优化方法后，程序总体提升仅有1.8%，说明程序在编译阶段就已尽可能地实现了内存对齐，但手动的内存对齐仍然能使程序性能提高。

3.4 HDF5文件I/O测试

在RMC程序中，HDF5原本使用串行I/O技术，本文采用并行I/O技术后，程序性能得到极大提升。HDF5并行I/O实现方法主要在于，在串行I/O代码的基础上，添加mpio属性，2.4节中阐述了具体方式。本文进行了2组测试：第1组使用16个进程，测试了400万粒子数、600万网格单元规模、72 GB续算文件的算例，测试结果如表3所示；第2组使用了64个进程，测试了320万粒子数、1 200万网格单元规模、146 GB续算文件的算例，测试结果如表4所示。

Table 3 HDF5 I/O test under 6 million grids

Table 4 HDF5 I/O test under 12 million grids

第1组测试中，读写时间分别缩短了90.12%和85.72%；第2组测试中，读写时间分别缩短了96.30%和95.18%。可见HDF5的并行I/O技术极为有效地提升了文件读写性能，而且网格规模越大，性能提升越明显。

3.5 总体性能测试

将以上4种优化方法全部应用到RMC程序中，对RMC程序进行优化后的整体性能测试。测试中同样固定进程数为32，每个进程创建10个线程，分别计算在1 200万网格情况下，200万粒子规模和100万粒子规模的程序总运行时间，测试结果如表5所示。从表5中可以看出，2种粒子规模下，程序的总体性能提升均达到了26.45%以上。

Table 5 Test results of overall performance

4 结束语

本文针对堆用蒙特卡罗分析程序(RMC)，先后开展了一系列的程序性能优化研究，采取了基于TCMalloc动态内存分配优化、OpenMP线程并行调度策略优化、vector内存对齐优化和基于HDF5的并行I/O优化等优化方法，并对优化后的程序进行了实际性能测试。测试结果表明，每种优化方法都使程序性能得到一定程度的提升，其中动态内存分配优化方法使程序性能整体提升了16.25%，表明在RMC程序中，频繁的内存分配是造成程序性能较低的主要原因之一，新的内存管理方法有效提高了多线程环境中内存管理的效率。基于HDF5的并行I/O优化极大地缩减了程序续算文件读写时间。相较串行I/O方式，进程数越多，并行I/O时间越短。因此，使用HDF5的并行I/O非常有效。综上所述，采用一系列的优化方法优化后，程序整体性能提高26.45%以上。