基于CPU-GPU异构并行的MOC中子输运计算并行效率优化研究

宋佩涛，张志俭，梁 亮，张 乾，赵 强

(哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001)

特征线方法(MOC)具有强大的几何处理能力和天然的并行特性，被广泛应用于反应堆物理计算。随着计算机计算能力的快速发展，越来越多的反应堆物理计算程序采用MOC作为其主要的输运计算工具，如DeCART[1]、nTRACER[2]、MPACT[3]、OpenMOC[4]和NECP-X[5]等程序。但即使采用基于CPU的高性能计算集群，堆芯规模的MOC输运计算仍然耗费大量的计算时间和计算资源。

近年来，基于GPU和CPU的异构系统得到了快速发展[6]。由于异构系统强大的浮点运算能力，借助英伟达(NVIDIA)推出的统一计算设备构架(CUDA)，CPU-GPU异构系统被广泛应用于科学计算领域。同时，在MOC输运计算领域，已开展了采用GPU加速MOC计算[7-11]及小规模的多GPU并行计算[12-13]相关研究。然而，目前基于CPU-GPU异构系统，针对MOC的大规模CPU-GPU协同计算研究较少。

本文采用性能分析模型，分析二维MOC异构并行算法的并行效率和性能瓶颈。同时针对显著影响异构并行效率的因素，采用输运计算与数据传递相互掩盖的技术，提升异构并行算法的整体性能和并行效率。

1 特征线方法

MOC是求解中子输运方程的一种数值方法，其求解思路是采用坐标变换，将偏微分方程转化为常微分方程。基于平源近似，省略能群标识，特征线形式的稳态多群中子输运方程为：

(1)

其中：φ(s)为中子角通量密度；s为中子穿行方向；Σt为总截面；Q为源项。

(2)

(3)

实际计算中，若干平行射线沿角度Ω穿过平源区n，引入角度离散p，采用数值积分，可求得平源区n的平均中子标通量密度φ为：

(4)

其中：ωp为求积权重；dp为Ωp方向的射线间距。

2 性能分析和效率优化

2.1 异构并行算法

结合目前CPU-GPU异构并行系统的特点，假设异构系统由N个异构节点组成，每个节点包含Nc个CPU核心及Ng个GPU，且通常Nc大于Ng。为充分利用异构系统的计算资源，采用MPI+OpenMP/CUDA的混合编程模型来实现二维MOC异构并行计算。基于该异构模型，采用空间区域分解和MPI规范，实现粗粒度并行；而计算耗时的输运扫描过程(细粒度并行)，则通过基于OpenMP/CUDA的特征线并行来实现。

图1示出基于空间区域分解和特征线并行的MPI+OpenMP/CUDA混合编程模型。假设单异构节点包含2个GPU(Ng=2)及14个CPU核心(Nc=14)。在单异构节点内，待求解区域被分成4个子区，对应启动Ng+Nx(Nx为采用MPI+OpenMP模型计算的进程数，Nx=2)个MPI进程，将每个子区分配给各进程进行计算。其中Ng个进程由Ng个CPU核心启动，每个进程通过CUDA环境驱动1个GPU，采用MPI+CUDA模型进行计算；另外Nx个进程调用剩余的CPU核心，启动Nc-Ng(12)个OpenMP线程，采用MPI+OpenMP模型进行计算，如图1所示，可为每个进程分配6个CPU核心参与计算。

在MPI+CUDA模型中，采用GPU线程并行实现细粒度并行。根据文献[14]研究结果，为获得最优并行性能，单个GPU线程执行单条特征线、单个能群的相关计算。如图1所示，GPU线程t执行特征线m第1个能群的相关计算，线程t+1则执行特征线m第2个能群的相关计算，依此类推。而对于MPI+OpenMP模型，为减少线程开辟和销毁带来的额外开销，单次外迭代对应进行1次OpenMP线程的开辟和销毁。相应的计算任务显式地分配给各线程，且单个线程负责特征线的所有能群和极角的输运计算。如图1所示，CPU线程s首先执行特征线n的相关计算，计算完成后，特征线n+线程数(线程数指的是当前进程中参与计算的OpenMP线程数)对应的计算被分配给线程n，依此类推。

