基于权窗的全局减方差方法中“长历史”问题的研究与改进

2017-12-19张恒权李佳祝庆军吴玉潇聂星辰

核技术 2017年12期

张恒权李佳祝庆军吴玉潇聂星辰

张恒权1李佳1祝庆军2吴玉潇3聂星辰1

1（中国科学技术大学核科学技术学院合肥 230027）2（中国科学院等离子体物理研究所合肥 230031）3（广西防城港核电有限公司防城港 538000）

在具有全局特性的蒙特卡罗粒子输运计算中，常采用基于权窗的全局减方差方法，由此引起的“长历史”问题使计算效率大大降低。为了提高蒙特卡罗计算效率，减小计算资源的浪费，本文对“长历史”的成因及解决方案进行了研究。通过研究中国聚变工程试验反应堆(Chinese Fusion Engineering Testing Reactor, CFETR)中子输运计算中“长历史”的事件记录文件，提出了“长历史”的成因与粒子连续在轮盘赌中存活有关，针对此成因提出了通过粒子在轮盘赌中存活的次数来限制粒子分裂的解决方法，并通过CFETR模型对这个方法进行了测试。结果显示，次数限取1−4时，并行效率达到80%，无“长历史”出现，品质因子(pFoMG)达到不采用次数限时的约3倍。综上，采用这种方法可以有效克服“长历史”问题，使全局减方差方法的计算效率明显提升。

全局减方差，全局权窗，长历史，轮盘赌，次数限

在一些大型装置，如中国聚变工程试验反应堆(Chinese Fusion Engineering Testing Reactor, CFETR)的蒙特卡罗输运计算中，常需要求解全空间中子或光子通量分布及停堆计量分布。为缩短得到精度较高的全局解的计算耗时，在计算中常采用两种加速方法：并行计算和减方差。通过并行计算的方法，用多个处理器同时进行粒子输运计算，从而在相同时间内计算更多的粒子数(Number of Particles, NPS)[1]，减小结果的相对误差。为了将粒子输运到全空间中，经常采用一种基于权窗的全局减方差方法(Global Variance Reduction, GVR)[2−6]。GVR方法主要使用粒子在相空间中的分布信息生成全局权窗，再利用全局权窗进行输运计算，多次迭代后便能得到精度较高的全局解。

用蒙特卡罗软件（如MCNP[1]）进行中子输运计算并使用GVR这一减方差方法时，经常会遇到“长历史”问题。计算一个“长历史”的用时远大于其他正常历史的用时，严重降低了并行效率，浪费了计算资源，使计算效率大大降低。关于“长历史”的成因，文献[7−8]认为是装置中一些细长的真空管道导致的。Turner[7]就此提出了一种动态调整权窗的解决方法。本文研究了CFETR模型中子输运计算中的“长历史”现象，通过观察“长历史”事件日志(Event Log)[1]，得出了“长历史”的另一个成因，由此提出了一种解决方法。通过测试，发现这个方法使计算效率有了很大的提升，解决了“长历史”问题。同时发现这个方法对粒子正常输运影响很小，也不会产生过度偏倚导致的对统计结果方差的影响。

1 GVR方法中“长历史”现象与成因分析

GVR方法使用中子在相空间的分布信息生成全局权窗，这个信息可以是通量、能量、径迹数等[2]，本文采用通量这一信息生成全局权窗。通常有两种方法估计中子通量的分布：一种是采用离散坐标的技术，如FW-CADIS4；另一种称为MAGIC (Method of Automatic Generation of Importances by Calculation)[4]，利用MCNP网格计数估计中子通量分布。评定MAGIC GVR方法计算效率的参数主要有并行效率(Parallel Efficiency, PE)、FoMG、 pFoMG、、Scoring[2−5,7]。其中：FoMG、pFoMG是全局品质因子，pFoMG的定义考虑到并行计算效率，是FoMG与PE的乘积；是网格计数相对误差的标准差，反映相对误差在全空间分布的不一致性；Scoring是网格计数率，反映计数不为0的网格数量占总网格数量的比例。

