馬黨偉，劉宙宇，趙 晨，賀 濤，張文鑫，曹良志，吳宏春，柴曉明

(1.中國(guó)核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

NECP-X[1-2]是由西安交通大學(xué)核工程計(jì)算物理實(shí)驗(yàn)室(NECP)自主開(kāi)發(fā)的三維全堆芯高保真輸運(yùn)計(jì)算程序系統(tǒng)，通過(guò)精確幾何處理、自主研發(fā)多群數(shù)據(jù)庫(kù)、全局-局部共振計(jì)算[3]、大規(guī)模并行二維/一維輸運(yùn)計(jì)算、雙重CMFD加速、pin-by-pin物理熱工耦合[2]及非均勻全堆瞬態(tài)計(jì)算[4]實(shí)現(xiàn)全堆芯高保真輸運(yùn)計(jì)算，獲得反應(yīng)堆全三維精細(xì)的穩(wěn)態(tài)和瞬態(tài)功率分布。其中輸運(yùn)計(jì)算為基于特征線(xiàn)方法(MOC)的二維/一維輸運(yùn)方法。

特征線(xiàn)方法具有強(qiáng)大的幾何處理能力，被廣泛應(yīng)用于反應(yīng)堆物理計(jì)算中。然而，當(dāng)利用特征線(xiàn)方法進(jìn)行三維全堆芯輸運(yùn)計(jì)算時(shí)，將面臨計(jì)算量極大、內(nèi)存要求極高的挑戰(zhàn)。結(jié)合當(dāng)前計(jì)算機(jī)硬件的發(fā)展，大規(guī)模并行計(jì)算成為重要的解決方法。目前發(fā)展中的全堆芯三維輸運(yùn)方法大部分是基于單一變量并行，如能群、空間、角度或特征線(xiàn)，其并行度均較有限，無(wú)法充分利用當(dāng)前最先進(jìn)的計(jì)算機(jī)資源，如韓國(guó)的CRX-3D[5]僅實(shí)現(xiàn)了角度的并行，DeCART[6]實(shí)現(xiàn)了軸向分層的并行，中國(guó)核動(dòng)力研究設(shè)計(jì)院的Tiger-3D程序[7]實(shí)現(xiàn)了空間的并行，西安交通大學(xué)的TOMMOC基于空間區(qū)域分解并行，AutoMOC-3D基于特征線(xiàn)并行等。本文針對(duì)NECP-X中二維/一維耦合輸運(yùn)計(jì)算方法實(shí)現(xiàn)基于空間、角度和特征線(xiàn)的多重并行，并在性能分析的基礎(chǔ)上對(duì)程序進(jìn)行效率優(yōu)化。

1 NECP-X中的并行策略

數(shù)值反應(yīng)堆物理計(jì)算程序NECP-X是基于二維/一維耦合方法進(jìn)行輸運(yùn)計(jì)算的，采用二維MOC處理徑向精細(xì)的幾何結(jié)構(gòu)，采用一維SN方法進(jìn)行軸向計(jì)算?？紤]到單純的MOC源迭代收斂速度較慢且計(jì)算量較大，NECP-X采用了粗網(wǎng)有限差分(CMFD)[1]和多重并行的方法加速收斂或計(jì)算。其中，多重并行不僅可提高計(jì)算速度，還可分解內(nèi)存壓力，實(shí)現(xiàn)大規(guī)模問(wèn)題的求解。NECP-X中的多重并行包括空間、角度和特征線(xiàn)的并行，多重并行的并行方案如圖1所示。

1.1 空間區(qū)域分解及效率優(yōu)化

