先進(jìn)中子學(xué)柵格程序KYLIN-Ⅱ輸運(yùn)模塊并行優(yōu)化開發(fā)

涂曉蘭，柴曉明，劉 東，蘆 韡，王 鑫，李勛昭，付元光，郭鳳晨

(1.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點(diǎn)實驗室，四川 成都 610213；2.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100094)

中國核動力研究設(shè)計院不僅承擔(dān)了壓水堆的設(shè)計研發(fā)，同時承擔(dān)大量新型反應(yīng)堆的研發(fā)，不同反應(yīng)堆的燃料組件的結(jié)構(gòu)、燃料成分復(fù)雜，為能滿足設(shè)計研發(fā)的需求，中國核動力研究設(shè)計院開發(fā)形成了先進(jìn)中子學(xué)柵格程序KYLIN-Ⅱ[1-3]，KYLIN-Ⅱ采用子群方法[4]求解共振核素的共振截面，采用特征線方法[5-11]進(jìn)行復(fù)雜幾何中子輸運(yùn)計算，并采用廣義粗網(wǎng)格有限差分加速方法(GCMFD)[12]來加速中子輸運(yùn)求解流程，采用基于改進(jìn)預(yù)估修正臨界-燃耗迭代方法(PPC)的切比雪夫方法[13]求解復(fù)雜燃耗鏈，同時，為方便用戶使用，開發(fā)形成了支持復(fù)雜組件的圖形化建模工具[14]和后處理顯示工具。針對1個典型的二維AFA3G組件，其網(wǎng)格規(guī)模為4 000，能群為190群，特征線模塊計算時間可達(dá)7 000 s，計算時間較長，需對輸運(yùn)模塊進(jìn)行并行優(yōu)化，提升組件計算的效率。本文首先介紹輸運(yùn)計算模型，然后對KYLIN-Ⅱ程序的特征線輸運(yùn)模塊進(jìn)行效率分析和并行特征分析，針對熱點(diǎn)計算模塊，提出相應(yīng)的并行策略及優(yōu)化策略，最后通過基準(zhǔn)題和典型的壓水堆問題驗證本文方法的正確性和并行優(yōu)化的效率。

1 二維幾何下特征線法求解輸運(yùn)方程

真實情況下中子輸運(yùn)問題是三維問題，對于軸向均勻的問題，可將計算問題退化為二維問題，然而在求解二維問題時，必須保證計算結(jié)果與三維計算結(jié)果等效。三維情況下特征線方法及求解公式見文獻(xiàn)[5]。

二維情況下沿射線的通量ψ求解公式為：

(1)

其中：g為能群編號；i為網(wǎng)格編號；m為方位角編號；n為極角編號；k為特征線編號；s′n為射線在二維情況下的長度；θn為極角；Q為源項；Σt為總截面。

網(wǎng)格i的標(biāo)通量計算公式為：

ψi,m,n,k(si,m,n,k))Δdi,m,k

(2)

其中：ωn為極角權(quán)重；ωm為方位角權(quán)重；Σt,i為總截面；Ai為網(wǎng)格面積；Δdi,m,k為射線段的寬度。

在求解中子輸運(yùn)問題時，各向異性散射會對計算結(jié)果有較大的影響[15]，KYLIN-Ⅱ程序截面庫包含了高階Pn散射截面，在特征線中子輸運(yùn)模塊中，考慮了各向異性散射，故式(2)中包含了各向異性散射源項，KYLIN-Ⅱ程序?qū)ι⑸湓错椫械纳⑸浣孛嬷械慕嵌葢?yīng)用勒讓德多項式進(jìn)行展開，展開后的散射源項公式為：

Qscat,g,i,k=

(3)

其中：Σs,n,i,g′→g為散射截面；Pl為勒讓德函數(shù)；Ψi,g′為角通量；Ω為角度。

射線布置及追蹤方式影響特征線方法處理復(fù)雜幾何問題的能力，KYLIN-Ⅱ程序綜合應(yīng)用了循環(huán)射線布置及追蹤技術(shù)以及反向延長追蹤技術(shù)，循環(huán)射線布置及追蹤技術(shù)用于處理特定的問題，計算效率高；反向延長追蹤技術(shù)用于處理任意幾何問題，幾何適應(yīng)性強(qiáng)。

2 特征線輸運(yùn)模塊并行優(yōu)化方案設(shè)計

針對圖1所示的典型壓水堆二維AFA3G組件例題，對KYLIN-Ⅱ程序的輸運(yùn)模塊進(jìn)行效率分析，KYLIN-Ⅱ程序特征線輸運(yùn)模塊中特征線掃描、高階源項更新較耗時，因此本文主要針對特征線掃描、高階源項更新進(jìn)行優(yōu)化。

圖1 AFA3G組件幾何網(wǎng)格劃分Fig.1 Mesh division of AFA3G assembly

2.1 高階散射源計算優(yōu)化策略

