异构系统三维高保真堆芯中子输运计算程序ALPHA研发进展

张 乾，梁 亮，宋佩涛，李 颂

(1.哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001；2.浙江大学 物理学系 浙江近代物理中心 先进核能理论与应用实验室，浙江 杭州 310030；3.西安核创能源科技有限公司，陕西 西安 710077；4.中国辐射防护研究院，山西 太原 030006；5.海军工程大学 核科学技术学院，湖北 武汉 430033)

近年来，随着数值反应堆概念的提出，三维全堆芯高保真中子输运计算成为了业界研究的热点。其中特征线方法(MOC)由于其灵活的几何处理能力以及天然的并行特性，被广泛应用于数值反应堆项目中的中子输运求解，如NECP-X[1]、MPACT[2]、DeCART[3]、PROTEUS[4]以及KYCORE[5]等。近年来，随着人工智能技术和计算机技术的迅速发展，装载图形处理单元(GPU)的异构计算系统具备强大的计算能力以及效能方面的优势，在工程计算领域得到了广泛应用。异构系统为全堆芯高保真中子输运计算的进一步工程化应用提供了新的思路，韩国首尔大学的初步研究表明，全堆芯中子输运程序nTRACER在移植到异构系统后，在单机环境下即可取得与小型传统集群相等的效率[6]。在工业界，GPU输运加速已被应用于上海核星核电科技有限公司的燃料组件计算程序ROBIN3中。另一方面，除中子输运外，韩国首尔大学实现了基于GPU大规模并行的燃耗计算[7]和核热耦合[8]；哈尔滨工程大学在ALPHA基础上初步实现了基于GPU并行的二维瞬态MOC程序开发[9]。可见，在GPU硬件下具备实现数值反应堆全面应用的潜力。

哈尔滨工程大学核动力仿真研究中心开发基于异构系统的三维全堆芯高保真中子输运计算程序ALPHA，初步实现异构系统上的全堆芯高保真中子输运计算。ALPHA程序以二维特征线中子输运GPU并行算法为核心求解器，基于性能分析模型对特征线扫描算法进行了深度优化；采用MPI+CUDA的混合编程模型，实现空间区域分解的粗粒度并行和特征线细粒度并行，并针对性地提出通信掩盖和异步数据拷贝优化方案来提升算法并行效率；在共振处理中，提出了细群-子群二级离散策略并采用多群求解核心适配异构系统；基于MOC-EX方法，实现了三维全堆芯中子输运异构并行计算，并在异构系统中实现了粗网有限差分加速。

1 方法模型

1.1 二维特征线GPU并行算法及性能优化

在充分了解异构系统的组成、GPU硬件特性和CUDA编程模型的基础上，将特征线方法最为耗时、并行度最高的输运扫描过程采用GPU的CUDA核心进行并行求解，将流程控制与逻辑处理部分由CPU完成，设计出二维特征线的GPU并行算法，如图1所示，其中，白色框图为CPU端执行的算法，橙色框图为CPU与GPU之间的数据传递，绿色框图为CUDA核心执行的操作。由图1可知，GPU主要负责源项计算核函数和输运扫描核函数的计算，其中，输运扫描核函数是特征线计算的主要计算消耗。输运扫描包括不同方位角下特征线、能群、极角的多层循环，可采用特征线并行、特征线-单能群并行、特征线-单能群-单极角并行3种并行模式；输运扫描核函数包括步特征线-指数函数、步特征线-线性插值及菱形差分；浮点运算精度可采用单精度、双精度及混合精度。针对这些技术选项，进行不同规模问题的精度与效率测试，提出了适用于GPU硬件的高效二维MOC中子输运并行方案。

1.2 MPI+CUDA混合编程与通信策略

分别基于空间区域分解和特征线并行，采用MPI和CUDA编程模型，粗粒度并行由空间区域分解结合MPI规范来提供，细粒度并行则采用CUDA编程模型调用GPU实现，计算框架如图2所示。

