中子輸運蒙特卡羅模擬的區(qū)域分解方法研究

梁金剛，王 侃，蔡 云,2，孫嘉龍

(1.清華大學 工程物理系，北京 100084；2.中國核動力研究設計院，四川 成都 610041)

近年來，在反應堆物理計算研究領域，求解中子輸運問題的蒙特卡羅(簡稱蒙卡)模擬方法越來越受到重視，正在成為研究的熱點和前沿。蒙卡方法作為一種基于概率論的隨機模擬方法，能有效借助計算機的速度與精度優(yōu)勢，在模擬粒子輸運方面表現(xiàn)出忠實物理過程程度高(如連續(xù)點截面)、受幾何限制小、對問題維度不敏感等諸多優(yōu)點，恰好符合反應堆中子輸運計算的幾何與能譜復雜、多維度、多變量等特點。因而，相較于中子輸運方程計算的確定論方法，蒙卡方法被認為是未來反應堆工程計算方法的一種重要選擇[1-2]。

然而，將蒙卡方法應用于反應堆全堆芯的大規(guī)模計算分析還有一些問題需進一步解決，其中，計算時間長和內(nèi)存占用大是兩個主要問題，著名的“Kord Smith挑戰(zhàn)”[3]及Bill Martin所做的修正[1]說明了這一點。針對計算時間長采取并行計算是有效的解決方法。由于蒙卡方法是對大量粒子進行逐一的跟蹤模擬，按照粒子數(shù)的并行容易實現(xiàn)，已有研究表明，粒子并行方法具有良好的并行效率和可擴展性[4]，可有效解決蒙卡方法計算時間長的限制。而對于內(nèi)存占用大的問題，目前還沒有較成熟的解決方法，區(qū)域分解方法被認為是一種可行的思路。

文獻[5]最早提出將區(qū)域分解應用于隱蒙卡的并行化，并研究了區(qū)域分解的實現(xiàn)方法。但區(qū)域分解在粒子輸運模擬中的應用并未得到推廣，在眾多著名的蒙卡程序中，鮮有程序具備區(qū)域分解功能。近年來，隨著蒙卡粒子輸運模擬的興起和反應堆大規(guī)模計算中內(nèi)存問題的凸顯，對區(qū)域分解的研究逐步增多。其中，美國LLNL對區(qū)域分解方法進行了一些研究，其開發(fā)的蒙卡程序Mercury[6]具備一定的區(qū)域分解功能。在國內(nèi)，對區(qū)域分解的研究屬于起步階段，北京應用物理與計算數(shù)學研究所基于JMCT程序進行了區(qū)域分解的研究[7]。

1 區(qū)域分解方法的基本原理

區(qū)域分解方法是指將研究對象從幾何上劃分為若干(子)區(qū)域，對不同的(子)區(qū)域進行計算，并通過建立(子)區(qū)域間的某種耦合關系使得對整個研究對象的求解結果正確。“分而治之”是區(qū)域分解的基本思想。它可將大模型問題化為小模型問題，是解決計算條件不足問題的一種有效方案。

區(qū)域分解作為數(shù)學方法，最早應用于偏微分方程的求解，其提出目的就是為了應對當時計算機計算條件的限制。求解偏微分方程的區(qū)域分解方法如圖1所示，一般按照區(qū)域劃分方式分為不重疊型和重疊型兩種。圖1中，Ω為計算區(qū)域，Γ為區(qū)域交界，下標1、2表示區(qū)域劃分為兩個子區(qū)域。在計算流程上，以各區(qū)域交替計算(如經(jīng)典的Schwarz交替法)為主要特點，即先計算子區(qū)域1的值(假設交界邊界面初值)，然后將新的邊界值傳遞至區(qū)域2，計算區(qū)域2的值，再次更新邊界值，重復上述步驟，直到兩個區(qū)域的解收斂到真實解。

區(qū)域分解方法的思想同樣可用于粒子輸運蒙卡模擬的過程。粒子在不同的子區(qū)域運動，當從一個子區(qū)域進入到另一子區(qū)域時，通過對粒子運動狀態(tài)進行存儲及區(qū)域間的傳遞，實現(xiàn)區(qū)域間的耦合計算。一般情況下，不同的子區(qū)域由不同的處理器計算，利用蒙卡模擬中粒子的相互獨立性，可實現(xiàn)不同子區(qū)域的并行計算，因而區(qū)域分解方法又稱為空間并行方法。

2 基于RMC的區(qū)域分解實現(xiàn)及測試

2.1 實現(xiàn)方法

