付元光，劉 鵬，李 瑞，王 鑫，鄧 力

(1.中國工程物理研究院 高性能數值模擬軟件中心，北京 100088；2.北京應用物理與計算數學研究所，北京 100088)

蒙特卡羅臨界計算是核臨界安全分析的重要手段，由于蒙特卡羅方法具有幾何描述精確、物理過程細致等特性，其臨界計算結果常被視為參考解。使用蒙特卡羅方法模擬復雜問題時，需采用大規(guī)模并行計算和海量粒子樣本確保問題收斂。研究發(fā)現，當并行規(guī)模較大時，裂變源分配方法和粒子權重調節(jié)策略將顯著影響程序的負載平衡特征。本工作對比幾種粒子分配與調權算法的性能，提出一種基于本地權重守恒的調權算法，通過實現動態(tài)負載平衡，以獲得較高的加速比與并行效率。

1 蒙特卡羅臨界計算流程

在蒙特卡羅臨界計算中，通常采用裂變源迭代算法計算核系統(tǒng)的有效增殖因數。在該算法中，從某代初始裂變源分布抽樣源粒子進行輸運模擬，當粒子發(fā)生裂變時，記錄當前位置作為下一代裂變源點用于下一代跟蹤。系統(tǒng)的有效增殖因數keff可通過連續(xù)兩代之間的裂變源總權重比值求得，通過多代迭代計算，keff均值和裂變源將逐漸收斂。然而在實際計算機程序的實現中，對于keff>1，會產生越來越多的裂變位置，直至占滿內存；對于keff<1，會產生越來越少的裂變位置，直至為0。為保證每代模擬中子的總數和總權重守恒，防止無限增殖或消亡，需對粒子權重或裂變位置數進行糾偏[1-3]。

一般，串行臨界計算的實現步驟[4-5]如下。在開始前需給定keff的初始猜測值k0、每代模擬粒子數N、循環(huán)代數It等參數以及一個初始裂變源分布，接著：

1)權重為1的N個中子逐個從初始裂變源發(fā)出并開始隨機游動。單個中子每次發(fā)生碰撞時記錄n個當前碰撞位置為新的裂變位置：

(1)

2)粒子發(fā)生碰撞或游動，會對keff的估計值產生貢獻，一般有碰撞、吸收、徑跡長度估計3種，以碰撞估計為例，單個中子每次碰撞對估計值的貢獻為：

(2)

其中，fi為核素i的核子密度。

4)在新一代的模擬中，M個中子將從M個新的裂變位置發(fā)射，并重復1～3的步驟。需注意，新一代中子的出生權重問題，理論上M是N的無偏估計，但由于存在統(tǒng)計漲落，M會稍微偏離N，這時需調節(jié)粒子出生權重為N/M，以保證每代模擬粒子總權重守恒為N。

2 粒子并行算法研究

2.1 統(tǒng)一調節(jié)粒子權重的算法

實際發(fā)現這種算法存在嚴重的負載不均衡問題。使用了一個一維兩區(qū)域單能量各向同性散射模型進行數值實驗，模型參數如圖2所示。將初始源設置在x=25 cm處，使用4個進程、不同樣本數、循環(huán)200代展開模擬，隨著模擬代數的增加，各進程上的裂變位置數變化情況如圖3所示。可看出，裂變位置數在不同進程上的不均衡性隨著每代模擬樣本的減少而增加，在30 000粒子樣本時，不同進程勉強能維持在同一水平；在3 000粒子樣本時，不均衡性有所加??；在300粒子樣本時，0號進程不斷減小直到為0，1號進程不斷增大，3號進程即將變?yōu)?，不均衡已非常顯著。可推斷，當每代模擬樣本量越少時，計算不均衡問題越嚴重。

圖2 數值實驗模型

a——每代模擬30 000；b——每代模擬3 000；c——每代模擬300

2.2 Master-Slave算法

如上所述，在初始雖然給每個進程均勻分配了裂變位置數，但隨著模擬進行，本地裂變位置數出現漲落，使得負載不均衡出現。Master-Slave算法通過在每個循環(huán)代結束、下一代開始前進行裂變位置的收集和再分配，以達到平衡各進程間權重的效果，代表程序是MCNP。若使用m+1個進程并行模擬，則會產生1個Master進程和m個Slave進程，Master進程負責數據收集與分發(fā)，Slave進程進行粒子并行模擬任務。每代模擬完畢后，各Slave進程將本地裂變位置發(fā)送給Master，Master搜集所有信息后平均分配成m份，再發(fā)送到各Slave進程上，保證每代各Slave進程的模擬權重相等。圖4給出一個N=3 000、m=5的示例。Master-Slave算法由于每代循環(huán)都需進行全局收集和分發(fā)通信操作，通信耗時會隨著裂變位置數、循環(huán)代數的增大而顯著增長。

