中子输运蒙特卡罗模拟的区域分解方法研究

梁金刚，王 侃，蔡 云,2，孙嘉龙

(1.清华大学 工程物理系，北京 100084；2.中国核动力研究设计院，四川 成都 610041)

近年来，在反应堆物理计算研究领域，求解中子输运问题的蒙特卡罗(简称蒙卡)模拟方法越来越受到重视，正在成为研究的热点和前沿。蒙卡方法作为一种基于概率论的随机模拟方法，能有效借助计算机的速度与精度优势，在模拟粒子输运方面表现出忠实物理过程程度高(如连续点截面)、受几何限制小、对问题维度不敏感等诸多优点，恰好符合反应堆中子输运计算的几何与能谱复杂、多维度、多变量等特点。因而，相较于中子输运方程计算的确定论方法，蒙卡方法被认为是未来反应堆工程计算方法的一种重要选择[1-2]。

然而，将蒙卡方法应用于反应堆全堆芯的大规模计算分析还有一些问题需进一步解决，其中，计算时间长和内存占用大是两个主要问题，著名的“Kord Smith挑战”[3]及Bill Martin所做的修正[1]说明了这一点。针对计算时间长采取并行计算是有效的解决方法。由于蒙卡方法是对大量粒子进行逐一的跟踪模拟，按照粒子数的并行容易实现，已有研究表明，粒子并行方法具有良好的并行效率和可扩展性[4]，可有效解决蒙卡方法计算时间长的限制。而对于内存占用大的问题，目前还没有较成熟的解决方法，区域分解方法被认为是一种可行的思路。

文献[5]最早提出将区域分解应用于隐蒙卡的并行化，并研究了区域分解的实现方法。但区域分解在粒子输运模拟中的应用并未得到推广，在众多著名的蒙卡程序中，鲜有程序具备区域分解功能。近年来，随着蒙卡粒子输运模拟的兴起和反应堆大规模计算中内存问题的凸显，对区域分解的研究逐步增多。其中，美国LLNL对区域分解方法进行了一些研究，其开发的蒙卡程序Mercury[6]具备一定的区域分解功能。在国内，对区域分解的研究属于起步阶段，北京应用物理与计算数学研究所基于JMCT程序进行了区域分解的研究[7]。

1 区域分解方法的基本原理

区域分解方法是指将研究对象从几何上划分为若干(子)区域，对不同的(子)区域进行计算，并通过建立(子)区域间的某种耦合关系使得对整个研究对象的求解结果正确。“分而治之”是区域分解的基本思想。它可将大模型问题化为小模型问题，是解决计算条件不足问题的一种有效方案。

区域分解作为数学方法，最早应用于偏微分方程的求解，其提出目的就是为了应对当时计算机计算条件的限制。求解偏微分方程的区域分解方法如图1所示，一般按照区域划分方式分为不重叠型和重叠型两种。图1中，Ω为计算区域，Γ为区域交界，下标1、2表示区域划分为两个子区域。在计算流程上，以各区域交替计算(如经典的Schwarz交替法)为主要特点，即先计算子区域1的值(假设交界边界面初值)，然后将新的边界值传递至区域2，计算区域2的值，再次更新边界值，重复上述步骤，直到两个区域的解收敛到真实解。

区域分解方法的思想同样可用于粒子输运蒙卡模拟的过程。粒子在不同的子区域运动，当从一个子区域进入到另一子区域时，通过对粒子运动状态进行存储及区域间的传递，实现区域间的耦合计算。一般情况下，不同的子区域由不同的处理器计算，利用蒙卡模拟中粒子的相互独立性，可实现不同子区域的并行计算，因而区域分解方法又称为空间并行方法。

2 基于RMC的区域分解实现及测试

2.1 实现方法

RMC[8]是由清华大学工程物理系核能科学与工程管理研究所反应堆工程计算分析实验室(REAL团队)自主开发的用于反应堆物理分析的三维粒子输运蒙卡程序。本文基于RMC平台进行区域分解功能的研究和开发。

区域分解的实现主要包括两部分内容：区域划分和区域间通信。前者指对模拟区域进行剖分，并将每个区域分配至不同的处理器，在区域划分方法和区域与处理器对应关系上进行探究以实现优化；后者指粒子穿越边界时，区域之间对粒子信息进行存储、传递的通信算法，蒙卡模拟中粒子穿出区域多、通信频繁、信息量大，不同的通信算法对区域分解的性能会有影响。