图2为MPI和CUDA混合并行框架显示，影响并行效率的主要因素包括：1)空间区域分解导致的迭代格式退化及迭代次数增加；2)计算子区间负载不均衡导致并行效率低下；3)主机内存与设备显存之间的数据传递，与输运计算存在数据依赖。其中，MPI通信及主机内存和设备显存之间的数据拷贝(数据传递)是影响异构并行效率的关键因素。因此，ALPHA提出并采用MPI通信掩盖和CUDA流技术，将数据传递与输运计算重叠执行，隐藏大部分数据传递时间，提升异构并行算法并行效率。采用CUDA流操作，将主机向设备的数据传递、GPU扫描计算及设备向主机的拷贝分割，将各部分对应的操作分别压入不同的CUDA流，从而分块实现输运计算与通信数据拷贝，可采用输运计算掩盖部分数据传递，减少数据传递过程执行时间，提高并行效率，如图3所示。

图2 基于空间区域分解和特征线并行的异构并行框架

图3 输运计算掩盖数据传递示意图

1.3 细群-子群二级离散策略与多群求解核心

超细群方法的基本思路为使用精细能群结构处理共振效应，再通过求解慢化方程归并能群。基于此思路，结合超细群方法和子群方法各自的原理及特点，可将超细群精度的优势和子群相对较高的计算效率相结合。首先对传统多群结构采取进一步的能量离散优化，再通过子群方法对优化后的共振群进行截面离散，最后求解共振群慢化通量，以归并有效共振截面用于输运计算。优化后的多群结构称为细群，在细群的基础上进行子群计算，形成细群-子群二级离散策略[10]。

细群-子群二级离散所采取的细群结构以SHEM-361能群结构为基础，并调整共振能量范围，扩充共振群的数量。在描述各核素共振峰的分布情况条件下，并使细群能量边界与常见多群数据库共振能群边界重合，将整体共振能量范围选取为1.855 4～9 118 eV，经优化后的细群结构在此范围内为289群。细群结构的快群和热群仍采用原SHEM-361能群结构，分别为56群和63群，因此细群总群数为408群。由于408群的细群结构并不适用于多群输运计算，与超细群方法相似，细群-子群计算同样对共振群进行了并群处理。根据式(1)建立逐群求解慢化通量的递推表达式，方程源项的其他能群变量均来自上游共振群或快群：

(1)

式中：Ω为中子运动方向；r为中子空间位置；g为能群编号；Σt为总截面；φ为中子角通量密度；Σt,g,l为核素l第g细群的总截面；Σs,g′→g,l为第g′群散射至第g群的散射截面。与超细群方法类似，细群在快群能量段同样采用了渐进能谱1/E，首个共振群的源项可通过快群的渐进能谱和散射矩阵求得。

细群-子群方法对每个共振核素采用8个等效子群进行固定源方程求解，子群输运计算将对计算系统内的每个共振核素的每个等效子群进行循环。对于Mr个共振核素，子群固定源方程对MOC输运求解器的调用次数为Mr×8次，每次调用均需要进行特征线扫描计算。基于MOC方法多群输运计算，可将Mr×8个子群合并，组成一个多群核心[11]。将该多群结构应用到子群固定源输运计算中，特征线扫描将在整个多群结构同时进行，实现对所有子群的固定源通量的同时求解。多群计算把对子群的循环转移至整个计算循环的最内侧，一次全局模块化特征线扫描即可对全部子群进行通量计算。整个多群结构同时进行通量收敛，因此全局迭代次数与逐群迭代次数最多的子群相同。相比于逐群计算频繁地进行特征线扫描，多群计算可显著减少输运计算量。两种循环逻辑的对比如图4所示，其中，r为共振核素编号；Mr为共振核素总数；i、I为子群编号和子群数目；h、H为子群分组编号和总组数；NG为分组内子群数目。数值结果[11]表明，多群核心在GPU硬件上的效率更优。