RMC[8]是由清華大學工程物理系核能科學與工程管理研究所反應堆工程計算分析實驗室(REAL團隊)自主開發(fā)的用于反應堆物理分析的三維粒子輸運蒙卡程序。本文基于RMC平臺進行區(qū)域分解功能的研究和開發(fā)。

區(qū)域分解的實現(xiàn)主要包括兩部分內(nèi)容：區(qū)域劃分和區(qū)域間通信。前者指對模擬區(qū)域進行剖分，并將每個區(qū)域分配至不同的處理器，在區(qū)域劃分方法和區(qū)域與處理器對應關系上進行探究以實現(xiàn)優(yōu)化；后者指粒子穿越邊界時，區(qū)域之間對粒子信息進行存儲、傳遞的通信算法，蒙卡模擬中粒子穿出區(qū)域多、通信頻繁、信息量大，不同的通信算法對區(qū)域分解的性能會有影響。

a——不重疊型；b——重疊型

粒子輸運蒙卡模擬的區(qū)域分解方法如圖2所示，每個區(qū)域由一個處理器模擬，當粒子穿出當前所在區(qū)域時，處理器立即保存粒子信息(包括位置、方向、速度及隨機數(shù)種子等)至緩存區(qū)，繼續(xù)模擬下一個粒子，直到當前區(qū)域的所有粒子模擬完畢，不同區(qū)域開始交換信息，將穿出區(qū)域的粒子發(fā)送到對應的區(qū)域(處理器)，然后進行下一輪模擬，直到所有粒子至歷史結束。

本文按照上述方法在RMC中添加區(qū)域分解模塊，在三維幾何下采用曲面組合邏輯的用戶輸入方式進行不重疊型的區(qū)域劃分，即用戶輸入?yún)^(qū)域邊界，根據(jù)布爾運算逐個指定子區(qū)域范圍。這種方式不影響計算模型原始輸入，同時具備幾何區(qū)域劃分的靈活性。在區(qū)域間通信方面，基于MPI并行庫進行粒子信息傳遞。為避免通信發(fā)生死鎖，在通信中根據(jù)區(qū)域編號規(guī)定信息發(fā)送、接收的順序。

2.2 測試

為了對實現(xiàn)的區(qū)域分解基本功能進行測試，設計了3個算例模型(立方體模型、球模型、組件模型)，如圖3所示，按照不同的區(qū)域劃分方式進行計算，將結果與未進行區(qū)域分解的RMC和MCNP的計算結果進行對比，并給出區(qū)域分解的并行加速效率。

1) 立方體模型

立方體模型采用均勻材料，對模型進行對稱的2區(qū)域(分別沿x、y、z3個方向)、4區(qū)域和8區(qū)域劃分，臨界計算結果列于表1。計算條件為模擬300代，前100代為非活躍代，每代20 000個中子。測試平臺為Windows系統(tǒng)，12核Intel Xeon CPU(型號X5670，2.93 GHz，下同)。

2) 球模型

表2列出球模型的計算結果，其中編號代表不同的劃分方式：so0為無區(qū)域分解；so1為2區(qū)域均分(沿x方向)；so2為4區(qū)域均分(x、y4個象限均分)；so3為2區(qū)域劃分(小球內(nèi)、外)；so4為2區(qū)域劃分(大球內(nèi)、外)；so5為6區(qū)域劃分(4個小球內(nèi)、大球內(nèi)小球外和大球外)。計算條件為模擬300代，前100代為非活躍代，每代50 000個中子。

圖2 粒子輸運蒙卡模擬的區(qū)域分解方法

a——立方體模型；b——球模型；c——組件模型

3) 組件模型

選取重水堆組件進行計算，按照圖4所示方式進行2～24區(qū)域劃分計算，其中2～12區(qū)域劃分為沿圓柱橫切面的角平分面進行劃分，24區(qū)域劃分是在12區(qū)域劃分的基礎上在圓柱軸向中間橫切。計算結果列于表3。計算條件為模擬300代，前100代為非活躍代，每代50 000個中子。

2.3 區(qū)域分解的結果及并行性能分析

上述3個測試模型分別采用不同的區(qū)域劃分方式，包括不同的子區(qū)域數(shù)目、子區(qū)域幾何類型、子區(qū)域劃分方向等多種情形。計算結果表明，臨界計算有效增殖因數(shù)keff與無區(qū)域分解的結果均在統(tǒng)計漲落誤差內(nèi)(相差均在2倍標準偏差以內(nèi))。

表1 立方體模型的計算結果

表2 球模型的計算結果

圖4 組件模型區(qū)域劃分方式

表3 組件模型的計算結果

