李 靜 宋 婧 龍鵬程 劉鴻飛 江 平

基于雙向遍歷的蒙特卡羅臨界計算并行方法

李 靜1宋 婧2龍鵬程2劉鴻飛2江 平1

1（合肥工業(yè)大學 數(shù)學學院 合肥 230009）
2（中國科學院核能安全技術研究所 中國科學院中子輸運理論與輻射安全重點實驗室 合肥 230031）

在基于蒙特卡羅粒子輸運方法的反應堆模擬中，如裂變堆、聚變裂變混合堆等，達到可接受的統(tǒng)計誤差需要大量的計算時間，這已成為蒙特卡羅方法的挑戰(zhàn)問題之一，需通過并行計算技術解決。為解決現(xiàn)有方法中通信死鎖的問題并保證負載均衡性，設計了基于雙向遍歷的臨界計算并行算法。該方法基于超級蒙特卡羅核計算仿真軟件系統(tǒng)SuperMC進行實現(xiàn)，以池式鈉冷快堆BN600基準模型進行驗證，并與MCNP進行對比。測試結果表明，串行和并行計算結果一致，且SuperMC并行效率高于MCNP。

臨界，并行計算，蒙特卡羅，SuperMC

近年來層次化并行、跨節(jié)點消息通信的多核處理器上的組合線程計算等計算機技術極大地促進了蒙特卡羅方法的發(fā)展與在反應堆中的實際應用[1]。然而，在反應堆尤其是復雜先進反應堆[2–6]應用中，蒙特卡羅方法存在巨大的挑戰(zhàn)。2012年，Martin[7]對全堆芯蒙特卡羅有效增殖因子(keff)模擬等面臨的挑戰(zhàn)進行了全面的分析，包括：(1) 達到可接受的統(tǒng)計誤差需要大量的計算時間；(2) 海量計算內存；(3) 慢的裂變源收斂；(4) 顯示方差不能反映真實方差；(5) 適應未來計算機架構，其中(1)和(2)需要依賴于并行計算技術解決。蒙特卡羅方法天然具有較好的并行粒度，基于蒙特卡羅的反應堆粒子輸運計算可分為固定源計算和臨界計算。固定源計算按源粒子依次模擬，相互獨立，一般能獲得較好的并行效率；基于指數(shù)迭代法的臨界計算需要逐代模擬，每代之間具有相關性，需進行節(jié)點及處理器之間的負載均衡才能得到較好的并行效率。

keff是臨界計算中最本質的物理量，目前求解時是將基本的中子輸運方程轉換成指數(shù)迭代形式的k特征值方程：

式中，L是泄露運算符；S是內散射運算符；T是碰撞運算符；M是裂變增殖運算符。采用蒙特卡羅方法進行模擬時，從初始假定的值keff(0)和源分布φ(0)開始，采用碰撞、吸收、徑跡長度估計及其組合估計等方法來估計某一代的的值，沿粒子徑跡抽樣下一代裂變源分布并作為下一代初始源迭代直至收斂。

臨界計算基于粒子并行的基本思路是：在某代粒子輸運過程中，各個計算單元抽樣并存儲裂變源，該代粒子模擬完后各計算單元將裂變源歸并到主計算單元再分配到各計算單元進行下一代計算。該傳統(tǒng)方法中，由于部分裂變源粒子是存儲在當前計算單元的，因此存在大量不必要的主計算單元與從計算單元之間的粒子通信，降低了并行效率。Romano等[8]進行了方法改進，保留本地產生的裂變源粒子，先計算并存儲平均后粒子的序號，據(jù)此進行計算單元之間的裂變源粒子傳遞。

本文在現(xiàn)有方法基礎上開展研究，并對其中的通信死鎖問題進行改進，介紹了基于雙向遍歷的臨界并行計算算法設計。為保障及檢驗并行的正確性，對串行隨機數(shù)按臨界特征值計算方法調整隨機數(shù)序列，介紹相應的方法，并采用基于池式鈉冷快堆BN600進行方法的驗證對比。