NECP-X空間區(qū)域分解是基于分布式存儲(chǔ)的并行編程技術(shù)MPI[8]進(jìn)行設(shè)計(jì)的，這對(duì)于緩解全堆芯問(wèn)題的內(nèi)存壓力非常重要。NECP-X的空間區(qū)域分解通過(guò)軸向和徑向兩次分解實(shí)現(xiàn)，如圖2所示。對(duì)于軸向區(qū)域分解，根據(jù)軸向?qū)訑?shù)和擬分軸向區(qū)域數(shù)，獲得每個(gè)計(jì)算核心的層數(shù)，若無(wú)法平均分配，則部分計(jì)算核心多計(jì)算一層。為了負(fù)載平衡，建議軸向區(qū)域數(shù)為軸向?qū)訑?shù)的約數(shù)。對(duì)于徑向區(qū)域分解，提供了兩種分解方式：1) 根據(jù)徑向的模塊數(shù)和區(qū)域數(shù)計(jì)算每個(gè)區(qū)域的模塊數(shù)；2) 由用戶(hù)人為指定區(qū)域分解的方式，這種分解方式具有較大的靈活性，但也給用戶(hù)的使用帶來(lái)了難度。用戶(hù)可根據(jù)問(wèn)題選擇合適的分解方式。軸向和徑向兩次分解確定了所有的三維子區(qū)域內(nèi)的軸向?qū)訑?shù)和徑向模塊數(shù)。

圖1 NECP-X的多重并行方案Fig.1 Multilevel parallel scheme of NECP-X

圖2 空間區(qū)域分解示意圖Fig.2 Space decomposition scheme

三維空間區(qū)域分解后，每個(gè)計(jì)算核心負(fù)責(zé)同一區(qū)域內(nèi)的二維MOC和一維SN計(jì)算。此外，每個(gè)計(jì)算核心僅與相鄰區(qū)域通信，因此每個(gè)計(jì)算核心的通信量不隨計(jì)算核心數(shù)線(xiàn)性增長(zhǎng)，顯著提高了并行擴(kuò)展性；并通過(guò)非阻塞通信減少二維/一維迭代計(jì)算過(guò)程中的通信等待。所以，空間并行理論上可獲得很高的并行效率和極大的并行度。

空間區(qū)域分解后，各子空間區(qū)域同時(shí)求解，從根本上改變了迭代流程，因此與串行計(jì)算相比，收斂過(guò)程也發(fā)生了改變。在空間區(qū)域分解中，各子區(qū)域在源迭代計(jì)算中的初始值均相同，只有在更新內(nèi)邊界條件時(shí)，各子區(qū)才會(huì)進(jìn)行信息交換。而在串行計(jì)算中，對(duì)于二維MOC計(jì)算，當(dāng)前模塊的初始通量值是由相鄰的已計(jì)算模塊的出射通量值傳遞而來(lái)的，即串行計(jì)算時(shí)通量更新較為及時(shí)。故空間區(qū)域分解的迭代格式相對(duì)于串行計(jì)算會(huì)有退化。

1.2 角度并行和特征線(xiàn)并行簡(jiǎn)介

角度區(qū)域分解是針對(duì)特征線(xiàn)方法和SN方法中的角度變量實(shí)施的并行策略，可在整個(gè)二維/一維輸運(yùn)計(jì)算中對(duì)角度進(jìn)行并行，實(shí)現(xiàn)一維計(jì)算和二維計(jì)算的并行度匹配。由于各角度的特征線(xiàn)段數(shù)不同，采用貪婪算法[9]的分配策略實(shí)現(xiàn)更好的角度區(qū)域分解負(fù)載平衡，如圖3所示。圖中不同矩形的面積代表不同角度的特征線(xiàn)段數(shù)量。

圖3 角度區(qū)域分解的貪婪算法Fig.3 Greedy algorithm for angle decomposition

特征線(xiàn)并行是對(duì)MOC計(jì)算的核心計(jì)算部分(特征線(xiàn)掃描部分)進(jìn)行的并行加速處理，通過(guò)基于共享內(nèi)存的編譯處理方案OpenMP[10]實(shí)現(xiàn)。在單純的特征線(xiàn)方法中，特征線(xiàn)的掃描部分會(huì)占到總計(jì)算時(shí)間的90%以上[11]，所以針對(duì)特征線(xiàn)掃描進(jìn)行加速很有必要。角度并行和特征線(xiàn)并行的實(shí)現(xiàn)詳見(jiàn)文獻(xiàn)[11]。

1.3 特征線(xiàn)并行的效率分析及優(yōu)化