a——逐群计算；b——多群求解核心

1.4 异构系统三维MOC-EX与粗网有限差分加速

ALPHA的三维特征线输运计算策略采用中子角通量密度轴向变量正交展开方法MOC-EX[12]，核心思想是将中子角通量密度轴向变量采用正交多项式展开，通过多项式正交特性，推导得到耦合的二维中子输运方程，最终求解一系列耦合的二维特征线方程，得到三维中子通量密度分布和本征值。将三维问题在轴向划分为若干层，如图5所示。

图5 三维MOC-EX分层求解示意图

(2)

图6 CMFD加速三维MOC-EX中子输运求解流程

采用源迭代法求解多群CMFD方程时，核心工作为线性方程组的求解，其中的系数矩阵为大型稀疏矩阵。迭代法是求解线性方程组最基本的数值方法，ALPHA采用GPU线程并行，实现了超松弛迭代(SOR)方法的CMFD并行求解。对线性方程组中的系数矩阵和向量进行红-黑网格分组后，直接对红色网格对应的分量进行并行求解，随后求解黑色网格分量时，根据相邻红-黑网格之间的相关性，采用已求得的红色网格分量进行求解。

2 数值结果

2.1 特征线GPU并行加速优化分析

基于C5G7二维基准题[13]，针对3种输运扫描核函数，采用NVIDIA官方提供的性能分析工具集，从指令执行数、DRAM效率、线程束签发效率、缓存命中率等多个分析指标的角度进行优化，具体优化细节参考文献[14]。分别统计双精度和单精度浮点运算情况下CPU和GPU的计算时间，表1为不同类型核函数及不同浮点运算精度情况对应的程序整体计算时间，同时给出了单GPU对CPU单核心的加速效果。由表1可知：双精度浮点运算情况下，单GPU对CPU单核心的加速效果达到35倍以上，3种输运扫描核函数的加速效果接近；采用单精度浮点运算时，单GPU对CPU单核心的加速效果达到100倍以上，其中，基于步特征线的指数计算方式加速效果最佳；在CPU硬件上，基于菱形差分和步特征线-线性插值的计算方式具有较高的效率，而在单精度GPU计算条件下，一般认为较为耗时的步特征线-指数计算反而取得较高的效率。这表明针对GPU硬件对中子输运计算方法进行适配和改进是有必要的。ALPHA的底层输运GPU并行求解核心，经过优化后展示出较高的效率，为ALPHA的整体效率表现奠定了基础。

表1 C5G7二维基准题计算时间

2.2 异构并行效率分析

将C5G7二维基准题扩展形成二维全堆芯算例，燃料组件采用10×10排列，外围设置了厚度为21.42 cm的反射层，在天河一号异构节点上进行MPI+CUDA混合编程模型的效率分析。应用通信掩盖和异步数据拷贝的优化方案后[15-16]，计算效率和并行效率的提升效果示于图7、8。数值结果表明，采用通信掩盖和异步数据拷贝后，大部分数据传递时间通过输运计算进行掩盖，程序执行时间有所降低。采用5节点(20个GK210 GPU)并行计算时，程序整体效率提升达到8%。对于异构并行，由于MPI通信和主机内存与设备显存间的数据拷贝显著影响整体并行效率，强并行效率随节点数的增加下降较快。而采用通信掩盖和异步数据拷贝后，强并行效率有明显提升，5异构节点并行计算对应的强并行效率达到95%。

图7 二维全堆芯问题采用通信掩盖和异步数据拷贝前后异构并行整体计算效率对比

图8 二维全堆芯问题MPI+CUDA/OpenMP异构并行强并行效率对比

2.3 共振计算精度验证与效率分析

本文计算所选取的三维问题为VERA基准题[17]。该基准题由Westinghouse典型单组件在轴向排布而成，其几何结构如图9所示。

a——常规组件；b——含Pyrex棒组件；c——轴向结构