由特征线并行的实现方式可知，本文采用了Jacobi格式的输运扫描算法，减弱了能群间的耦合关系，相对于Gauss-Seidel格式的输运扫描算法，会增加输运计算外迭代次数。对于二维C5G7 1/4堆芯基准题，外迭代次数增量约为12.8%。但文献[14]研究表明，GPU并行计算需通过设置大量并行任务来掩盖GPU访存延迟，而同时将特征线和能群并行，能够保证足够多的并行任务，获得最优的计算性能。因此基于文献[14]的研究结果，在CPU-GPU异构并行计算中依然采用Jacobi格式的输运扫描算法。

图1 基于空间区域分解和特征线并行的混合MPI+OpenMP/CUDA编程模型Fig.1 Hybrid MPI+OpenMP/CUDA programming model with spatial domain decomposition and ray parallelization

为保证异构并行计算中CPU与GPU计算之间的负载平衡，本文采用文献[15]所提出的动态任务分配模型，根据硬件的实际计算性能实现最优任务分配。在该模型中，总计算任务被任意初始分配给MPI进程；随后执行预处理过程，进行若干次输运扫描，并获取各进程输运计算时间；假设每个MPI进程计算能力在整个计算过程中保持为常数，则进程的输运计算时间与所分配的任务量呈正比，可采用预处理过程中初始分配的计算任务和获取的计算时间，求解最优任务分配方案，保证各进程之间的负载平衡。

2.2 计算流程

图2示出二维MOC异构并行计算流程：1) 并行初始化过程中执行空间区域分解，程序将计算任务分配到p个MPI进程；2) CPU端执行本地初始化，生成截面、特征线等信息；3) 对于采用MPI+CUDA并行的进程(进程#1)，将计算所需的信息(截面、特征线信息等)从CPU一次性拷贝到GPU，随后在GPU端执行源项计算、输运扫描和反应率及通量残差计算，再将反应率、通量残差及边界角通量密度从GPU通过CPU拷贝到系统内存；4) 对于采用MPI+OpenMP并行的进程(进程#2)，则依次在CPU进行源项计算、输运扫描和反应率及残差计算，其中输运扫描过程采用OpenMP并行实现；5) 空间子区边界角通量密度交换由MPI通信实现；6) 在主进程计算keff并进行收敛判断。

图2 二维MOC异构并行计算流程Fig.2 Flowchart of heterogeneous parallelization of 2D MOC

2.3 性能分析模型

结合图2所示的计算流程可知，对于异构并行计算，MPI+CUDA模式是影响异构并行性能的主要方面。因此，识别影响MPI+CUDA并行效率的因素并进行相应的优化是提升异构并行效率的主要途径。而影响MPI+CUDA并行效率的主要因素包括：1) 由于引入空间区域分解，边界角通量密度更新滞后，导致迭代次数增加；2) 负载不均衡导致的性能下降；3) CPU和GPU间的数据拷贝过程对性能的影响；4) MPI通信对并行效率的影响。本文采用4个指标量化上述因素对并行效率的影响。程序计算时间由式(5)描述。

Ttotal(p)=T(p)+Tdata_copy(p)+TMPI_comm(p)

(5)

其中：p为进程编号；Ttotal为各进程整体计算时间；Tdata_copy为CPU与GPU间数据拷贝时间；TMPI_comm为MPI通信时间；T为各进程除数据传递(MPI通信和CPU与GPU间数据拷贝)之外的计算时间。

则上述4个影响因素可由式(6)～(9)量化：

(6)

(7)

(8)

(9)

其中：Np为并行进程数；Tserial为串行计算时间。

因此，总并行效率则为上述4个指标的乘积，即：