a——不重叠型；b——重叠型

粒子输运蒙卡模拟的区域分解方法如图2所示，每个区域由一个处理器模拟，当粒子穿出当前所在区域时，处理器立即保存粒子信息(包括位置、方向、速度及随机数种子等)至缓存区，继续模拟下一个粒子，直到当前区域的所有粒子模拟完毕，不同区域开始交换信息，将穿出区域的粒子发送到对应的区域(处理器)，然后进行下一轮模拟，直到所有粒子至历史结束。

本文按照上述方法在RMC中添加区域分解模块，在三维几何下采用曲面组合逻辑的用户输入方式进行不重叠型的区域划分，即用户输入区域边界，根据布尔运算逐个指定子区域范围。这种方式不影响计算模型原始输入，同时具备几何区域划分的灵活性。在区域间通信方面，基于MPI并行库进行粒子信息传递。为避免通信发生死锁，在通信中根据区域编号规定信息发送、接收的顺序。

2.2 测试

为了对实现的区域分解基本功能进行测试，设计了3个算例模型(立方体模型、球模型、组件模型)，如图3所示，按照不同的区域划分方式进行计算，将结果与未进行区域分解的RMC和MCNP的计算结果进行对比，并给出区域分解的并行加速效率。

1) 立方体模型

立方体模型采用均匀材料，对模型进行对称的2区域(分别沿x、y、z3个方向)、4区域和8区域划分，临界计算结果列于表1。计算条件为模拟300代，前100代为非活跃代，每代20 000个中子。测试平台为Windows系统，12核Intel Xeon CPU(型号X5670，2.93 GHz，下同)。

2) 球模型

表2列出球模型的计算结果，其中编号代表不同的划分方式：so0为无区域分解；so1为2区域均分(沿x方向)；so2为4区域均分(x、y4个象限均分)；so3为2区域划分(小球内、外)；so4为2区域划分(大球内、外)；so5为6区域划分(4个小球内、大球内小球外和大球外)。计算条件为模拟300代，前100代为非活跃代，每代50 000个中子。

图2 粒子输运蒙卡模拟的区域分解方法

a——立方体模型；b——球模型；c——组件模型

3) 组件模型

选取重水堆组件进行计算，按照图4所示方式进行2～24区域划分计算，其中2～12区域划分为沿圆柱横切面的角平分面进行划分，24区域划分是在12区域划分的基础上在圆柱轴向中间横切。计算结果列于表3。计算条件为模拟300代，前100代为非活跃代，每代50 000个中子。

2.3 区域分解的结果及并行性能分析

上述3个测试模型分别采用不同的区域划分方式，包括不同的子区域数目、子区域几何类型、子区域划分方向等多种情形。计算结果表明，临界计算有效增殖因数keff与无区域分解的结果均在统计涨落误差内(相差均在2倍标准偏差以内)。

表1 立方体模型的计算结果

表2 球模型的计算结果

图4 组件模型区域划分方式

表3 组件模型的计算结果

在区域分解情况下，粒子随机数种子的分配相对于无区域分解时发生了变化，导致最终结果并非完全一致，但这种分配的变化是随机的，不会影响结果的收敛与准确性(蒙卡方法本身具有统计性)。相对误差总是在2倍标准偏差以内，说明误差是因随机数不同产生的统计性涨落导致的，因此这种结果是可靠的。通过进一步的设计可实现两者结果完全一致，即区域分解结果的可重复性。

3个测试均采取处理器与区域一一对应的方式，区域分解后，处理器并行计算，因而整体缩短了计算时间。对比不同划分方式的加速效果，可以看出，区域划分均衡情况下，具有较好的并行效率，如立方体模型及组件模型(图5，加速比为串行计算与并行计算之比)；反之，区域划分不均衡时，并行效率较差，如球模型。另一方面，随着区域数目增多，区域间通信增多，并行效率将降低。负载均衡与通信性能是影响区域分解算法性能的关键因素。

图5 组件模型对称区域划分下的并行加速比

3 区域分解结果可重复性研究

由于随机数对源中子的分配顺序发生变化，使得粒子运动径迹改变，以致计算结果相对于不分解计算出现偏差。为验证程序的正确性，设计出在输入一致的条件下区域分解的计算结果与无分解结果一致的算法是有必要的。