三維輸運(yùn)計(jì)算時(shí)，程序早期的數(shù)值結(jié)果表明特征線(xiàn)并行的并行效率隨線(xiàn)程數(shù)的增加下降較快，如圖4所示，10個(gè)線(xiàn)程并行時(shí)并行效率已低于60%。為此本文針對(duì)特征線(xiàn)并行部分進(jìn)行分析與優(yōu)化。

圖4 特征線(xiàn)并行效率Fig.4 Efficiency of characteristic ray parallel

特征線(xiàn)掃描部分的偽代碼如下：

1:do i_sweep=1,n_sweep #內(nèi)迭代循環(huán)開(kāi)始

2:do iang=stt_ang,stp_ang #幅角循環(huán)開(kāi)始

3:calculate source item #源項(xiàng)計(jì)算

4:do iray=1,n_ray #特征線(xiàn)掃描計(jì)算循環(huán)開(kāi)始

5:calculate fluxes #特征線(xiàn)掃描計(jì)算

6:end do #特征線(xiàn)掃描計(jì)算循環(huán)結(jié)束

7:tally angular flux #統(tǒng)計(jì)角通量

8:tally current & leakage #流與泄漏項(xiàng)計(jì)算

9:calculate scalar flux #標(biāo)通量計(jì)算，線(xiàn)程歸約

10:end do #幅角循環(huán)結(jié)束

11:update boundary condition #更新邊界條件

12:end do #內(nèi)迭代循環(huán)結(jié)束

圖5為串行計(jì)算特征線(xiàn)掃描各部分的時(shí)間比例，其中藍(lán)色為特征線(xiàn)掃描計(jì)算，紅色為源項(xiàng)計(jì)算。可明顯看出，特征線(xiàn)掃描計(jì)算約占總時(shí)間的90%以上。特征線(xiàn)并行通過(guò)在該部分添加OpenMP的編譯制導(dǎo)語(yǔ)句實(shí)現(xiàn)多線(xiàn)程的動(dòng)態(tài)任務(wù)分配。

圖5 串行計(jì)算的時(shí)間分布Fig.5 Time distribution of serial calculation

通過(guò)敏感性分析，并行起始區(qū)間為內(nèi)迭代循環(huán)開(kāi)始處。在此條件下，研究特征線(xiàn)并行的負(fù)載平衡、通信開(kāi)銷(xiāo)等問(wèn)題。本研究利用TAU(tuning and analysis utilities)[12]和PAPI(performance application programming interface)庫(kù)[13]統(tǒng)計(jì)并行化后的特征線(xiàn)掃描部分每段的計(jì)算時(shí)間，特征線(xiàn)并行的弱可擴(kuò)展性[14]測(cè)試表明，隨線(xiàn)程數(shù)的增加，特征線(xiàn)掃描計(jì)算的并行效率一直保持在94%以上，由于采用了動(dòng)態(tài)調(diào)度的分配方案，10個(gè)線(xiàn)程依然保持著良好的負(fù)載平衡，如圖6所示，圖中左側(cè)數(shù)據(jù)為每個(gè)線(xiàn)程的特征線(xiàn)掃描計(jì)算時(shí)間。然而，非并行部分(如源項(xiàng)計(jì)算、標(biāo)通量計(jì)算等)計(jì)算時(shí)間理論上應(yīng)隨線(xiàn)程數(shù)線(xiàn)性增加，測(cè)試結(jié)果(表1)卻顯示為超線(xiàn)性增加，因而造成了整體并行效率的急劇下降。

圖6 10個(gè)線(xiàn)程時(shí)特征線(xiàn)掃描時(shí)間Fig.6 Time of characteristic ray sweeping at 10 threads

表1 不同線(xiàn)程時(shí)特征線(xiàn)并行時(shí)間及效率Table 1 Characteristic ray parallel time and efficiency for various thread numbers

結(jié)合特征線(xiàn)理論模型分析程序源代碼后，通過(guò)消除源項(xiàng)計(jì)算中的重復(fù)性計(jì)算以及將原來(lái)的串行計(jì)算改為多線(xiàn)程并行計(jì)算，同時(shí)對(duì)并行區(qū)域內(nèi)的其他部分通過(guò)優(yōu)化數(shù)據(jù)結(jié)構(gòu)提高緩存命中率，最終實(shí)現(xiàn)了并行效率的明顯提升。優(yōu)化后的特征線(xiàn)并行時(shí)間及效率列于表2?？煽闯?，源項(xiàng)計(jì)算時(shí)間由53.59 s降低至11.52 s，時(shí)間占比降到10%以下。標(biāo)通量的計(jì)算時(shí)間隨線(xiàn)程數(shù)的增加也出現(xiàn)非線(xiàn)性增加，該部分的源代碼為：