1 臨界并行計算算法設計

在計算單元之間進行裂變源粒子傳輸時，需要通過MPI (Multi Point Interface)進行通信，通信接口采用阻塞模式，相鄰計算單元間在通信時進程必須掛起，等待數(shù)據(jù)傳輸完成后，兩個進程才能進行其它工作。按Romano[8]的算法對每個計算單元進行粒子傳輸時，當某個計算單元對相鄰節(jié)點提出的粒子收發(fā)需求不能滿足時，該計算單元以及和它通信的計算單元的進程將同時掛起，發(fā)生死鎖。這里以案例說明該情況：假設在使用10個核完成某代keff計算后，產生了總共100個下代源粒子，每個進程源粒子分布依次為：21、5、0、0、10、28、10、3、2、21。1號進程向2號進程傳遞11個源粒子，2號進程向3號進程傳遞6個源粒子，此時3號進程擁有6個源粒子，不滿足10個粒子的需求，當其向4號進程索取4個源粒子時，3號進程和4號進程將進入死鎖狀態(tài)，程序將停留在該位置。

為避免該問題，本文改為雙向遍歷的方式，設計的臨界計算并行的主要算法仍以上述案例為例，該流程的示意圖如圖1所示。

(1) 將存儲的裂變源粒子廣播到所有的計算單元。如果把所有的裂變源粒子存儲看作是一個分布在多個計算單元之間的大數(shù)組，那么這一步實際上是確定某個處理單元(ith)在該數(shù)組中裂變源粒子的開始序號(ai)和結束序號(bi)。

(2) 每個處理單元對存儲的裂變源粒子進行模擬，在每個處理單元創(chuàng)建一個單獨的數(shù)組，用于存儲該處理單元模擬粒子過程中抽樣產生的下一代裂變源粒子，并更新該處理器中裂變源粒子在整個數(shù)組中的編號ai與bi。

(3) 從左向右遍歷所有計算單元，并使計算單元僅向右邊計算單元發(fā)送粒子，僅從左邊計算單元接收粒子。即判斷：如果bi大于(i+1)N/p，則向右邊鄰近的處理單元發(fā)送bi?(i+1)N/p個裂變源粒子，其中N是該代總的裂變源粒子數(shù)；p是處理單元的個數(shù)；i是處理單元的序號；如果ai大于iN/p，則從左邊鄰近的處理單元接收ai?(i+1)N/p個裂變源粒子。

(4) 從右向左遍歷所有計算單元，并使計算單元僅向左邊計算單元發(fā)送粒子，僅從右邊計算單元接收粒子。即判斷：如果ai小于iN/p，則需要發(fā)送iN/p?ai個裂變源粒子到左邊鄰近的處理器單元；如果bi小于(i+1)N/p，則需要從右邊鄰近的處理單元接收(i+1)N/p?bi個裂變源粒子。

圖1 基于雙向遍歷的臨界并行計算流程示意圖Fig.1 Flow chart of criticality parallel calculation based on bi-directional traversal.

2 串并行一致的臨界計算中隨機數(shù)調整方法

為保證臨界并行計算結果串行和并行一致，本文采用乘同余法和斐波納契法組合方法產生隨機數(shù)，并用跳躍法確定每段隨機數(shù)初始種子。乘同余法的遞推公式如下：

式中，t0為初始種子；M1為模，通常設定為2的整數(shù)次冪；Tn為0?1之間均勻分布的隨機數(shù)。

斐波納契方法產生隨機數(shù)遞推公式如下：

式中，f0和f1為隨機數(shù)初始種子；M2為模；Fn為0?1區(qū)間均勻分布的隨機數(shù)。

產生的組合隨機數(shù)序列可表示如下，其隨機數(shù)周期為258：