KYLIN-Ⅱ程序應(yīng)用式(3)計算高階散射源，從公式可知能群g、方向k的散射源來自其他所有能群g′和方向k′角通量的貢獻(xiàn)，將是一個關(guān)于g、k、g′、k′的4重循環(huán)，循環(huán)次數(shù)較多，使得計算效率較低，同時必須存儲角通量數(shù)據(jù)以備計算，但會使用大量內(nèi)存，如網(wǎng)格數(shù)10 000、能群數(shù)200、方向數(shù)100，用雙精度存儲角通量會耗去約1.5 GB內(nèi)存；且在能群并行中，需通信角通量，而由于數(shù)據(jù)量大，會顯著增加通信時間，降低整體并行效率，因此本文考慮將角通量應(yīng)用球諧函數(shù)展開。

角通量應(yīng)用球諧函數(shù)展開的公式為：

(4)

(5)

式(3)中的勒讓德函數(shù)應(yīng)用球諧函數(shù)表示：

(6)

將式(4)和(6)代入式(3)，利用球諧函數(shù)的正交性，可得：

(7)

利用正交性可得通量矩的計算公式：

(8)

因此在特征線掃描的過程中不需存儲角通量，僅需存儲通量矩即可，與角通量相比，通量矩的散射階遠(yuǎn)小于角度離散數(shù)，對于L階散射，計算共需通量矩的個數(shù)為(L+1)2，一般情況下，(L+1)2較K小1～2個量級，可節(jié)省內(nèi)存，同時減少散射源計算的循環(huán)次數(shù)。

cosj(φk-φk′)

(9)

可見結(jié)果一定是實數(shù)，沒有虛部。因此構(gòu)造一個函數(shù)Zlj，滿足：

(10)

則有：

(11)

通過構(gòu)造Zlj函數(shù)，在計算過程中不需計算和存儲復(fù)數(shù)。

2.2 基于MPI的能群并行策略

由于輸運(yùn)模塊完全通過外迭代驅(qū)動收斂，即取消了散射源內(nèi)迭代，所有能群1次特征線掃描結(jié)束后同時更新散射源和裂變源，因此能群之間沒有數(shù)據(jù)依賴，且各能群之間的計算量嚴(yán)格相等，可采用能群并行，并建立能群并行通信域，同時在高階散射源更新過程中也涉及能群的循環(huán)，因此在特征線掃描和高階散射源更新部分均應(yīng)用MPI進(jìn)行能群的并行，并行流程如圖2所示，首先將宏觀截面等數(shù)據(jù)從主CPU廣播到其他CPU；然后進(jìn)行任務(wù)分配，當(dāng)核數(shù)等于能群數(shù)或核數(shù)被能群數(shù)整除時，將能群均分到核上，達(dá)到理想的負(fù)載均衡；特征線掃描完成后對網(wǎng)格標(biāo)通量、通量矩以及粗網(wǎng)凈中子流進(jìn)行數(shù)據(jù)通信；GCMFD求解過程中調(diào)用PETSC函數(shù)庫進(jìn)行并行求解，求解結(jié)束后需對粗網(wǎng)格通量和粗網(wǎng)格源項進(jìn)行數(shù)據(jù)通信。

圖2 特征線模塊計算流程圖Fig.2 Calculation flow of method of characteristics module

對特征線掃描部分進(jìn)行性能分析，E指數(shù)(式(1))計算時間占比較高，約占特征線掃描時間的40%，因此本文對E指數(shù)進(jìn)行優(yōu)化，即在第1次計算時，將E指數(shù)進(jìn)行存儲，由于在特征線掃描過程中應(yīng)用MPI并行模式對能群并行，E指數(shù)也將分能群存儲，即每個CPU只存儲當(dāng)前分配能群的E指數(shù)，因此并行計算時不會有太大的內(nèi)存消耗。

3 數(shù)值驗證

3.1 C5G7 UO2組件基準(zhǔn)題

C5G7 UO2組件基準(zhǔn)題是從C5G7 2D基準(zhǔn)題中提取的UO2組件，幾何布置如圖3所示，它的群常數(shù)參見OECD發(fā)布的基準(zhǔn)題文檔。計算時的網(wǎng)格規(guī)模為6 936，能群數(shù)為7，臨界計算keff收斂準(zhǔn)則是10-5，通量收斂準(zhǔn)則是10-4，特征線方法計算條件為射線間距為0.01，0～2π區(qū)間方位角數(shù)為40，0～π/2區(qū)間方位角滿足式(12)條件進(jìn)行布置，0～π區(qū)間6個高斯勒讓德極角。

m=1,2,…,N

φm∈(0,π/2)

(12)

其中：a為問題長度；b為問題寬度；N為0～π/2區(qū)間方位角數(shù)目。

圖3 C5G7基準(zhǔn)題的UO2組件Fig.3 Geometry of C5G7 UO2 assembly

