反應(yīng)堆中子學(xué)屏蔽計(jì)算軟件Hydra-SN3D的改進(jìn)

明平洲,唐松乾,安 萍,蘆 韡,劉 東,余紅星

(1. 中國核動力研究設(shè)計(jì)院設(shè)計(jì)所,四川 成都 610213；2. 中國核動力研究設(shè)計(jì)院核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室,四川 成都 610213)

中子、光子、電子等一些粒子可以導(dǎo)出波爾茲曼微分—積分守恒方程，統(tǒng)一被稱為輸運(yùn)方程。核反應(yīng)堆工程領(lǐng)域的輻射屏蔽輸運(yùn)計(jì)算主要表現(xiàn)為中子、光子等的物理參數(shù)求解。輸運(yùn)分析計(jì)算程序的核心，涵蓋了幾何、物理方法、核數(shù)據(jù)等多方面。現(xiàn)有的反應(yīng)堆中子學(xué)輸運(yùn)程序使用分治策略，例如提供堆芯級別(Core Level)和柵元級別(Lattice level)這種主要分類計(jì)算形式[1]。Hydra-SN3D軟件是由中國核動力研究設(shè)計(jì)院和西安交大研制的SN方法輸運(yùn)計(jì)算軟件，使用源迭代法求解離散后的穩(wěn)態(tài)粒子輸運(yùn)方程，在多種差分格式基礎(chǔ)上使用KBA并行算法實(shí)現(xiàn)空間網(wǎng)格的并行掃描計(jì)算。該軟件內(nèi)部計(jì)算核心為類似于DORT和TORT等知名軟件的核心輸運(yùn)求解器，外圍則逐步集成較多的輻射屏蔽計(jì)算功能，例如中子/光子耦合計(jì)算功能[2-3]。

1 數(shù)值方法

1.1 SN方法

中子或光子的輻射效應(yīng)通常由波爾茲曼輸運(yùn)方程的離散縱標(biāo)(SN, Discrete Ordinates)形式進(jìn)行建模，它的主要特點(diǎn)在于空間角度的離散。

Ωm·Img+(σA+σS)Img

(1)

輻射通量Img(x,y,z)按照能群、角方向進(jìn)行離散化，公式表明沿著一個特定方向Ωm輻射通量發(fā)生的變化，求解過程便是定量獲得此變化趨勢。該近似方法是把變量Ω直接離散的數(shù)值方法，即只對選定的若干個離散方向Ωm對中子輸運(yùn)方程求解。從中子輸運(yùn)方程求出φ(r,E,Ωm)后，關(guān)于方向Ω的有關(guān)積分則用數(shù)值積分來近似表示。

(2)

式中,求積系數(shù)ωm、離散方向及其數(shù)目取決于計(jì)算精度的要求，這便構(gòu)成離散縱標(biāo)法(SN方法，又寫作SN)。此處下標(biāo)N表示方向向量在某個縱坐標(biāo)方向上(例如直角坐標(biāo)系的x方向坐標(biāo)軸)的離散點(diǎn)數(shù)目。在反應(yīng)堆物理和屏蔽計(jì)算中，對中子角通量密度分布各向異性比較嚴(yán)重或非均勻性比較強(qiáng)烈的問題(例如柵元或燃料組件)的計(jì)算中都廣泛應(yīng)用[4]。文獻(xiàn)歸納SN數(shù)值方法求解中子輸運(yùn)方程的主要問題為：

1)角度方向的離散、數(shù)目及求積組的選取；

2)針對空間、方向變量的中子輸運(yùn)方程的離散化方法及離散方程組的獲得；

3)離散方程的合理求解，包括加速收斂等。

以簡化形式書寫多群近似中子輸運(yùn)方程：

Ω·φg(r,Ω)+∑t,gφg(r,Ω)=Qg(r,Ω)

(3)

在選定求積組(Hydra-SN3D軟件默認(rèn)使用全對稱求積組，也可以輸入其他類型的求解組)之后，對每一個選定的離散方向附近區(qū)域進(jìn)行積分，可以得到對應(yīng)的離散方程：

[Ω·φ(r,Ω)]m+∑tφm(r)=Qm(r)

(4)