圖4 Master-Slave算法示例

2.3 Nearest Neighbor算法

為減少每代循環(huán)間的通信數據量，以OpenMC為代表的程序提出Nearest Neighbor算法[6]。該算法拋棄了Master進程，而在各Slave進程上點對點通信必要的數據。若使用m個進程并行模擬N個粒子，則每個進程初始分配到N/m個粒子，1代結束后，每個進程上產生的粒子不等于N/m，將新的粒子排序編號，對j號進程而言，設其本地所有的粒子中第1個和最后1個編號分別是aj、bj，若ajjN/m則從j-1號進程接收aj-jN/m個粒子；類似地，若bjj(N+1)/m，則將bj-j(N+1)/m個粒子發(fā)送給j+1號進程。圖5給出一個N=3 000、m=4的示例。Nearest Neighbor算法能保留一部分本地的裂變位置，避免了所有位置的收集分發(fā)操作，使用點對點通信，有效減少了通信數據量。

圖5 Nearest Neighbor算法示例

2.4 調節(jié)本地出生粒子權重的算法

圖6 改進的粒子并行的源迭代算法

圖7 不同進程上裂變位置數隨循環(huán)代的變化

本算法的優(yōu)點在于通信量為上述諸算法中最小，由式(2)可知，每代的keff仍需全局規(guī)約求和，除此之外本算法沒有其他需要通信的數據，而keff數據量較粒子的裂變位置信息數據量基本可忽略不計。本算法的不足之處在于，無法保證串并行一致，對比之前的幾種算法，可發(fā)現它們均能保證每個進程上的粒子出生權重完全相同，這樣只要保證每個粒子分配到的隨機數序列與串行完全相同，那么無論不同進程間如何分配粒子，都會得到與串行相同的結果，而這對于常用的乘加同余隨機數發(fā)生器是易實現的；對于本算法，每個進程在每次循環(huán)代的出生權重各不相同，結合keff估計的累加式(2)，可發(fā)現權重的差異影響keff的計算結果，即使每個粒子使用與串行完全相同的隨機數序列，也無法還原串行下每個粒子累加時的權重，進而造成串并行的不一致。盡管如此，對蒙特卡羅隨機模擬而言，只要計算結果可收斂，即使串并行不一致也可接受。

3 數值結果

3.1 壓水堆單柵元測試

用1個帶有全反射邊界條件的壓水堆棒柵算例[8]，模擬3×106個粒子，采用3 000粒子×1 000代、30 000粒子×100代、300 000粒子×10代 3種模式測試上述算法的加速比，結果列于表1～4，其中計算耗時主要包括粒子的跟蹤與通信時間，去除了初始化、截面讀取等其他耗時。

表1 統(tǒng)一調節(jié)粒子權重算法并行測試結果

由表1可看出，對于統(tǒng)一調節(jié)粒子權重算法，加速比很差，可推測由于負載不均的存在，進程存在大量等待現象。由表2可看出，Master-Slave算法加速比整體有一定程度提升，每代粒子數越少、代數越多，加速比越差，對于3 000粒子×1 000代的情況，在120核時加速比甚至出現了下降，這主要因為代數越多，通信的次數越多，通信的調用越頻繁，額外開銷越大，加速效果下降越明顯。由表3可看出，Nearest Neighbor算法加速比又有進一步提升，對于少樣本、多代數仍存在加速比較低的現象，但較Master-Slave算法已改善很多。從表4可看出，本工作提出的新算法加速性能整體優(yōu)于Nearest Neighbor和Master-Slave。各算法的共同現象是少樣本、多代數設置下加速比最低，這是由于源迭代算法總是在每代之間進行通信，越多的代數意味著越多的通信操作，因此在蒙特卡羅臨界計算粒子并行模擬時，最好使用多樣本、少代數的設置。

表2 Master-Slave算法并行測試結果(單柵元)

表4 本地權重調節(jié)算法并行測試結果(單柵元)