图10为3A和3B算例的轴向归一化功率分布计算结果。ALPHA程序可准确计算三维轴向功率分布，对绝大多数燃料层的功率相对计算偏差均在1%以内。表2为径向积分功率和有效增殖因数相对计算偏差，两个算例整体均具有较高的计算精度。使用4张GK210 GPU并行计算，MOC多群输运计算模块采用CMFD进行加速。表3为VERA基准题两个算例的共振计算和总计算效率。共振计算时间占总时间比例为10%左右。

表2 VERA三维基准题计算结果

表3 VERA三维基准题计算效率

a——3A算例；b——3B算例

2.4 MOC-EX在全堆芯问题中的验证与效率分析

为验证ALPHA对于小型压水堆问题的计算能力，本文首先选取KAIST三维堆芯问题[18]进行宏观输运问题的验证与效率分析。为验证程序计算精度，采用MIT开发的Monte Carlo程序OpenMC[19]对KAIST三维基准题进行模拟，获取基准解。采用ALPHA程序计算KAIST三维问题提棒与插棒工况，keff和归一化棒功率分布相对偏差列于表4。插棒工况下的堆芯归一化棒功率及相对偏差(绝对值)分布如图11所示。KAIST三维基准题数值结果表明：ALPHA程序计算结果能保证良好的计算精度，对于提棒和插棒工况，keff和棒功率分布均与Monte Carlo参考解吻合较好。

表4 KAIST三维基准题keff和中子通量密度分布结果

a——归一化棒功率分布；b——棒功率相对偏差分布

表5为ALPHA程序采用不同硬件的计算效率。相对于CPU计算，GPU加速效果明显，表格后三列给出了GPU相对于CPU的加速比。1080Ti GPU相对于i9-7900 CPU单核心，整体加速比为37.5。对于KAIST提棒问题，计算时间从CPU的6.7 h降低到GPU的10 min。GPU在加速MOC和CMFD计算中均表现出良好的加速效果。MOC和CMFD计算的1080Ti GPU加速比分别为42.5和15.2。相对于CPU计算，ALPHA的GPU加速效果明显。

表5 KAIST三维基准题不同硬件计算时间对比

为测试ALPHA在更接近实际堆芯精细计算需求的压水堆堆芯的表现，基于C5G7三维基准题，构造了典型压水堆1/8堆芯47群算例，如图12所示，宏观截面来源于OpenMC。

图12 压水堆1/8堆芯算例径向组件布置

燃料栅元径向网格划分采用3/2×4形式(燃料3环、慢化剂2环、辐向4个扇形分区)，径向慢化剂区域采用0.252 cm×0.252 cm的矩形网格划分。堆芯轴向活性区平均划分为40层，顶部反射层平均划分为15层。射线间距选择0.05 cm，一个卦限内选取8方位角和2极角，采用Tabuchi-Yamamoto极角求积组，keff和中子通量密度收敛准则分别为10-6和5×10-4。分别采用5异构节点(20个GK210 GPU)和121个E5-2690 CPU核心并行求解，keff及计算时间列于表6。可知，对于压水堆1/8堆芯47群算例以及本文所采用的计算硬件，CPUs/GPUs异构并行整体计算效率达到CPU并行的10倍。

表6 压水堆1/8堆芯47群算例计算效率及存储空间占用量

3 结论

基于异构并行思路，哈尔滨工程大学核动力仿真技术研究中心开发了异构系统中子输运计算程序ALPHA，不同于传统的数值反应堆应用，为了适配异构计算平台，ALPHA程序在二维特征线并行、共振计算、异构并行、三维特征线输运及加速等方面均采取了特殊的计算策略和细致的优化。当前数值结果表明，ALPHA程序已具备在单机或异构超算开展初步数值反应堆中子学计算的能力，在保证精度的前提下，具备较高的计算效率和良好的可扩展性，并探索了数值反应堆轻量化、工程化应用的可行性。