實(shí)際上對于不同的坐標(biāo)系統(tǒng)，式(4)的求解具有不同的離散表達(dá)式。常使用的是直角坐標(biāo)系(x,y,z)和圓柱坐標(biāo)系(r,φ,z)進(jìn)行確定論的SN方法推導(dǎo)。

1.2 外推模型

空間變量在中子輸運(yùn)方程中含有空間變量的導(dǎo)數(shù)又含有對方向向量的導(dǎo)數(shù)，所以空間變量兼具曲幾何和直角坐標(biāo)系的特點(diǎn)，在離散過程中將產(chǎn)生對應(yīng)的外推模型[5]。

(1)菱形差分

針對每一個網(wǎng)格，以二維情況為例，(i,j)網(wǎng)格中守恒方程含有φi,j,m,φi±1/2,j,m,φi,j±1/2,m等五個未知函數(shù)，曲幾何還需要加上φi,j,m±1/2兩個未知函數(shù)，共7個未知數(shù)。菱形差分便是補(bǔ)充式(5)所示的關(guān)系式，便于使用待定系數(shù)法來求解這些未知函數(shù)。

(5)

其他維度或方向上有類似的菱形差分格式，公式表明網(wǎng)格入射中子角通量密度與出射中子角通量密度具備線性關(guān)系的連續(xù)性。

(2)帶權(quán)菱形差分

菱形差分近似認(rèn)為網(wǎng)格中心的中子通量密度和網(wǎng)格表面中子通量密度按照直線變化，帶權(quán)菱形差分則使用更為一般的線性關(guān)系來表示它們之間的關(guān)系。

φi,j,m=0.5(1+a)·φi+1/2,j,m+ 0.5(1-a)·φi-1/2,j,m

(6)

(3)階躍近似

對于每個維度或方向上，對帶權(quán)菱形差分使用特定的系數(shù)，此時對應(yīng)于階躍近似。

(7)

式中:μm——方向余弦。

(4)θ加權(quán)菱形差分

該模型得到的結(jié)果更接近線性結(jié)果，使用可調(diào)節(jié)幾何外推參數(shù)γi(θ)來代替一般線性關(guān)系中的系數(shù)。

φi,j,m=γi(θ)·φi+1/2,j,m+(1-γi(θ))·φi-1/2,j,m,0.5≤γi(θ)≤γi(θ)≤1.0

(8)

θ(0≤θ≤1)作為用戶輸入?yún)?shù)，用于調(diào)整網(wǎng)格內(nèi)通量分布的形狀。θ在DORT/TORT程序內(nèi)的建議值為0.9。Hydra-SN3D軟件同樣將權(quán)值默認(rèn)取為0.9，也可以由用戶輸入文件進(jìn)行讀取。

1.3 核數(shù)據(jù)

在Hydra-SN3D軟件內(nèi)核數(shù)據(jù)指代的是問題獨(dú)立的多群輸運(yùn)常數(shù)截面數(shù)據(jù)，Hydra-SN3D軟件存放截面數(shù)據(jù)的文件遵循ANISN格式，它由TRANSX接口程序?qū)镜亩嗳撼?shù)庫進(jìn)一步處理得到計(jì)算所需的輸運(yùn)截面。實(shí)際計(jì)算過程中這些核數(shù)據(jù)參數(shù)按照FORTRAN編程語言提供的UNFORMATTED格式進(jìn)行組織，可以對基準(zhǔn)例題提供的輸運(yùn)截面參數(shù)進(jìn)行手動填寫或者使用TRANSX接口程序自動生成。Hydra-SN3D軟件具備上散射輸運(yùn)計(jì)算和輸入高階散射截面的功能[6]。

2 Hydra-SN3D的改進(jìn)

根據(jù)Hydra-SN3D軟件的實(shí)際使用情況，進(jìn)行以下兩個方面的改進(jìn)：一方面是提升輸入文件在讀取過程中的容錯性；另一方面是進(jìn)一步提升計(jì)算效率。

2.1 前處理

Hydra-SN3D軟件的輸入文件按照FIDO格式進(jìn)行提供，用戶填入內(nèi)容相對比較自由，但缺乏容錯性的支持，當(dāng)用戶提供明顯錯誤的輸入時軟件將正常執(zhí)行而不會報(bào)錯退出。對Hydra-SN3D軟件的改進(jìn)策略之一便是重新設(shè)計(jì)輸入文件形式，使用Python腳本編程語言書寫前處理模塊，判斷輸入數(shù)據(jù)的不合理內(nèi)容，提升軟件本身的容錯性(見圖1)。