do ireg=stt_fsr_1+1，stt_fsr

tphi(：，ireg，ithd)=tphi(：，ireg，ithd)+

phi_bar(：，ireg)*weights(iang)

enddo

其中：ithd表示每個(gè)線(xiàn)程的編號(hào)；tphi表示標(biāo)通量。TAU的測(cè)試結(jié)果表明，該部分每個(gè)線(xiàn)程的浮點(diǎn)運(yùn)算次數(shù)不隨線(xiàn)程數(shù)而改變，其理論計(jì)算時(shí)間應(yīng)保持不變。但時(shí)間的測(cè)試結(jié)果卻為超線(xiàn)性增加，與理論不符。由于該部分的計(jì)算時(shí)間占比較小，且多種測(cè)試分析均未獲得該現(xiàn)象的原因，目前未進(jìn)行進(jìn)一步的優(yōu)化。

表2 優(yōu)化后的特征線(xiàn)并行時(shí)間及效率Table 2 Characteristic ray parallel time and efficiency after optimization

2 數(shù)值驗(yàn)證

本文首先分別驗(yàn)證NECP-X的空間、角度及特征線(xiàn)并行在小規(guī)模問(wèn)題下的并行效率，然后分別以真實(shí)三維單組件問(wèn)題和大規(guī)模并行可擴(kuò)展性分析為例說(shuō)明NECP-X多重并行的正確性和并行能力。

分別對(duì)空間、角度和特征線(xiàn)的并行效率進(jìn)行測(cè)試。為使每個(gè)核的計(jì)算量盡可能地相等，即達(dá)到較好的負(fù)載平衡，本測(cè)試選取了4×4×4的三維壓水堆組件，整個(gè)組件可分為64個(gè)柵元模塊，如圖2所示，每個(gè)柵元分為320個(gè)平源區(qū)。為使測(cè)試結(jié)果真實(shí)可靠，每個(gè)算例測(cè)試進(jìn)行100次特征值迭代計(jì)算。計(jì)算條件為在每個(gè)卦限內(nèi)8個(gè)幅角、3個(gè)極角，特征線(xiàn)寬度為0.03 cm，使用7群宏觀截面。其中，特征線(xiàn)并行的測(cè)試?yán)}為上述三維組件的四分之一，且使用69群宏觀截面，其余參數(shù)保持一致。測(cè)試算例所使用計(jì)算機(jī)的相關(guān)參數(shù)如下：處理器為Intel(R) Xeon(R) CPU E5-2640 v4@2.40 GHz(79)、計(jì)算核心個(gè)數(shù)為20、Ubuntu14.0.4操作系統(tǒng)、編譯器為GNU Fortran 4.8.4。

2.1 空間并行的數(shù)值驗(yàn)證

受限于計(jì)算機(jī)計(jì)算核心的個(gè)數(shù)，對(duì)于空間并行效率的測(cè)試分徑向和軸向區(qū)域分解分別進(jìn)行。為保證負(fù)載平衡，在徑向區(qū)域分解時(shí)分別使用1、2、4、8、16個(gè)核進(jìn)行計(jì)算，在軸向區(qū)域分解時(shí)分別使用1、2、4個(gè)核進(jìn)行計(jì)算。在保證計(jì)算核個(gè)數(shù)不超過(guò)20的情況下，對(duì)徑向和軸向同時(shí)區(qū)域分解進(jìn)行測(cè)試。3種情況下的測(cè)試結(jié)果大體相近，其中徑向區(qū)域分解的并行加速比示于圖7。

圖7 NECP-X空間并行加速比Fig.7 Space parallel speedup of NECP-X

