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

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

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

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

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

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

(1)

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

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

ψ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)問題時(shí)，各向異性散射會(huì)對(duì)計(jì)算結(jié)果有較大的影響[15]，KYLIN-Ⅱ程序截面庫包含了高階Pn散射截面，在特征線中子輸運(yùn)模塊中，考慮了各向異性散射，故式(2)中包含了各向異性散射源項(xiàng)，KYLIN-Ⅱ程序?qū)ι⑸湓错?xiàng)中的散射截面中的角度應(yīng)用勒讓德多項(xiàng)式進(jìn)行展開，展開后的散射源項(xiàng)公式為：

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ù)用于處理特定的問題，計(jì)算效率高；反向延長追蹤技術(shù)用于處理任意幾何問題，幾何適應(yīng)性強(qiáng)。

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

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

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

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

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

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

(4)

(5)

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

(6)

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

(7)

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

(8)

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

cosj(φk-φk′)

(9)

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

(10)

則有：

(11)

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

2.2 基于MPI的能群并行策略

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

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

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

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

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

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

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

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

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

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

針對(duì)AFA3G燃料組件進(jìn)行正確性測(cè)試和性能測(cè)試，AFA3G燃料組件含有呈17×17方形排列的289個(gè)柵元，包括264根燃料棒、24根鋯合金導(dǎo)向管和1根鋯合金測(cè)量管。燃料棒柵距為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%時(shí)的計(jì)算，導(dǎo)向管和測(cè)量管內(nèi)部全部填充輕水慢化劑。

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

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

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

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

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

4 結(jié)論

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