圖1 外部前處理的改進(jìn)策略Fig.1 Improvement strategy of external pre-processing

新的輸入文件格式命名為NAC(Numerical Arithmetic Context)，采用類似于互聯(lián)網(wǎng)應(yīng)用常使用的JSON方式的統(tǒng)一格式，利用層級關(guān)系來表達(dá)文本信息。輸入文件格式強(qiáng)調(diào)以數(shù)據(jù)塊名稱進(jìn)行區(qū)分，在每個數(shù)據(jù)塊內(nèi)可以包含各種基本數(shù)據(jù)類型及其值，例如單變量或數(shù)組變量等。為了提升易用性，每個數(shù)據(jù)塊的數(shù)據(jù)格式至多支持三級結(jié)構(gòu)，每級結(jié)構(gòu)的數(shù)據(jù)塊名稱不能重合，不同級的數(shù)據(jù)塊名稱可以重合，訪問時按照A∶B∶C這種方式來顯式說明，表1的示例中數(shù)據(jù)塊Test1中變量var1有兩個，訪問時分別為Test1∶var1或者Test1∶var5∶var1。

表1 NAC輸入文件格式示例

各個數(shù)據(jù)塊在使用Python編程實(shí)現(xiàn)時可以按行順序讀取，也可以按照變量名稱選擇性讀取，正則表達(dá)式等機(jī)制則應(yīng)用于關(guān)鍵字或數(shù)值的匹配。為了與數(shù)值計(jì)算部分有效連接，在Hydra-SN3D軟件的改進(jìn)過程中使用Python外處理過程讀取新的NAC格式的SN計(jì)算輸入文件，共計(jì)約40個變量。這些變量的值分為單數(shù)據(jù)或者數(shù)組數(shù)據(jù)兩種情況，讀取完成后然后進(jìn)行容錯性判斷，隔離基本的輸入錯誤或者一些不合理的輸入(例如截面數(shù)值小于0，堆芯材料布置異常)，然后自動生成舊的FIDO格式輸入文件，同時鏈接Hydra-SN3D軟件執(zhí)行運(yùn)算。這種方案提升了計(jì)算用戶建立輸入文件的效率，且部分錯誤在執(zhí)行Hydra-SN3D軟件之前便由新的前處理進(jìn)行了隔離。與此同時，NAC輸入文件格式還具備耦合前處理的能力，即將數(shù)值計(jì)算部分作為庫的形式進(jìn)行耦合，使用Python-SWIG將Python讀取的輸入文件數(shù)據(jù)傳遞給數(shù)值計(jì)算的編程語言部分(通常為C/C++或FORTRAN)，實(shí)現(xiàn)內(nèi)耦合的技術(shù)目標(biāo)，可以應(yīng)用于堆芯物理-熱工耦合數(shù)值計(jì)算[7]。

2.2 多級并行

SN方法在結(jié)構(gòu)化網(wǎng)格的條件下有著著名的KBA并行掃描算法[2,8]，Hydra-SN3D對波爾茲曼輸運(yùn)方程在方向和角度離散化之后的數(shù)學(xué)方程進(jìn)行合并，構(gòu)造所有通量矩的線性系統(tǒng)來進(jìn)行掃描計(jì)算。

(9)

φ=Dψ

(10)

L——差分輸運(yùn)算符；

M——通量矩離散化算符矩陣；

S——散射矩陣。

待求解對象ψ代表角通量;

φ為通量矩。XYZ三維結(jié)構(gòu)化網(wǎng)格情況下KBA掃描計(jì)算的實(shí)質(zhì)是求解L-1，它將三維網(wǎng)格進(jìn)行二維劃分，總的網(wǎng)格在每個方向上為(I,J,K)，每個區(qū)域(Ia,Jb,K)映射至每個MPI進(jìn)程進(jìn)行并行求解。整個掃描計(jì)算過程對每個卦限進(jìn)行各個離散方向的網(wǎng)格掃描求解通量矩和角通量信息。假定整個結(jié)構(gòu)化網(wǎng)格被區(qū)域分解為5×2，在第一卦限內(nèi)進(jìn)程之間具有如圖2所示的波前網(wǎng)格區(qū)域掃描結(jié)構(gòu)。進(jìn)程9首先進(jìn)行計(jì)算，進(jìn)程4和進(jìn)程8接收到進(jìn)程9發(fā)送的邊界數(shù)據(jù)后可以并行計(jì)算。

