先進中子學柵格程序KYLIN-Ⅱ輸運模塊并行優化開發

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

(1.中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都 610213；2.北京應用物理與計算數學研究所，北京 100094)

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

1 二維幾何下特征線法求解輸運方程

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

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

(1)

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

網格i的標通量計算公式為：

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

(2)

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

在求解中子輸運問題時，各向異性散射會對計算結果有較大的影響[15]，KYLIN-Ⅱ程序截面庫包含了高階Pn散射截面，在特征線中子輸運模塊中，考慮了各向異性散射，故式(2)中包含了各向異性散射源項，KYLIN-Ⅱ程序對散射源項中的散射截面中的角度應用勒讓德多項式進行展開，展開后的散射源項公式為：

Qscat,g,i,k=

(3)

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

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

2 特征線輸運模塊并行優化方案設計

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

圖1 AFA3G組件幾何網格劃分Fig.1 Mesh division of AFA3G assembly

2.1 高階散射源計算優化策略

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

角通量應用球諧函數展開的公式為：

(4)

(5)

式(3)中的勒讓德函數應用球諧函數表示：

(6)

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

(7)

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

(8)

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

cosj(φk-φk′)

(9)

可見結果一定是實數，沒有虛部。因此構造一個函數Zlj，滿足：

(10)

則有：

(11)

通過構造Zlj函數，在計算過程中不需計算和存儲復數。

2.2 基于MPI的能群并行策略

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

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

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

3 數值驗證

3.1 C5G7 UO2組件基準題

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

m=1,2,…,N

φm∈(0,π/2)

(12)

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

圖3 C5G7基準題的UO2組件Fig.3 Geometry of C5G7 UO2 assembly

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

表1 C5G7 UO2組件計算結果Table 1 Calculation result of C5G7 UO2 assembly

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

3.2 AFA3G燃料組件基準題

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

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

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

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

表2 AFA3G燃料組件串行優化計算結果Table 2 Serial optimization result of AFA3G fuel assembly

表3 AFA3G燃料組件并行優化計算結果Table 3 Parallel optimization result of AFA3G fuel assembly

4 結論

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