对于采用通量信息生成全局权窗的MAGIC GVR方法，基本原理是权窗下限值th,i正比于通量，这种全局减方差方法已应用在了国际热核聚变实验反应堆(International Thermonuclear Experimental Reactor, ITER)及CFETR模型全局减方差计算中[2−4]。在低通量区域，毗连区域的通量衰减很快，当直接采用正比于通量的权窗具有很高的梯度变化，使中子过度分裂，而大大增加了一个中子历史的计算时间，所以当通量的相对误差Re超过限值时将关闭权窗或者统一地设定权窗值[2]。本文将相对误差限值设为0.5，当Re大于0.5时，权窗统一设置为一个固定值。

将MAGIC GVR方法应用于CFETR三维模型的中子通量分布计算中，进行多次迭代来得到精度较高的全空间中子通量分布。为缩短计算耗时，采用了并行计算的方法运行MCNP程序，计算中共用21个核（1个核充当master，20个核充当slaves）。CFETR三维模型的几何视图如图1所示。网格在、、三个方向上划分，、方向各有50个网格，方向有80个网格，网格尺寸为10 cm×24 cm×24 cm。用通量的网格计数结果生成全局权窗的公式，如式(1)所示，当Re大于等于0.5时，权窗统一设置为0.001。一共进行了4次迭代，计算用时、平均误差等结果见表1。可以看出，随着迭代次数的增加，一个历史平均用时增大，网格计数率Scoring增大，平均相对误差减小。网格计数率在第4次迭代时达到97.22%，接近100%，说明此次使用的全局权窗已经可以使中子输运到权窗网格覆盖的绝大部分空间。在并行效率方面，第4次迭代时并行效率达到最低，只有0.1189，通过定位随机数最多的历史编号，并用MCNP单独运行这一个历史，计算用时418.2 min，是第4次迭代总共用时1340.14 min的近1/3、一个历史平均用时0.0118 min的35000倍，可见“长历史”在此次迭代中是很明显的。

导致“长历史”问题的原因尚不明确，但有理由怀疑“长历史”是粒子过度分裂导致[8]。在具有真空细长管道的深穿透问题中，管道离粒子源最远一端的权窗比离源最近一端的权窗低很多，如果一个粒子恰好从近端进入真空管道，沿管道飞到远端，由于粒子在远端的权重远高于权窗，所以会过度分裂，产生大量粒子，使完成这个历史输运计算的时间增大。通过观察第四次迭代中“长历史”事件日志，发现粒子并未穿过包层之间的真空管道，而是在包层里正常输运。此外CFETR包层之间的真空管道两端的权窗差异不显著，约200倍的差距，如果粒子恰好穿过管道，也只会分裂为近200个粒子，不会使这个历史的用时比平均用时高35000倍。由此可以认为真空管道在这个问题中不是引起“长历史”的唯一原因。

表1 MAGIC GVR方法的计算结果

本文通过对“长历史”事件日志的观察提出了引起“长历史”的另一个原因。包层外存在上一次迭代的通量计数相对误差大于0.5的区域，生成权窗时直接将这些区域的权窗下限th设为0.001（见式(1)），而通过通量计数结果生成的th大部分在10−10以下，粒子在恰好th=0.001的区域多次在轮盘赌中存活，权重变得很大，然后输运到th<10−10的区域，通过多次分裂降低权重，引起“长历史”。记录第4次迭代“长历史”中粒子权重这一先升后降的事件日志如图2所示。第一个方框内的数列显示了粒子多次在轮盘赌中存活，从而权重从1.823×10−10上升到1.778×10−3。第二个方框内的数列显示了粒子输运到th<10−10的区域，多次分裂，权重由1.689×10−3下降到7.827×10−10。这个过程产生的过度分裂是显著的，足以引起“长历史”。当然粒子多次在轮盘赌中存活的概率也是极低的，需要计算大量的历史才能观察到这一现象。

为了克服“长历史”问题，对生成全局权窗的公式做一些改进，如式(2)所示。第一，对通量计数之比取0.9的指数，使屏蔽区权窗的梯度降低，减少粒子的分裂；第二，通量计数相对误差大于0.5的区域的权窗下限由0.001提升至10000，使粒子的权重由于轮盘赌上升至0.001后仍需继续轮盘赌，增大这些区域通过高权窗终结粒子的概率。