圖2 多進(jìn)程通信示意圖Fig.2 Schematic of multi-process communication

按照文獻(xiàn)給出的數(shù)據(jù)，理論的并行效率與較多輸入?yún)?shù)相關(guān)，例如每個卦限內(nèi)離散角的數(shù)量變化由S4到S32，相應(yīng)的并行效率變化范圍可以達(dá)到0.75到0.99[3]。這里對Hydra-SN3D軟件的KBA并行掃描算法增加OpenMP并行算法，借鑒文獻(xiàn)[8]給出的MPI-OpenMP-SIMD多級并行在SN方法中的應(yīng)用，本算法編程實(shí)現(xiàn)時將依序?qū)Χ鄠€卦限內(nèi)XY方向上的掃描計(jì)算循環(huán)增加OpenMP編譯指導(dǎo)語句，即在已有MPI多進(jìn)程運(yùn)行的基礎(chǔ)上添加多線程并行。以第一卦限內(nèi)XY方向上的掃描計(jì)算為例，算法描述見表2。

表2 多線程并行算法的描述

在上述計(jì)算過程中，菱形權(quán)重差分外推模型的計(jì)算子函數(shù)T具備獨(dú)立性，僅僅與網(wǎng)格序號(i,j,k)相關(guān)。通量矩變量φup、φpsi等存在著先后依賴關(guān)系，因此徑向(i,j)的掃描具備流水線結(jié)構(gòu)。引入計(jì)算是否完成的標(biāo)志位二維數(shù)組，其維度為(j,i)，每個維度的大小與掃描計(jì)算所需的上一個維度數(shù)據(jù)φup,j相同。每個線程在實(shí)際執(zhí)行計(jì)算之前判斷標(biāo)志位done，當(dāng)標(biāo)志位表示數(shù)據(jù)已就緒，則馬上執(zhí)行計(jì)算；當(dāng)標(biāo)志位表示數(shù)據(jù)還未提供，則線程阻塞。計(jì)算子函數(shù)T在Hydra-SN3D軟件內(nèi)是FORTRAN編程語言書寫，為了實(shí)現(xiàn)數(shù)據(jù)一致性，多級并行方案首先進(jìn)行C/C++編程語言的等價轉(zhuǎn)換，調(diào)整數(shù)據(jù)結(jié)構(gòu)易于多線程并行；然后將外推模型的計(jì)算進(jìn)行局部擴(kuò)展，輔助展開流水線結(jié)構(gòu)的計(jì)算語句。整個多級并行的編程實(shí)現(xiàn)命名為NAC4R，它將在數(shù)值驗(yàn)算部分給出性能的對比。

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

對Hydra-SN3D軟件的改進(jìn)涉及程序代碼的轉(zhuǎn)換、結(jié)構(gòu)的調(diào)整和多級并行編程，這里使用三個計(jì)算問題來對軟件進(jìn)行完整驗(yàn)算，確保計(jì)算與改進(jìn)之前的結(jié)果一致，且同時驗(yàn)證前處理的正確實(shí)現(xiàn)。整體運(yùn)行環(huán)境使用單個節(jié)點(diǎn)具備20個處理器核心的集群計(jì)算機(jī)。如表3所示。

表3 計(jì)算例題說明

在集群上對Hydra-SN3D(MPI多進(jìn)程并行)運(yùn)行上述三個例題，圖3中顯示了它們的計(jì)算時間開銷，縱軸計(jì)算時間取LOG10對數(shù)來反映加速趨勢。

圖3 源Hydra-SN3D軟件的多進(jìn)程運(yùn)行統(tǒng)計(jì)Fig.3 Statistics of multi-process operation of Hydra-SN3D