為說明本工作算法在串并行不一致方面的影響，統(tǒng)計了3 000粒子×1 000代情況下使用不同進程數計算keff的結果，列于表5?？煽闯觯畲笾蹬c最小值的差距為62 pcm，說明算法雖無法保證串并行一致，但不影響收斂性。

表5 使用不同進程數計算keff結果對比

3.2 壓水堆2×2組件測試

參考C5G7-2D基準題[9]構造了一個2×2組件排列的模型，其中燃料棒尺寸、柵距、反射層等幾何參數以及邊界條件設計和C5G7-2D完全相同，略有差別的是所有棒束(包括導向管位置)均使用二氧化鈾。采用每代10 000粒子×1 000代、100 000粒子×100代、1 000 000粒子×10代 3種模式測試加速比，對于統(tǒng)一調節(jié)權重算法，由上節(jié)可看出其效率過差，因此沒有再進行測試，僅就Master-Slave、Nearest Neighbor和新算法進行測試，結果列于表6～8。

表6 Master-Slave算法并行測試結果(2×2組件)

表7 Nearest Neighbor算法并行測試結果(2×2組件)

表8 本地權重調節(jié)算法并行測試結果(2×2組件)

可看出，類似上節(jié)的測試，新算法整體上加速效果優(yōu)于前兩者；Master-Slave算法在少粒子、多代數設置下仍會出現加速比快速下降的現象；值得注意的是，Nearest Neighbor算法在多粒子、少代數設置下性能較Master-Slave差，原因分析如下：當進程數增大時，每個核分配的計算任務變少，計算時間縮短，通信時間所占比重增大，由于OpenMC的串行速度最快，隨著進程變多，其通信占比提升更明顯，上述測試的樣本量對于OpenMC略顯不足，增大樣本量可提升OpenMC性能。

3.3 全堆大規(guī)模測試

在天河2超算上使用BEAVRS模型[10]測試本工作算法的并行效率，BEAVRS是由MIT提出的壓水堆堆芯模擬驗證基準題，具有193個燃料組件和17×17棒束結構，對BEAVRS進行精細描述需百萬量級柵元，適用于較大規(guī)模的并行測試。BEAVRS具有兩個循環(huán)周期的實測數據，盡管計算這些內容需多物理耦合，但出于并行效率測試的目的，本工作僅對其初始冷態(tài)零功率工況進行蒙特卡羅輸運模擬。所有涉及到的測試時間均為去除數據初始化步驟(輸入解析，截面讀取、廣播，內存分配等)后余下的源迭代流程的執(zhí)行時間。

對BEAVRS進行了pin-by-pin精細建模，燃料棒軸向劃分16層，模型總柵元數達到540余萬。首先測試了弱可擴展并行效率，通過調節(jié)粒子樣本數以控制程序計算量，樣本數的設置、計算結果、并行效率等列于表9。然后測試了強可擴展并行效率，固定每代使用107粒子、循環(huán)500代，計算結果列于表10?？梢钥闯觯鄬τ?28進程，4 800進程的弱可擴效率達到92.54%，強可擴效率達到81.47%，說明算法具有較好的并行性能。keff最大值與最小值差距為35 pcm，說明算法不會影響計算的收斂。

表9 弱可擴展并行效率測試結果

表10 強可擴展并行效率測試結果

4 結論

在蒙特卡羅臨界計算粒子并行模擬中，隨機漲落的存在會顯著影響并行的負載均衡特性，如果不平衡每個進程上的粒子數，將會發(fā)生嚴重的負載不均衡現象，特別當每代模擬數目較小、代數較多時，負載均衡更差。Master-Slave算法通過粒子的全局收集與再分配，實現各進程上模擬權重的平衡，但會引入較多通信操作。Nearest Neighbor算法通過保留部分本地粒子、在相鄰結點間點對點通信部分粒子，有效降低了通信數據量。本工作提出一種基于本地權重期望值守恒調節(jié)出生粒子權重的方法，可以保證隨著循環(huán)代數增加，各進程上的裂變位置數和本地權重不出現顯著偏離，負載均衡較好，且數據通信操作更少?；趬核寻魱?、組件問題，通過不同粒子數和循環(huán)代數設置，測試了各算法的加速性能，本工作提出的算法較Master-Slave和Nearest Neighbor算法獲得進一步提升，特別對于少粒子、多循環(huán)的參數設置，算法提升最為明顯。在天河2超算上測試了算法的并行效率，BEAVRS冷態(tài)零功率計算4 800進程相對128進程的弱可擴和強可擴效率可以達到92.54%和81.47%以上。