在區(qū)域分解情況下，粒子隨機數(shù)種子的分配相對于無區(qū)域分解時發(fā)生了變化，導致最終結果并非完全一致，但這種分配的變化是隨機的，不會影響結果的收斂與準確性(蒙卡方法本身具有統(tǒng)計性)。相對誤差總是在2倍標準偏差以內(nèi)，說明誤差是因隨機數(shù)不同產(chǎn)生的統(tǒng)計性漲落導致的，因此這種結果是可靠的。通過進一步的設計可實現(xiàn)兩者結果完全一致，即區(qū)域分解結果的可重復性。

3個測試均采取處理器與區(qū)域一一對應的方式，區(qū)域分解后，處理器并行計算，因而整體縮短了計算時間。對比不同劃分方式的加速效果，可以看出，區(qū)域劃分均衡情況下，具有較好的并行效率，如立方體模型及組件模型(圖5，加速比為串行計算與并行計算之比)；反之，區(qū)域劃分不均衡時，并行效率較差，如球模型。另一方面，隨著區(qū)域數(shù)目增多，區(qū)域間通信增多，并行效率將降低。負載均衡與通信性能是影響區(qū)域分解算法性能的關鍵因素。

圖5 組件模型對稱區(qū)域劃分下的并行加速比

3 區(qū)域分解結果可重復性研究

由于隨機數(shù)對源中子的分配順序發(fā)生變化，使得粒子運動徑跡改變，以致計算結果相對于不分解計算出現(xiàn)偏差。為驗證程序的正確性，設計出在輸入一致的條件下區(qū)域分解的計算結果與無分解結果一致的算法是有必要的。

圖6 區(qū)域分解下新一代源中子產(chǎn)生順序的變化

要實現(xiàn)結果的可重復，必須保證所有粒子的運動歷史一致，進而給每一個源中子(每一代)分配的初始隨機數(shù)種子必須一致。在源中子順序一致的情況下，通過跳躍法預先產(chǎn)生分段的隨機數(shù)種子，逐一分配給每一個源中子可確保粒子的運動歷史不變，這與常見的粒子并行算法結果可重復性使用的方法相同。然而，在區(qū)域分解算法中，由于不同區(qū)域的粒子可同時運動，穿出區(qū)域、尚未死亡的粒子需要等待，被傳遞后才能繼續(xù)運動，以至新產(chǎn)生的下一代源中子的順序被打亂。圖6示出區(qū)域分解下新一代源中子產(chǎn)生順序的變化。由圖6可見，A中子在運動中，先后在左、右區(qū)域各產(chǎn)生1個新中子(分別記為1#、2#)，B中子在右區(qū)域產(chǎn)生1個新中子(記為3#)。在無區(qū)域分解情況下，新中子的產(chǎn)生順序為1#、2#、3#。當區(qū)域被分為左、右兩個子區(qū)域進行計算時，A、B中子同時開始運動(分別在兩個處理器上)，新一代源中子的產(chǎn)生順序為1#、3#、2#，亦即2#與3#中子的存放順序發(fā)生了變化。

為保證下一代計算結果不變，需使源中子的隨機數(shù)分配不變，因而需將順序變化的源中子重新排序。為此，對于每個新中子，可記錄致其產(chǎn)生的舊中子序號和已由此舊中子碰撞產(chǎn)生的新中子數(shù)目，據(jù)此便可對所有源中子進行排序，使得源中子順序與無區(qū)域分解情況一致。實際計算中，也可根據(jù)新源中子的順序對隨機數(shù)種子序列進行排序使得兩者符合一致。

在實現(xiàn)結果可重復性后，區(qū)域分解與不區(qū)域分解的程序，兩者模擬的粒子歷史完全一致，因而程序的統(tǒng)計結果是完全一致的，包括keff、通量等。表4列出按照上述原理對程序修改后，對立方體模型2區(qū)域分解下的臨界計算得到的keff。計算條件為模擬300代，前100代為非活躍代，每代5 000個中子。由表4可見，2區(qū)域與單區(qū)域的keff計算結果完全一致。需注意的是，在確保結果一致的計算中增加了排序過程，因而計算時間有所增加。

表4 區(qū)域分解實現(xiàn)結果的可重復性

4 小結

對于中子輸運的蒙卡模擬方法，區(qū)域分解是解除內(nèi)存限制的值得探索的方向。本文介紹了區(qū)域分解的基本思想，分析其實現(xiàn)方法，并基于RMC開發(fā)區(qū)域分解功能，進行了算例測試和性能分析，總結出影響區(qū)域分解并行性能的關鍵因素，即負載均衡和通信性能。最后研究了區(qū)域分解結果可重復性的實現(xiàn)方法，提出通過對源中子排序以保證粒子運動歷史的一致性。

參考文獻：

[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.