某代模擬時，每個計算單元上裂變源粒子序號相連，對于每個源粒子及其次級粒子使用的隨機數(shù)最大序列長度固定。開始計算之前需要初始化每個源粒子的隨機數(shù)種子，即每個計算單元使用的第一個隨機數(shù)。本文使用跳躍法產生隨機數(shù)種子：乘同余法隨機數(shù)利用跳躍法產生隨機數(shù)種子時，使用式(2)、(3)產生，此時可用式(7)代替進行求模的位運算。斐波納契方法隨機數(shù)利用跳躍法產生隨機數(shù)種子時，先使用式(8)斐波納契通項公式算出第n項和n?1項隨機數(shù)值，再代入式(9)得到。最后使用式(6)得到組合隨機數(shù)的種子。

3 池式鈉冷快堆BN600基準模型驗證

上述方法在超級蒙特卡羅核模擬軟件SuperMC[9–13]中進行實現(xiàn)，并采用BN600基準模型進行驗證。BN600[14]是俄羅斯的池式鈉冷商用快堆，其熱功率為1 470MW，電功率為600MW。本文采用國際原子能機構(International Atomic Energy Agency, IAEA)發(fā)布的BN600基準模型進行程序驗證，圖2為BN600在SuperMC中的軸向視圖和徑向視圖。BN600的306個氧化鈾燃料組件依富集度布置在4個區(qū)：低富集度的內外兩個區(qū)(Low Emission Zone, LEZ)、中間的中等富集度區(qū)(Middle Emission Zone, MEZ)和外面的高富集度區(qū)(High Emission Zone, HEZ)。在MEZ和HEZ之間還有一個MOX燃料區(qū)（90個燃料組件）。有19根補償棒(Shim rod, SHR)插入到堆芯的中心平面，6根安全棒(Safety control rod, SCR)提出堆芯平面以上5.5cm。堆芯徑向外側有300個不銹鋼反射層組件，最外側是102個碳化硼屏蔽組件。其詳細的幾何布置、各區(qū)材料、幾何尺寸描述見文獻[15]。

圖2 BN600模型在SuperMC中的徑向視圖與軸向視圖 (a) 徑向視圖，(b) 軸向視圖Fig.2 Radial and axial view of BN600 model in SuperMC. (a) Radial view, (b) Axial view

使用2.2GHz 4核處理器，2×1 000M以太網(wǎng)，單節(jié)點32核的高性能計算集群進行測試，每代模擬十萬個粒子，共計算500代，選擇相同的數(shù)據(jù)庫進行對比測試。SuperMC串并行計算結果一致，keff均為1.01708(0.00005)，MCNP串并行結果一致，keff均為1.01709(0.00004)。其中，MCNP并行設計模式是主從設計模式(Master/Slave)：Master節(jié)點只負責任務平均分發(fā)（幾何信息、計數(shù)信息、材料信息等）和結果匯總；Slave節(jié)點負責計算各自的任務。Master-Slave Algorithm基本滿足臨界并行計算的需要，但存在兩個問題：一是無效的通信耗時太多；二是需要耗費很大的緩沖空間，這些問題會隨著計算任務和計算核數(shù)的增加更加明顯[15]。相比Master/Slave算法，雙向遍歷算法減少了大量的無效通信傳輸，減少了每代計算完成后的平均下代粒子的耗時。并且隨著計算負載和計算核的增加，該算法仍能保證很好的性能。表1為SuperMC和MCNP在不同核數(shù)下的并行效率對比，得出基于本文方法的SuperMC并行效率高于MCNP。

表1 BN600例題SuperMC與MCNP并行效率對比Table 1 Comparison of parallel efficiency of SuperMC and MCNP with BN600.

4 結語

本文對臨界計算并行方法開展研究，設計了基于雙向遍歷的并行方法，并采用俄羅斯池式鈉冷商用快堆BN600模型進行校驗。keff計算結果串并行一致，說明了本方法與實現(xiàn)的正確性；SuperMC計算的并行效率高于MCNP，說明本方法的有效性及更好的可擴展性。

致謝 感謝中國科學院核能安全技術研究所FDS團隊其他成員對本文工作的支持幫助。

1 Brown F B. Recent advances and future prospects for Monte Carlo[J]. Progress in Nuclear Science and Technology, 2011, 2: 1–4