從圖7可看出，空間并行在使用4個(gè)核時(shí)，總時(shí)間、一維計(jì)算時(shí)間和二維計(jì)算時(shí)間的加速比均接近于4；16個(gè)核并行時(shí)總時(shí)間的加速比也在15左右，一維計(jì)算時(shí)間加速比下降較快，是因?yàn)橐痪S計(jì)算時(shí)間在單核時(shí)也僅為0.29 s，16個(gè)核時(shí)已下降到0.03～0.02 s，存在較大的測(cè)量不確定性，不具有參考性。

2.2 角度并行的數(shù)值驗(yàn)證

測(cè)試?yán)}為在每個(gè)卦限內(nèi)8個(gè)幅角，故在0～180°范圍內(nèi)共有16個(gè)幅角。為使分配到每個(gè)核的計(jì)算任務(wù)盡量相同，在對(duì)角度并行進(jìn)行測(cè)試時(shí)分別使用了1、2、4、8、16個(gè)核進(jìn)行計(jì)算，測(cè)試結(jié)果如圖8所示，由于一維計(jì)算時(shí)間較少，不具有參考性，故未在圖中示出。從前文介紹的貪婪算法可知，當(dāng)使用2個(gè)核并行時(shí)，2個(gè)核的計(jì)算量是完全相同的，此時(shí)的并行加速比約為2，也證實(shí)了貪婪算法的高效；隨著核數(shù)的增加，每個(gè)核的計(jì)算量逐漸出現(xiàn)不同，當(dāng)使用16個(gè)核并行時(shí)，分配到每個(gè)核的特征線(xiàn)段數(shù)目存在較大不同，同時(shí)通信時(shí)間增加，此時(shí)并行加速比有所下降，但也在14左右。

圖8 NECP-X角度并行加速比Fig.8 Angle parallel speedup of NECP-X

2.3 特征線(xiàn)并行的數(shù)值驗(yàn)證

在驗(yàn)證特征線(xiàn)并行的加速效果時(shí)，本文采用69群的宏觀截面，因?yàn)樵趯?shí)際計(jì)算中，大部分問(wèn)題都使用了69群甚至更多能群的宏觀截面，這樣特征線(xiàn)掃描部分的計(jì)算時(shí)間占比也更趨于真實(shí)。

本文選擇1個(gè)線(xiàn)程到10個(gè)線(xiàn)程分別測(cè)試其加速比，測(cè)試結(jié)果如圖9所示?？煽闯?，當(dāng)線(xiàn)程數(shù)達(dá)到4時(shí)，并行加速比依然在3.5以上，并行效率接近90%；隨著線(xiàn)程數(shù)的逐漸增加，非并行部分計(jì)算時(shí)間逐漸增加，同時(shí)各線(xiàn)程對(duì)共有數(shù)據(jù)的讀取存在競(jìng)爭(zhēng)關(guān)系，造成并行加速比逐漸下降。在實(shí)際使用中，特征線(xiàn)的并行度一般不超過(guò)4個(gè)線(xiàn)程。

圖9 NECP-X特征線(xiàn)并行加速比Fig.9 Characteristic ray parallel speedup of NECP-X

2.4 多重并行的數(shù)值驗(yàn)證

NECP-X多重并行能力的測(cè)試?yán)}為美國(guó)橡樹(shù)嶺國(guó)家實(shí)驗(yàn)室在2014年發(fā)布的VERA(virtual environment for reactor applications)基準(zhǔn)題[14]中的3A問(wèn)題。該問(wèn)題為一典型的壓水堆17×17三維燃料組件在壽期初熱態(tài)零功率工況下的穩(wěn)態(tài)計(jì)算，詳細(xì)的問(wèn)題描述可參考文獻(xiàn)[14]。該三維組件在軸向上共有8個(gè)格架，同時(shí)考慮到軸向的上下端塞以及上下反射層區(qū)域等，在軸向上被劃分為66層，其模型示意圖如圖10所示。

圖10 VERA問(wèn)題3A參考模型Fig.10 VERA problem 3A reference model