這里也對每個例題的通量分布進(jìn)行可視化，數(shù)值驗(yàn)算過程逐一對比了串行結(jié)果和并行結(jié)果，由于串行運(yùn)行和并行運(yùn)行的結(jié)果相同，因此僅繪制出串行運(yùn)行情況的通量分布來表明未省略該驗(yàn)算步驟。如圖4所示。

(a)KUCA例題的通量分布(g=1)

(b)LWR-pool例題的通量分布(g=1)

(c)VENUS2例題的通量分布(g=10)

分析1：原始的Hydra-SN3D軟件使用MPI多進(jìn)程并行，它在各種規(guī)模例題下存在著不同的加速特點(diǎn)。從計(jì)算結(jié)果來看，單個MPI進(jìn)程運(yùn)行時Hydra-SN3D軟件的效率較低，結(jié)構(gòu)化網(wǎng)格導(dǎo)致的對稱性使得多進(jìn)程并行可以通過不同的(Ia,Ib)組合來取得更好的多進(jìn)程并行增益。KUCA例題使用100個進(jìn)程的加速比為62.4，LWR-pool例題使用64個進(jìn)程的加速比為45.6，VENUS2例題使用100個進(jìn)程的加速比為80.1。

多級并行方案NAC4R存在著兩個并行層次，第一個層次是對MPI并行部分重新組織數(shù)據(jù)結(jié)構(gòu)，更快速地實(shí)現(xiàn)區(qū)域分解的MPI多進(jìn)程并行；另一個層次是在徑向掃描計(jì)算過程中引入OpenMP多線程并行，增強(qiáng)局部計(jì)算效率。這里對KUCA例題和LWR-pool例題這兩個結(jié)構(gòu)化網(wǎng)格對象的計(jì)算效率進(jìn)行統(tǒng)計(jì)，觀察OpenMP對計(jì)算效率的進(jìn)一步改善情況。如表4所示。

表4 多級并行的計(jì)算時間統(tǒng)計(jì)

分析2：通過程序改寫和數(shù)據(jù)的重新抽象和組織，單進(jìn)程的計(jì)算效率得到了較大提升，整體性能提升了一個數(shù)量級。對于每個基準(zhǔn)例題，對每個進(jìn)程使用8個線程實(shí)施的多線程并行均起到了加速的作用，KUCA例題和LWR-pool例題相比源Hydra-SN3D軟件的運(yùn)行結(jié)果，加速比分別達(dá)到132.8和205.4。此外統(tǒng)計(jì)兩個例題使用一個進(jìn)程聯(lián)合8個線程的計(jì)算效率，可以得到各自的線程加速效果為1.8和1.6，這表明計(jì)算資源充足的情況下流水線掃描結(jié)構(gòu)的多級并行能夠取得更高的性能增益，且能夠使用更多的計(jì)算資源。

驗(yàn)算過程中也對比了串行結(jié)果和并行結(jié)果，引入OpenMP并行之后的多級并行結(jié)果與原始的Hydra-SN3D軟件運(yùn)行結(jié)果保持一致。

4 結(jié)論

在Hydra-SN3D的應(yīng)用過程中持續(xù)對該軟件進(jìn)行發(fā)展和改進(jìn)，一方面對已有的前處理部分進(jìn)行了改進(jìn)，重新設(shè)計(jì)輸入文件及其前處理流程，提升計(jì)算的容錯性，自動檢查計(jì)算用戶輸入?yún)?shù)的合理性；另一方面對Hydra-SN3D的編程語言實(shí)施等價變換和結(jié)構(gòu)改進(jìn)，提升單進(jìn)程的運(yùn)行效率。提出并實(shí)施多級并行算法，驗(yàn)算過程論證了整體計(jì)算效率的提升情況(均取得一個數(shù)量級的提升)。隨著后續(xù)堆芯數(shù)值計(jì)算研究工作的開展，將借鑒解決此類特征問題的思路，引入形式化方法，對更多堆芯數(shù)值算法進(jìn)行優(yōu)化、改進(jìn)。此外，流水線并行掃描結(jié)構(gòu)存在于較多的堆芯數(shù)值計(jì)算過程中，其加速效果存在著IJ/(I+J+1)的理論推導(dǎo)，如何在實(shí)際軟件編程實(shí)現(xiàn)中逼近這一理論上限和取得更好的增益將是下一步的技術(shù)研究重點(diǎn)。