由这个新的全局权窗作为权窗输入文件，进行第4次迭代，设定CPU时间限CTME为3000 min，得到计算结果。将此结果与原来的第4次迭代的结果进行比较，如表2中原始权窗与新权窗两行数据所示。从表2中数据可以看出，采用新权窗后，相同CPU时间里计算的历史数量大大增加，FoMG变为原来的2.15倍，pFoMG变为原来的6.67倍，相对误差标准差也有所减小，网格计数率Scoring基本没变，计算效率有了明显的提升。

“长历史”方面，通过定位使用随机数最多的历史编号，并用MCNP单独运行这一个历史，计算用时38.6 min。虽然减小到原来418.2 min的约1/10，但仍然是一个历史平均用时9.2×10−4min的42000倍，整个计算用时422.49 min的1/10。并行计算效率提升至0.36，仍然不高。由这些判定，“长历史”问题仍然存在。可能粒子多次在轮盘赌中存活，权重上升到一定值（远小于10000，但远高于10−10）之后，通过输运离开th=10000的区域，进入到th<10−10的区域，通过权重先升后降产生过度分裂。总之，新权窗使计算效率有了明显提升，但未能有效克服“长历史”问题。“长历史”问题应从粒子在轮盘赌中存活的次数这个角度来解决。

图2 事件日志记录的粒子权重先升后降的现象

表2 三种情况GVR计算的参数比较

2 “长历史”问题的解决方案

就前面分析的“长历史”成因，提出一种由粒子在轮盘赌中存活次数来限制粒子分裂的解决“长历史”问题的方法。该方法的流程图如图3所示，其中：S表示源粒子的权重；C1是MCNP截断卡(Cutoffs Card)[1]的一个参数；MXSPLN和WSURVN是MCNP权窗参数卡(Weight Window Parameter Card, WWP)[1]的参数。MXSPLN规定粒子至少以1/MXSPLN存活的概率参与轮盘赌。为防止粒子连续在轮盘赌中存活从而权重变得很大，对粒子连续在轮盘赌中存活的次数WRN (Weight Roulette Number)进行统计，如果达到限定值，则禁止粒子参与分裂为MXSPLN粒子的权窗分裂。如果粒子以1/MXSPLN的概率在轮盘赌中存活，则WRN增加1。如果粒子输运到权窗为0的区域，参与权截断，以1/2的概率存活，则WRN增加0.5。采用这个方法，当粒子权重上升到一定程度后，就无法下降了，从而克服了粒子权重先升后降的问题，阻止粒子过度分裂。

图3 由粒子在轮盘赌中存活次数限制粒子分裂方法的流程图

3 结果与讨论

将上述功能加入MCNP程序中，用修改过的程序进行CFETR通量分布计算第4次迭代，全局权窗采用式(2)计算得到的新权窗。粒子在轮盘赌中存活的次数限在输入文件中输入，取1−10这10个限值进行测试，每次测试运行相同的粒子数NPS，观察“长历史”是否存在及计算效率的变化情况，并与不采用此功能得到的计算结果对比（即次数限为¥，计算结果的各项参数见表2中新的权窗）。采用次数限1−10计算结果的各项参数变化情况见图4。

从图4(a)可以看出，随着次数限的增大，CTM、时钟时间乘以处理器数量(Computer Time)逐渐增大。并行效率在次数限小于5时达到约0.8，明显好于不采用次数限的0.36，且在计算的过程没有明显的“卡死”现象，即没有出现“长历史”。在次数限为3的情况下，通过定位使用随机数最多的历史编号，并用MCNP单独运行这一个历史，计算用时2.03 min，减小到不采用次数限时38.6 min的约1/20，是一个历史平均用时6.6×10−4min的3000倍，整个计算用时135.45min的约3/200，已不再占有整个计算用时较大的比例，整个计算过程没有出现明显的“长历史”。当次数限达到5之后，并行效率明显下降。当次数限达到9时，并行效率下降到0.36，接近不采用次数限的并行效率。

