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

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

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

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

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

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

臨界，并行計(jì)算，蒙特卡羅，SuperMC

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

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

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

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

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

1 臨界并行計(jì)算算法設(shè)計(jì)

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

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

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

(2) 每個(gè)處理單元對存儲(chǔ)的裂變源粒子進(jìn)行模擬，在每個(gè)處理單元?jiǎng)?chuàng)建一個(gè)單獨(dú)的數(shù)組，用于存儲(chǔ)該處理單元模擬粒子過程中抽樣產(chǎn)生的下一代裂變源粒子，并更新該處理器中裂變源粒子在整個(gè)數(shù)組中的編號(hào)ai與bi。

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

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

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

2 串并行一致的臨界計(jì)算中隨機(jī)數(shù)調(diào)整方法

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

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

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

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

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

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

3 池式鈉冷快堆BN600基準(zhǔn)模型驗(yàn)證

上述方法在超級(jí)蒙特卡羅核模擬軟件SuperMC[9–13]中進(jìn)行實(shí)現(xiàn)，并采用BN600基準(zhǔn)模型進(jìn)行驗(yàn)證。BN600[14]是俄羅斯的池式鈉冷商用快堆，其熱功率為1 470MW，電功率為600MW。本文采用國際原子能機(jī)構(gòu)(International Atomic Energy Agency, IAEA)發(fā)布的BN600基準(zhǔn)模型進(jìn)行程序驗(yàn)證，圖2為BN600在SuperMC中的軸向視圖和徑向視圖。BN600的306個(gè)氧化鈾燃料組件依富集度布置在4個(gè)區(qū)：低富集度的內(nèi)外兩個(gè)區(qū)(Low Emission Zone, LEZ)、中間的中等富集度區(qū)(Middle Emission Zone, MEZ)和外面的高富集度區(qū)(High Emission Zone, HEZ)。在MEZ和HEZ之間還有一個(gè)MOX燃料區(qū)（90個(gè)燃料組件）。有19根補(bǔ)償棒(Shim rod, SHR)插入到堆芯的中心平面，6根安全棒(Safety control rod, SCR)提出堆芯平面以上5.5cm。堆芯徑向外側(cè)有300個(gè)不銹鋼反射層組件，最外側(cè)是102個(gè)碳化硼屏蔽組件。其詳細(xì)的幾何布置、各區(qū)材料、幾何尺寸描述見文獻(xiàn)[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é)點(diǎn)32核的高性能計(jì)算集群進(jìn)行測試，每代模擬十萬個(gè)粒子，共計(jì)算500代，選擇相同的數(shù)據(jù)庫進(jìn)行對比測試。SuperMC串并行計(jì)算結(jié)果一致，keff均為1.01708(0.00005)，MCNP串并行結(jié)果一致，keff均為1.01709(0.00004)。其中，MCNP并行設(shè)計(jì)模式是主從設(shè)計(jì)模式(Master/Slave)：Master節(jié)點(diǎn)只負(fù)責(zé)任務(wù)平均分發(fā)（幾何信息、計(jì)數(shù)信息、材料信息等）和結(jié)果匯總；Slave節(jié)點(diǎn)負(fù)責(zé)計(jì)算各自的任務(wù)。Master-Slave Algorithm基本滿足臨界并行計(jì)算的需要，但存在兩個(gè)問題：一是無效的通信耗時(shí)太多；二是需要耗費(fèi)很大的緩沖空間，這些問題會(huì)隨著計(jì)算任務(wù)和計(jì)算核數(shù)的增加更加明顯[15]。相比Master/Slave算法，雙向遍歷算法減少了大量的無效通信傳輸，減少了每代計(jì)算完成后的平均下代粒子的耗時(shí)。并且隨著計(jì)算負(fù)載和計(jì)算核的增加，該算法仍能保證很好的性能。表1為SuperMC和MCNP在不同核數(shù)下的并行效率對比，得出基于本文方法的SuperMC并行效率高于MCNP。

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

4 結(jié)語

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

致謝 感謝中國科學(xué)院核能安全技術(shù)研究所FDS團(tuán)隊(duì)其他成員對本文工作的支持幫助。

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

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

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

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

2015-01-29，

2015-03-20