图6 区域分解下新一代源中子产生顺序的变化

要实现结果的可重复，必须保证所有粒子的运动历史一致，进而给每一个源中子(每一代)分配的初始随机数种子必须一致。在源中子顺序一致的情况下，通过跳跃法预先产生分段的随机数种子，逐一分配给每一个源中子可确保粒子的运动历史不变，这与常见的粒子并行算法结果可重复性使用的方法相同。然而，在区域分解算法中，由于不同区域的粒子可同时运动，穿出区域、尚未死亡的粒子需要等待，被传递后才能继续运动，以至新产生的下一代源中子的顺序被打乱。图6示出区域分解下新一代源中子产生顺序的变化。由图6可见，A中子在运动中，先后在左、右区域各产生1个新中子(分别记为1#、2#)，B中子在右区域产生1个新中子(记为3#)。在无区域分解情况下，新中子的产生顺序为1#、2#、3#。当区域被分为左、右两个子区域进行计算时，A、B中子同时开始运动(分别在两个处理器上)，新一代源中子的产生顺序为1#、3#、2#，亦即2#与3#中子的存放顺序发生了变化。

为保证下一代计算结果不变，需使源中子的随机数分配不变，因而需将顺序变化的源中子重新排序。为此，对于每个新中子，可记录致其产生的旧中子序号和已由此旧中子碰撞产生的新中子数目，据此便可对所有源中子进行排序，使得源中子顺序与无区域分解情况一致。实际计算中，也可根据新源中子的顺序对随机数种子序列进行排序使得两者符合一致。

在实现结果可重复性后，区域分解与不区域分解的程序，两者模拟的粒子历史完全一致，因而程序的统计结果是完全一致的，包括keff、通量等。表4列出按照上述原理对程序修改后，对立方体模型2区域分解下的临界计算得到的keff。计算条件为模拟300代，前100代为非活跃代，每代5 000个中子。由表4可见，2区域与单区域的keff计算结果完全一致。需注意的是，在确保结果一致的计算中增加了排序过程，因而计算时间有所增加。

表4 区域分解实现结果的可重复性

4 小结

对于中子输运的蒙卡模拟方法，区域分解是解除内存限制的值得探索的方向。本文介绍了区域分解的基本思想，分析其实现方法，并基于RMC开发区域分解功能，进行了算例测试和性能分析，总结出影响区域分解并行性能的关键因素，即负载均衡和通信性能。最后研究了区域分解结果可重复性的实现方法，提出通过对源中子排序以保证粒子运动历史的一致性。

参考文献：

[1] MARTIN W R. Advances in Monte Carlo methods for global reactor analysis[C]∥Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering. USA: American Nuclear Society, 2007.

[2] BROWN F B, MARTIN W R. Reactor physics analysis with Monte Carlo[C]∥PHYSOR2010. USA: American Nuclear Society, 2010.

[3] SMITH K. Reactor core methods[C]∥Nuclear Mathematical and Computational Sciences. USA: American Nuclear Society, 2003.

[4] 丘意书，佘顶，范潇，等. 堆用蒙特卡罗程序RMC的全堆计算研究[J]. 核动力工程，2013，34(S1)：1-4.

QIU Yishu, SHE Ding, FAN Xiao, et al. Analysis of full-core calculation of RMC[J]. Nuclear Power Engineering, 2013, 34(S1): 1-4(in Chinese).

[5] URBATSCH T J, EVANS T M. Parallel implicit Monte Carlo in C++[C]∥ISCOPE’98/International Symp. on Computing in Object Oriented Parallel Environments. USA: Los Alamos National Laboratory, 1998.

[6] BRANTLEY P S, DAWSON S A, MCKINLEY M S, et al. Recent advances in the Mercury Monte Carlo particle transport code[C]∥International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2013). USA: American Nuclear Society, 2013.

[7] LI Gang, ZHANG Baoyin, DENG Li. Domain decomposition of combinatorial geometry Monte Carlo transport code JMCT[J]. ANS Transactions, 2013, 19: 1 425-1 427.

[8] WANG Kan, LI Zeguang, SHE Ding, et al. RMC: A Monte Carlo code for reactor physics analysis[C]∥Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo (SNA+MC). France: French Alternative Energies and Atomic Energy Commission, 2013.