从图4(b)可以看出，随着次数限的增大，通量网格计数的平均相对误差先减小后增大。在次数限较小时，粒子的权窗分裂受到很强的限制，可能出现粒子还未输运到屏蔽层外，在轮盘赌中存活的次数已达到限值，从而不能进行分裂为MXSPLN粒子的权窗分裂，这样使权重较高的粒子输运到低权窗区域，增大此区域计数的统计误差。当次数限增大时，粒子权窗分裂受到的限制减弱，故计数相对误差会减小。当次数限增大到6时，计数相对误差有小幅上升，可能是由于粒子在屏蔽区外权重先升后降的现象不能得到有效控制，粒子有一定概率由于在轮盘赌中存活，权重变得很大，再通过分裂产生大量粒子，使一定区域的通量计数偏离正常情况下的计数，增大了统计误差。总体来讲，粒子连续多次在轮盘赌中存活本身就是一个小概率事件，采用次数限对相对误差的影响总体较小，从图4(b)可以看出，通量计数的相对误差变化不大，也说明了这一点。

将原始权窗、新权窗、新权窗加设定次数限为3这三种情况下计算结果的参数总结于表2。从表2中的数据比较可以看出，采用式(2)得到的全局权窗后，Parallel Efficiency、FoMG、pFoMG变大，减小，Scoring基本没变，计算效率比之前有很大提升，在用新的权窗的基础上设定次数限后，Parallel Efficiency、FoMG、pFoMG变大，、Scoring基本没变，计算效率又得到了很大的提升。将这三种情况下变化明显的参数(Parallel Efficiency、FoMG、pFoMG)共同绘制在图5中。可以看出，采用新权窗和设定次数限后，这三个参数均提升明显，并行效率达到了约80%，克服了“长历史”问题。

图4 CFETR模型GVR计算各项参数随次数限变化

图5 三种情况GVR计算的典型参数比较

4 结语

针对采用GVR方法的中子输运计算中遇到的“长历史”问题，首先通过“长历史”的事件记录文件分析了其成因，即在计数相对误差大于0.5而被统一设置权窗下限的区域中，粒子连续在轮盘赌中存活，从而权重变得很大，再输运到相对误差小于0.5的低权窗区域，产生过度分裂。对权窗生成程序做了改进，利用新的权窗进行中子输运计算，发现品质因子(pFoMG)达到权窗改进前的5.6倍，但“长历史”仍然存在。接着，根据对“长历史”成因的分析，提出了由粒子在轮盘赌中存活的次数来限制粒子分裂的方法。在对权窗改进的基础上利用这个新方法，发现当次数限取1−4时并行效率达到了80%，“长历史”消失，品质因子(pFoMG)达到不采用次数限的约3倍，是不采用次数限且权窗改进前的约17倍。这说明“长历史”问题得到有效的解决，计算效率得到显著提高。由于粒子多次在轮盘赌中存活本身就是一个小概率事件，故这个方法对粒子正常输运影响较小，是否采用这个方法以及次数限的取值多少对减方差的效果影响较小，品质因子的增大主要是由于运行一定历史数(NPS)所需时间的减小和并行效率的增大导致的。

本文分析的“长历史”成因与模型的具体结构没有直接关系，不局限于文中提到的CFETR模型。通过采用本文提出的解决方法，可以有效解决“长历史”问题，在采用基于权窗的全局减方差计算中具有一定的普适性。

本文没有详细研究在计数相对误差大于0.5区域设置的权窗下限值的变化对计算效率和“长历史”问题的影响，也许存在巧妙的设置权窗的方法，也可以解决“长历史”问题。当然，除了本文分析的成因，可能还存在其他的原因导致“长历史”。对于“长历史”问题，今后需要更多的探索。

1 X-5 Monte Carlo Team. A general Monte Carlo N-particle transport code: user’s guide[R]. Version 5. Vol.II. USA: Los Alamos National Laboratory, 2003.