(10)

2.4 输运计算和数据传递相互掩盖

考虑到图2所示的数据传递(CPU与GPU间的数据拷贝及MPI通信)显著影响程序并行效率，本文采用CUDA流、输运计算和MPI通信异步执行等方案来实现输运计算掩盖数据传递。对于MPI+CUDA并行，本文采用CUDA流技术，将CPU与GPU间的数据拷贝和输运计算分成若干部分，各对应部分组成CUDA流，将“分块”的数据拷贝和输运计算重叠执行，如图3所示。同时，考虑到CUDA函数为异步执行，即CPU端调用CUDA函数后，该函数在GPU端执行，随后计算流程控制权返回CPU。因此在GPU计算期间，与之绑定的CPU核心可执行CPU与GPU间的数据拷贝及MPI通信，实现数据传递与GPU计算的相互掩盖。对于MPI+OpenMP并行，则从原来执行输运计算的t个线程中分出1个线程，专门执行MPI通信，另外的t-1个线程依然执行输运计算，实现MPI通信与输运计算的相互掩盖。相比原执行方式，实现输运计算掩盖数据传递之后，可提高计算效率，如图3所示标注了计算时间节省。

图3 输运计算掩盖数据传递示意图Fig.3 Scheme of transport sweep and data movement overlapping

3 数值验证

3.1 精度验证

本文采用二维C5G7基准题[16]验证异构并行程序计算精度，图4示出该基准题的1/4堆芯布置，每个燃料组件为17×17栅元排列。整个计算区域被划分为51×51个栅元，每个栅元采用如图5所示的网格划分。特征线初始间距为0.03 cm，每个象限选择14个方位角和3个极角，使用Tabuchi-Yamamoto极角求积组。keff和通量收敛准则分别设置为10-6和10-5。计算平台采用天津超算中心TianHe-ES异构计算系统，单个节点硬件参数列于表1。

本征值和归一化功率分布误差结果列于表2，参考解来源于MCNP。表2中：MAX Error为最大棒功率相对误差；AVG Error为平均棒功率相对误差；RMS Error为棒功率均方根误差；MRE Error为平均相对误差。由表2可见，异构并行计算得到的本征值和功率分布误差在可接受范围内，除最小栅元功率外，其他结果相对误差均在MCNP统计不确定度范围内。结果表明，本文所提出的并行算法具备良好的计算精度。

图4 二维C5G7基准题1/4堆芯布置Fig.4 Quarter core configuration of 2D C5G7 benchmark problem

图5 栅元组成及网格划分Fig.5 Fuel pin layout and mesh

项目参数规格CPU2×Intel XeonCPU E5-1690 v4,主频2.60 GHz,14核GPU2×Tesla K80,即4×GK 210,主频824 MHz,计算能力3.7操作系统RedHat 7.2编译环境Intel MPI5.1.3,CUDA 9.2

3.2 异构并行算法效率分析

为分析异构并行算法的计算效率和性能瓶颈，将二维C5G7基准题进行扩展，构造10×10燃料组件的全堆芯问题，堆芯布置如图6所示。采用本文提出的性能分析模型测试了MPI+CUDA并行算法在采用输运计算掩盖数据传递前的强并行效率及算法各部分对并行效率的影响，如图7所示。由图7可知：1) 随并行进程数的增多，MPI+CUDA并行算法强并行效率逐步降低，当采用20进程并行时，算法强并行效率约为78%；2) MPI通信(EMPI_comm)是造成MPI+CUDA并行算法并行效率低下的主要因素，而CPU和GPU间的数据拷贝(Edata_copy)是影响并行效率的另一因素；3) 相对而言，迭代次数增加(Eincr_iter)对整体并行效率的影响较小，且由于图7所示的并行算例基本保持负载均衡，因此负载不均衡(Eload_imb)对并行效率的影响相对较小。图8示出MPI通信和CPU与GPU间的数据拷贝在MPI+CUDA并行中的时间占比，在20个进程并行时，这两部分时间占比分别达到了11%和5%。可见，提升异构并行算法并行效率的主要途径为减小数据传递过程(MPI通信和CPU与GPU间的数据拷贝)对并行效率的影响。