2 Wu Y, FDS Team. Conceptual design activities of FDS series fusion power plants in China[J]. Fusion Engineering and Design, 2006, 81(23–24): 2713–2718

3 Wu Y, Qian J, Yu J, et al. The fusion-driven hybrid system and its material selection[J]. Journal of Nuclear Materials, 2002, 307–311: 1629–1636

4 Wu Y, FDS Team. Fusion-based hydrogen production reactor and its material selection[J]. Journal of Nuclear Materials, 2009, 386–388: 122–126

5 Wu Y, Jiang J, Wang M, et al. A fusion-driven subcritical system concept based on viable technologies[J]. Nuclear Fusion, 2011, 51(10): 103036

6 Wu Y, FDS Team. Conceptual design of the China fusion power plant FDS-II[J]. Fusion Engineering and Design, 2008, 83(10–12): 1683–1689

7 Martin W R. Challenges and prospects for whole-core Monte Carlo analysis[J]. Nuclear Engineering and Technology, 2012, 44(2): 151–160

8 Romano P K, Forget B. Parallel fission bank algorithms in Monte Carlo criticality calculations[J]. Nuclear Science and Engineering, 2012, 170(2): 125–135

9 Wu Y C, Song J, Zheng H Q, et al. CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC[J]. Annals of Nuclear Energy, 2014, 8: 58–65

10 Wu Y, FDS Team. CAD-based interface programs for fusion neutron transport simulation[J]. Fusion Engineering and Design, 2009, 84(7–11): 1987–1992

11 Li Y, Lu L, Ding A, et al. Benchmarking of MCAM4.0 with the ITER 3D model[J]. Fusion Engineering and Design, 2007, 82: 2861–2866

12 Hu H, Wu Y, Chen M, et al. Benchmarking of SNAM with the ITER 3D model[J]. Fusion Engineering and Design, 2007, 82: 2867–2871

13 Wu Y, Xie Z, Fischer U. A discrete ordinates nodal method for one-dimensional neutron transport calculation in curvilinear geometries[J]. Nuclear Science and Engineering, 1999, 133(3): 350–357

14 IAEA. BN-600 hybrid core benchmark analyses[R]. Austria: IAEA-TECDOC-1623, 2010

15 Brown F B, Goorley J T, Sweezy J E. MCNP5 parallel processing workshop[R]. ANS Mathematics & Computation Topical Meeting, 2003

CLC TL329.2

Parallel computing method of Monte Carlo criticality calculation based on bi-directional traversal

LI Jing1SONG Jing2LONG Pengcheng2LIU Hongfei2JIANG Ping1
1(School of Mathematics, Hefei University of Technology, Hefei 230009, China)
2(Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei 230031, China)

Background: It requires much computational time with acceptable statistics errors in reactor simulations including fission reactors and fusion-fission hybrid reactors, which has become one challenge of the Monte Carlo method. Purpose: In this paper, an efficient parallel computing method was presented for resolving the communication deadlock and load balancing problem of current methods. Methods: The parallel computing method based on bi-directional traversal of criticality calculation was implemented in super Monte Carlo simulation program (SuperMC) for nuclear and radiation process. The pool-type sodium cooled fast reactor BN600 was proposed for benchmarking and was compared with MCNP. Results: Results showed that the parallel method and un-parallel methods were in agreement with each other. Conclusion: The parallel efficiency of SuperMC is higher than that of MCNP, which demonstrates the accuracy and efficiency of the parallel computing method.

Criticality, Parallel computing, Monte Carlo, SuperMC

TL329.2

10.11889/j.0253-3219.2015.hjs.38.060602

中國科學院戰(zhàn)略性先導科技專項(No.XDA03040000)、國家ITER 973計劃(No.2014GB112001)、國家自然科學基金(No.91026004)資助

李靜，女，1989年出生，2015年于合肥工業(yè)大學獲碩士學位，計算數(shù)學專業(yè)

宋婧，E-mail: jing.song@fds.org.cn

2015-01-29，

2015-03-20