2 聂星辰, 李佳, 赵平辉, 等. 蒙特卡罗模拟CFETR中子输运计算中的全局减方差方法应用及对比[J]. 核技术, 2016, 39(3): 030501. DOI: 10.11889/j.0253-3219.2016. hjs.39.030501. NIE Xingchen, LI Jia, ZHAO Pinghui,. Application and comparison of global variance reduction methods employed in Monte Carlo neutron transport for CFETR[J]. Nuclear Techniques, 2016, 39(3): 030501. DOI: 10.11889/ j.0253-3219.2016.hjs.39.030501.

3 Van Wijk A J, van den Eynde G, Hoogenboom J E. An easy to implement global variance reduction procedure for MCNP[J]. Annals of Nuclear Energy, 2011, 38: 2496−2503. DOI: 10.1016/j.anucene.2011.07.037.

4 Davis A, Turner A. Comparison of global variance reduction techniques for Monte Carlo radiation transport simulations of ITER[J]. Fusion Engineering and Design, 2011, 86: 2698−2700. DOI: 0.1016/j.fusengdes.2011.01. 059.

5 Naish J, Fox F, Ghani Z,. Radiation mapping at JET and ITER using advanced computational acceleration techniques and tools[J]. Nuclear Technology, 2015, 192: 299−307. DOI: 10.13182/nt14-132.

6 聂星辰, 李佳, 赵平辉, 等. 深穿透屏蔽计算中MCNP 减方差技巧应用及比较[J]. 核电子学与探测技术, 2016, 36(7): 729−733, 741. DOI: 10.3969/j.issn.0258-0934. 2016.07.016. NIE Xingchen, LI Jia, ZHAO Pinghui,. Application and comparison of MCNP variance reduction techniques for deep penetration shielding calculation[J]. Nuclear Electronics & Detection Technology, 2016, 36(7): 729−733, 741. DOI: 10.3969/j.issn.0258-0934.2016.07. 016.

7 Turner A, Davis A. Improving computational efficiency of Monte-Carlo simulations with variance reduction[C]. International Conference on Mathematics and Computational Methods Applied to Nuclear Science & Engineering (M&C 2013), Sun Valley, Idaho, USA, 2013.

8 Wasastjerna F. Using MCNP for fusion neutronics[D]. Finland: VTT Technical Research Centre of Finland, 2008.

A solution to “long history” in the method of GVR based on weight window

ZHANG Hengquan1LI Jia1ZHU Qingjun2WU Yuxiao3NIE Xingchen1

1(School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230027, China) 2(Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, China) 3(Guangxi Fangchenggang Nuclear Power Co., Ltd., Fangchenggang 538000, China)

The“long history” induced by the method of global variance reduction (GVR) based on weight window would reduce the efficiency of the Monte Carlo particle transport calculation with global characteristic.To enhance the efficiency and minimize the waste of computing resource, this article proposes a solution to “long history” based on the study of its causes.According to the event log of “long history” in Chinese fusion engineering testing reactor (CFETR) neutron transport calculation, the cause of the “long history” is brought up, which is in relation to particle’s successively surviving in the Russian roulette games. A solution to the “long history” by limiting particle’s splitting according to the number of surviving games (weight roulette number, WRN), is came up with, which is tested in CFETR model.The result shows that when the limited number is between 1 and 4, the parallel efficiency reaches 80%, with no “long history” present and global figure of merit (pFoMG) being around 3 times of that in calculation without the limited number.The solution proposed prevented the “long history” and thus greatly enhanced the efficiency of GVR.

GVR, Global weight window, Long history, Russian roulette, Limited number

ZHANG Hengquan, male, born in 1995, undergraduate, major in nuclear engineering and technology

LI Jia, E-mail: lijia@ustc.edu.cn

2017-05-05,

2017-06-30

TL62

10.11889/j.0253-3219.2017.hjs.40.120602

张恒权，男，1995年出生，目前就读于中国科学技术大学核科学技术学院攻读学士学位，核工程与核技术专业

李佳，E-mail: lijia@ustc.edu.cn

2017-05-05，

2017-06-30

Supported by the Fundamental Research Funds for the Central Universities (No.WK2140000011)