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

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

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

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

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

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

區(qū)域分解方法是指將研究對象從幾何上劃分為若干(子)區(qū)域，對不同的(子)區(qū)域進(jìn)行計算，并通過建立(子)區(qū)域間的某種耦合關(guān)系使得對整個研究對象的求解結(jié)果正確?！胺侄沃笔菂^(qū)域分解的基本思想。它可將大模型問題化為小模型問題，是解決計算條件不足問題的一種有效方案。

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

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

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

2.1 實現(xiàn)方法

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

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

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

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

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

2.2 測試

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

1) 立方體模型

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

2) 球模型

表2列出球模型的計算結(jié)果，其中編號代表不同的劃分方式：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) 組件模型

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

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

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

表1 立方體模型的計算結(jié)果

表2 球模型的計算結(jié)果

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

表3 組件模型的計算結(jié)果

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

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

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

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

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

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

要實現(xiàn)結(jié)果的可重復(fù)，必須保證所有粒子的運動歷史一致，進(jìn)而給每一個源中子(每一代)分配的初始隨機(jī)數(shù)種子必須一致。在源中子順序一致的情況下，通過跳躍法預(yù)先產(chǎn)生分段的隨機(jī)數(shù)種子，逐一分配給每一個源中子可確保粒子的運動歷史不變，這與常見的粒子并行算法結(jié)果可重復(fù)性使用的方法相同。然而，在區(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#。當(dāng)區(qū)域被分為左、右兩個子區(qū)域進(jìn)行計算時，A、B中子同時開始運動(分別在兩個處理器上)，新一代源中子的產(chǎn)生順序為1#、3#、2#，亦即2#與3#中子的存放順序發(fā)生了變化。

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

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

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

4 小結(jié)

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

參考文獻(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.