表1所列為并行優(yōu)化后的計算結(jié)果，由于其截面不含高階散射，因此此基準(zhǔn)題角通量球諧函數(shù)展開不起作用，其kinf計算結(jié)果小數(shù)點(diǎn)后前12位數(shù)與串行計算結(jié)果保持一致，說明本文的E指數(shù)優(yōu)化和能群并行計算精度很高。應(yīng)用E指數(shù)優(yōu)化具有較高的計算效率，將特征線掃描速度提升75%。并行核數(shù)為7，加速比為6.37，并行的效率為90%，加E指數(shù)優(yōu)化后的并行效率提升更好，較大減少了特征線掃描的時間，對閉循環(huán)特征線的掃描過程效率提升明顯。

表1 C5G7 UO2組件計算結(jié)果Table 1 Calculation result of C5G7 UO2 assembly

表1列出的并行前、后迭代次數(shù)保持一致，由于輸運(yùn)模塊計算過程中取消了散射源內(nèi)迭代，所有能群一次特征線掃描結(jié)束后同時更新散射源和裂變源，能群之間沒有數(shù)據(jù)依賴，因此能群的并行并不會導(dǎo)致迭代的退化。

3.2 AFA3G燃料組件基準(zhǔn)題

針對AFA3G燃料組件進(jìn)行正確性測試和性能測試，AFA3G燃料組件含有呈17×17方形排列的289個柵元，包括264根燃料棒、24根鋯合金導(dǎo)向管和1根鋯合金測量管。燃料棒柵距為1.26 cm，燃料棒芯塊半徑為0.409 6 cm，燃料包殼外半徑為0.475 cm；導(dǎo)向管內(nèi)半徑為0.560 5 cm，外半徑為0.622 5 cm；中心導(dǎo)向管內(nèi)半徑為0.572 5 cm，外半徑為0.622 5 cm。本算例僅考慮熱態(tài)、不含可燃毒物的情況下，燃料富集度為3.0%時的計算，導(dǎo)向管和測量管內(nèi)部全部填充輕水慢化劑。

計算時的網(wǎng)格規(guī)模為4 105，能群數(shù)為190，臨界計算keff收斂準(zhǔn)則是10-5，通量收斂準(zhǔn)則是10-4，并行計算的核數(shù)為190，特征線方法計算條件為射線間距為0.01，0～2π區(qū)間方位角數(shù)目為40，0～π/2區(qū)間方位角滿足式(12)條件進(jìn)行布置，0～π區(qū)間6個高斯勒讓德極角。

表2列出了串行優(yōu)化計算結(jié)果，kinf計算結(jié)果基本一致，角通量球諧函數(shù)展開對高階散射源更新效率提升非常好，將高階散射源速度提升87%，E指數(shù)優(yōu)化提升了特征線掃描的效率，將特征線掃描速度提升77%，E指數(shù)優(yōu)化和角通量球諧函數(shù)展開同時應(yīng)用，整體計算時間效率得到了較好的提升，特征線模塊的整體速度提升83%。

表3列出了并行計算結(jié)果，kinf計算結(jié)果基本一致，并行前、后迭代次數(shù)保持一致，若只是對能群并行計算效率并不高，加E指數(shù)后特征線掃描并行效率有了更好的提升，加角通量球諧函數(shù)展開對散射源更新的并行效率有了更大的提升，綜合應(yīng)用E指數(shù)優(yōu)化和球諧函數(shù)展開后的并行效果較好，使原計算時間7 148.76 s縮短為38.85 s。

表2 AFA3G燃料組件串行優(yōu)化計算結(jié)果Table 2 Serial optimization result of AFA3G fuel assembly

表3 AFA3G燃料組件并行優(yōu)化計算結(jié)果Table 3 Parallel optimization result of AFA3G fuel assembly

4 結(jié)論

本文針對KYLIN-Ⅱ程序輸運(yùn)模塊中的特征線掃描和高階散射源計算應(yīng)用MPI進(jìn)行了能群并行，應(yīng)用角通量球諧函數(shù)展開方法對高階散射源計算進(jìn)行了優(yōu)化，應(yīng)用通量矩取代了角通量，提高了并行通信效率以及高階散射源的計算效率，并針對特征線掃描計算的E指數(shù)進(jìn)行了優(yōu)化存儲，提高了特征線計算效率，由于每個CPU只存儲當(dāng)前分配能群的E指數(shù)，因此并行計算時不會有太大的內(nèi)存消耗。并對多個基準(zhǔn)題進(jìn)行了驗證，結(jié)果表明，E指數(shù)優(yōu)化對特征線掃描加速效果顯著，角通量球諧函數(shù)展開對高階散射源計算加速效果顯著，通過能群并行再結(jié)合E指數(shù)優(yōu)化和角通量球諧函數(shù)展開，輸運(yùn)計算模塊加速效果顯著，在相同計算精度下，大幅提高了計算速度。