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

(1. 中國核動(dòng)力研究設(shè)計(jì)院設(shè)計(jì)所,四川 成都 610213；2. 中國核動(dòng)力研究設(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)程序使用分治策略，例如提供堆芯級(jí)別(Core Level)和柵元級(jí)別(Lattice level)這種主要分類計(jì)算形式[1]。Hydra-SN3D軟件是由中國核動(dòng)力研究設(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)行離散化，公式表明沿著一個(gè)特定方向Ωm輻射通量發(fā)生的變化，求解過程便是定量獲得此變化趨勢。該近似方法是把變量Ω直接離散的數(shù)值方法，即只對(duì)選定的若干個(gè)離散方向Ωm對(duì)中子輸運(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表示方向向量在某個(gè)縱坐標(biāo)方向上(例如直角坐標(biāo)系的x方向坐標(biāo)軸)的離散點(diǎn)數(shù)目。在反應(yīng)堆物理和屏蔽計(jì)算中，對(duì)中子角通量密度分布各向異性比較嚴(yán)重或非均勻性比較強(qiáng)烈的問題(例如柵元或燃料組件)的計(jì)算中都廣泛應(yīng)用[4]。文獻(xiàn)歸納SN數(shù)值方法求解中子輸運(yùn)方程的主要問題為：

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

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

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

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

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

(3)

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

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

(4)

實(shí)際上對(duì)于不同的坐標(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ù)又含有對(duì)方向向量的導(dǎo)數(shù)，所以空間變量兼具曲幾何和直角坐標(biāo)系的特點(diǎn)，在離散過程中將產(chǎn)生對(duì)應(yīng)的外推模型[5]。

(1)菱形差分

針對(duì)每一個(gè)網(wǎng)格，以二維情況為例，(i,j)網(wǎng)格中守恒方程含有φi,j,m,φi±1/2,j,m,φi,j±1/2,m等五個(gè)未知函數(shù)，曲幾何還需要加上φi,j,m±1/2兩個(gè)未知函數(shù)，共7個(gè)未知數(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)階躍近似

對(duì)于每個(gè)維度或方向上，對(duì)帶權(quán)菱形差分使用特定的系數(shù)，此時(shí)對(duì)應(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)行組織，可以對(duì)基準(zhǔn)例題提供的輸運(yùn)截面參數(shù)進(jìn)行手動(dòng)填寫或者使用TRANSX接口程序自動(dòng)生成。Hydra-SN3D軟件具備上散射輸運(yùn)計(jì)算和輸入高階散射截面的功能[6]。

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

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

2.1 前處理

Hydra-SN3D軟件的輸入文件按照FIDO格式進(jìn)行提供，用戶填入內(nèi)容相對(duì)比較自由，但缺乏容錯(cuò)性的支持，當(dāng)用戶提供明顯錯(cuò)誤的輸入時(shí)軟件將正常執(zhí)行而不會(huì)報(bào)錯(cuò)退出。對(duì)Hydra-SN3D軟件的改進(jìn)策略之一便是重新設(shè)計(jì)輸入文件形式，使用Python腳本編程語言書寫前處理模塊，判斷輸入數(shù)據(jù)的不合理內(nèi)容，提升軟件本身的容錯(cuò)性(見圖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)一格式，利用層級(jí)關(guān)系來表達(dá)文本信息。輸入文件格式強(qiáng)調(diào)以數(shù)據(jù)塊名稱進(jìn)行區(qū)分，在每個(gè)數(shù)據(jù)塊內(nèi)可以包含各種基本數(shù)據(jù)類型及其值，例如單變量或數(shù)組變量等。為了提升易用性，每個(gè)數(shù)據(jù)塊的數(shù)據(jù)格式至多支持三級(jí)結(jié)構(gòu)，每級(jí)結(jié)構(gòu)的數(shù)據(jù)塊名稱不能重合，不同級(jí)的數(shù)據(jù)塊名稱可以重合，訪問時(shí)按照A∶B∶C這種方式來顯式說明，表1的示例中數(shù)據(jù)塊Test1中變量var1有兩個(gè)，訪問時(shí)分別為Test1∶var1或者Test1∶var5∶var1。

表1 NAC輸入文件格式示例

各個(gè)數(shù)據(jù)塊在使用Python編程實(shí)現(xiàn)時(shí)可以按行順序讀取，也可以按照變量名稱選擇性讀取，正則表達(dá)式等機(jī)制則應(yīng)用于關(guān)鍵字或數(shù)值的匹配。為了與數(shù)值計(jì)算部分有效連接，在Hydra-SN3D軟件的改進(jìn)過程中使用Python外處理過程讀取新的NAC格式的SN計(jì)算輸入文件，共計(jì)約40個(gè)變量。這些變量的值分為單數(shù)據(jù)或者數(shù)組數(shù)據(jù)兩種情況，讀取完成后然后進(jìn)行容錯(cuò)性判斷，隔離基本的輸入錯(cuò)誤或者一些不合理的輸入(例如截面數(shù)值小于0，堆芯材料布置異常)，然后自動(dòng)生成舊的FIDO格式輸入文件，同時(shí)鏈接Hydra-SN3D軟件執(zhí)行運(yùn)算。這種方案提升了計(jì)算用戶建立輸入文件的效率，且部分錯(cuò)誤在執(zhí)行Hydra-SN3D軟件之前便由新的前處理進(jìn)行了隔離。與此同時(shí)，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 多級(jí)并行

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

(9)

φ=Dψ

(10)

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

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

S——散射矩陣。

待求解對(duì)象ψ代表角通量;

φ為通量矩。XYZ三維結(jié)構(gòu)化網(wǎng)格情況下KBA掃描計(jì)算的實(shí)質(zhì)是求解L-1，它將三維網(wǎng)格進(jìn)行二維劃分，總的網(wǎng)格在每個(gè)方向上為(I,J,K)，每個(gè)區(qū)域(Ia,Jb,K)映射至每個(gè)MPI進(jìn)程進(jìn)行并行求解。整個(gè)掃描計(jì)算過程對(duì)每個(gè)卦限進(jìn)行各個(gè)離散方向的網(wǎng)格掃描求解通量矩和角通量信息。假定整個(gè)結(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)，例如每個(gè)卦限內(nèi)離散角的數(shù)量變化由S4到S32，相應(yīng)的并行效率變化范圍可以達(dá)到0.75到0.99[3]。這里對(duì)Hydra-SN3D軟件的KBA并行掃描算法增加OpenMP并行算法，借鑒文獻(xiàn)[8]給出的MPI-OpenMP-SIMD多級(jí)并行在SN方法中的應(yīng)用，本算法編程實(shí)現(xiàn)時(shí)將依序?qū)Χ鄠€(gè)卦限內(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)格序號(hào)(i,j,k)相關(guān)。通量矩變量φup、φpsi等存在著先后依賴關(guān)系，因此徑向(i,j)的掃描具備流水線結(jié)構(gòu)。引入計(jì)算是否完成的標(biāo)志位二維數(shù)組，其維度為(j,i)，每個(gè)維度的大小與掃描計(jì)算所需的上一個(gè)維度數(shù)據(jù)φup,j相同。每個(gè)線程在實(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í)并行方案首先進(jìn)行C/C++編程語言的等價(jià)轉(zhuǎn)換，調(diào)整數(shù)據(jù)結(jié)構(gòu)易于多線程并行；然后將外推模型的計(jì)算進(jìn)行局部擴(kuò)展，輔助展開流水線結(jié)構(gòu)的計(jì)算語句。整個(gè)多級(jí)并行的編程實(shí)現(xiàn)命名為NAC4R，它將在數(shù)值驗(yàn)算部分給出性能的對(duì)比。

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

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

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

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

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

這里也對(duì)每個(gè)例題的通量分布進(jìn)行可視化，數(shù)值驗(yàn)算過程逐一對(duì)比了串行結(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é)果來看，單個(gè)MPI進(jìn)程運(yùn)行時(shí)Hydra-SN3D軟件的效率較低，結(jié)構(gòu)化網(wǎng)格導(dǎo)致的對(duì)稱性使得多進(jìn)程并行可以通過不同的(Ia,Ib)組合來取得更好的多進(jìn)程并行增益。KUCA例題使用100個(gè)進(jìn)程的加速比為62.4，LWR-pool例題使用64個(gè)進(jìn)程的加速比為45.6，VENUS2例題使用100個(gè)進(jìn)程的加速比為80.1。

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

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

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

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

4 結(jié)論

在Hydra-SN3D的應(yīng)用過程中持續(xù)對(duì)該軟件進(jìn)行發(fā)展和改進(jìn)，一方面對(duì)已有的前處理部分進(jìn)行了改進(jìn)，重新設(shè)計(jì)輸入文件及其前處理流程，提升計(jì)算的容錯(cuò)性，自動(dòng)檢查計(jì)算用戶輸入?yún)?shù)的合理性；另一方面對(duì)Hydra-SN3D的編程語言實(shí)施等價(jià)變換和結(jié)構(gòu)改進(jìn)，提升單進(jìn)程的運(yùn)行效率。提出并實(shí)施多級(jí)并行算法，驗(yàn)算過程論證了整體計(jì)算效率的提升情況(均取得一個(gè)數(shù)量級(jí)的提升)。隨著后續(xù)堆芯數(shù)值計(jì)算研究工作的開展，將借鑒解決此類特征問題的思路，引入形式化方法，對(duì)更多堆芯數(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)。