表2 本征值及功率分布的结果Table 2 Result of Eigenvalue and power distribution

图6 虚拟堆芯组件布置Fig.6 Arrangement of mock core assembly

图7 MPI+CUDA并行算法强并行效率Fig.7 Strong scaling efficiency for MPI+CUDA parallel algorithm

图8 MPI+CUDA并行中的时间占比Fig.8 Proportion of execution time in MPI+CUDA parallelization

为提升异构并行算法的并行效率，本文实现了输运计算与数据传递相互掩盖。针对图6所示的全堆芯问题，采用1～5个异构节点，测试了异构并行算法采用输运计算掩盖数据传递前后的效率对比。表3列出程序计算收敛所需的迭代次数。图9示出程序整体计算效率。由表3可知，由于MPI通信和输运计算同时执行，造成迭代格式退化，迭代次数增多，迭代次数增量在3%左右。但由于数据传递时间在异构并行计算中占有显著比例，由图9可知，实现输运计算与数据传递相互掩盖后，程序整体效率有所提升。采用5节点计算时，程序整体效率提升达到8%。

表3 输运计算掩盖数据传递前后迭代次数对比Table 3 Comparison of iteration number with and without transport sweep and data movement overlapping

图9 程序整体计算效率Fig.9 Overall calculating efficiency of program

图10示出MPI+OpenMP/CUDA异构并行模式下的以计算节点为参考的强并行效率(单节点计算效率为100%)。由图10可知，未采用输运计算掩盖数据传递时，强并行效率随节点数的增加下降较快，这是由于在这种情况下，MPI通信和CPU与GPU间的数据拷贝会显著影响程序的整体并行效率。而采用输运计算掩盖数据传递后，强并行效率有明显提升，在5节点并行时，强并行效率达到95%。

针对图6所示的问题，采用多CPU节点进行并行计算，其中CPU节点与异构节点除GPU以外的部分配置相同。图11示出CPU-GPU异构节点和CPU节点的性能对比。其中，异构并行计算基于MPI+OpenMP/CUDA(输运计算掩盖数据传递)并行模式完成；基于MPI的并行计算则分别采用5、10、15、20个CPU节点并行执行。由图11可知：CPU-GPU异构节点计算性能明显优于CPU节点，4个CPU-GPU异构节点计算能力优于20个CPU节点计算能力。

图10 MPI+OpenMP/CUDA强并行效率对比Fig.10 Strong scaling efficiency comparison of MPI+OpenMP/CUDA

图11 CPU-GPU异构节点与CPU节点计算效率对比Fig.11 Calculating efficiency comparison of CPU-GPU heterogeneous node and CPU node

4 结论

本文基于空间区域分解和特征线并行，采用MPI+OpenMP/CUDA编程模型，实现了适用于CPU-GPU异构系统的二维MOC异构并行算法。同时，采用性能分析模型，针对影响异构并行算法并行效率的影响因素进行了分析，并针对性地进行了优化，实现了输运计算和数据传递的相互掩盖，缓解了MPI通信和CPU与GPU间的数据传递对程序整体计算效率的影响。数值结果表明：程序具备良好的计算精度；MPI+OpenMP/CUDA异构并行模式下，数据传递(MPI通信和CPU与GPU间的数据拷贝)是影响并行效率的主要因素；实现输运计算与数据传递相互掩盖能有效提升程序整体性能；5节点(包含20个GPU)并行时，程序整体效率提升达到8%，强并行效率达到95%。同时，CPU-GPU异构节点性能明显优于CPU节点，4异构节点计算能力优于20个CPU节点。未来工作将侧重于基于大规模CPU-GPU异构系统的粗网有限差分(CMFD)加速算法研究。