本次驗(yàn)證NECP-X多重并行使用的計(jì)算機(jī)是“天河二號(hào)”超級(jí)計(jì)算機(jī)，每次測(cè)試均為節(jié)點(diǎn)獨(dú)占式。計(jì)算條件為每個(gè)卦限內(nèi)使用8個(gè)幅角、3個(gè)極角，特征線(xiàn)寬度為0.03 cm。為保證負(fù)載平衡，本文選擇空間軸向33個(gè)區(qū)域、角度2個(gè)區(qū)域、特征線(xiàn)2個(gè)線(xiàn)程并行的并行測(cè)試組合，測(cè)試結(jié)果列于表3?；鶞?zhǔn)題3A特征值的參考值為1.175 72，NECP-X在p0散射的計(jì)算條件下的計(jì)算結(jié)果為1.176 36，相差64 pcm，并行計(jì)算結(jié)果與串行計(jì)算結(jié)果相差不到10 pcm，說(shuō)明多重并行的計(jì)算結(jié)果正確。需要說(shuō)明的是，本次測(cè)試為強(qiáng)可擴(kuò)展性[11]效率測(cè)試，由于上管座等位置處與燃料段的網(wǎng)格劃分不一致，所以負(fù)載并不平衡，將會(huì)影響并行效率。另外，測(cè)試使用了粗網(wǎng)有限差分(CMFD)加速，但CMFD目前并未使用角度和特征線(xiàn)的并行資源，所以在使用角度和特征線(xiàn)并行時(shí)并行效率會(huì)有明顯下降，后期將對(duì)CMFD如何高效利用MOC中角度和特征線(xiàn)的并行資源進(jìn)行研究，這將顯著提升總并行效率。

表3 三維組件的多重并行計(jì)算Table 3 Multilevel parallel calculation of 3-D assembly

2.5 大規(guī)模并行可擴(kuò)展性分析

圖11 NECP-X空間并行弱可擴(kuò)展性測(cè)試Fig.11 Space parallel weak scaling test of NECP-X

對(duì)NECP-X的大規(guī)模(上萬(wàn)核)并行能力進(jìn)行了弱可擴(kuò)展性[15]分析。該分析采用的測(cè)試?yán)}為徑向每層225(15×15)個(gè)17×17組件，采用四分之一組件模塊化，故每層有225個(gè)模塊；軸向上隨著問(wèn)題規(guī)模的增大而增加層數(shù)，故總的模塊數(shù)為225的整數(shù)倍。測(cè)試計(jì)算機(jī)為“天河二號(hào)”超級(jí)計(jì)算機(jī)，計(jì)算條件與2.4節(jié)相同，使用69群宏觀截面。在每次測(cè)試中保證每個(gè)核計(jì)算1個(gè)模塊，輸運(yùn)計(jì)算迭代次數(shù)為30次，不使用CMFD，并計(jì)算了弱可擴(kuò)展性的加速比，結(jié)果如圖11所示。在10 125核大規(guī)模并行的情況下，并行加速比高達(dá)7 000以上，并行效率約為70%，體現(xiàn)了NECP-X較好的并行可擴(kuò)展性。

3 結(jié)論

本文介紹了數(shù)值反應(yīng)堆物理計(jì)算程序NECP-X的空間并行、角度并行和特征線(xiàn)并行的多重并行策略?；趨^(qū)域分解的空間并行可分解計(jì)算機(jī)的內(nèi)存壓力，極大地?cái)U(kuò)大了可計(jì)算問(wèn)題的規(guī)模，同時(shí)具有較高的并行效率；基于貪婪算法的角度并行可充分地考慮負(fù)載平衡，達(dá)到較高的加速效果；基于OpenMP的特征線(xiàn)并行在動(dòng)態(tài)調(diào)度的分配方案下可較好地加速計(jì)算，同時(shí)優(yōu)化了非并行區(qū)域的計(jì)算。計(jì)算結(jié)果表明，角度和特征線(xiàn)并行可在空間并行的基礎(chǔ)上進(jìn)一步擴(kuò)大并行規(guī)模，加速計(jì)算；空間并行可利用上萬(wàn)核的并行資源實(shí)現(xiàn)高效的三維全堆芯計(jì)算。NECP-X的多重并行策略可充分利用計(jì)算資源，實(shí)現(xiàn)不同規(guī)模問(wèn)題的高效并行計(jì)算。