郭海兵，黃洪文，馬紀(jì)敏，丁文杰

(中國工程物理研究院 核物理與化學(xué)研究所，四川 綿陽 621900)

中子輸運方程是反應(yīng)堆物理分析的基礎(chǔ)，它通過描述大量微觀粒子運動所遵循的微分-積分關(guān)系，確定粒子在時間、能量、幾何空間和速度相空間等7個維度上的分布[1]。求解中子輸運方程獲得中子角通量密度，便能計算出中子與原子核的各類反應(yīng)率，進而得到系統(tǒng)的物理性能。通過數(shù)值離散偏微分方程并解得離散點上近似值的確定論方法，具有計算速度快、可獲取物理量的精細(xì)場分布、可高效多物理耦合等優(yōu)點。確定論方法采用分群方式離散能量變量，對方向變量的處理則衍生了多種方法，目前研究和應(yīng)用最廣泛的是離散縱標(biāo)(SN)法。

對空間變量的離散可分為結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格兩種，結(jié)構(gòu)化網(wǎng)格在對復(fù)雜結(jié)構(gòu)的高保真建模方面存在明顯缺陷?；诜墙Y(jié)構(gòu)網(wǎng)格的有限元方法在中子輸運求解中應(yīng)用后，復(fù)雜結(jié)構(gòu)的高保真建模和計算得以解決[2-4]，但計算量相對較大。因而在全堆芯計算和大空間輻射屏蔽問題中，必須借助現(xiàn)代大規(guī)模并行計算技術(shù)才能實現(xiàn)大型復(fù)雜結(jié)構(gòu)問題的高保真模擬分析[5-6]。

采用共享內(nèi)存的多線程并行方式(OpenMP)容易實現(xiàn)SN方法下各獨立離散方向的并行求解，但受制于單臺計算機的CPU數(shù)量和內(nèi)存資源，這種方式的并行規(guī)模及允許網(wǎng)格量有限。對幾何空間進行區(qū)域分解，并采用獨享內(nèi)存的多進程并行方式(MPI)計算，則能大幅擴展并行規(guī)模[7-8]。

并行自適應(yīng)非結(jié)構(gòu)網(wǎng)格應(yīng)用框架(JAUMIN)是北京應(yīng)用物理與計算數(shù)學(xué)研究所針對科學(xué)計算中的非結(jié)構(gòu)網(wǎng)格應(yīng)用而開發(fā)的數(shù)值模擬支撐平臺，通過封裝高性能數(shù)據(jù)結(jié)構(gòu)、集成成熟數(shù)值算法、屏蔽大規(guī)模并行和網(wǎng)格自適應(yīng)技術(shù)，加速可使用現(xiàn)代高性能計算機的應(yīng)用程序研制[9-11]?；谧钚《擞邢拊蚐N方法[12]的中子輸運程序ENTER是自主開發(fā)的用于支持新型反應(yīng)堆和次臨界包層設(shè)計的數(shù)值模擬軟件，銜接成熟的CAD建模-網(wǎng)格劃分前處理軟件和數(shù)據(jù)可視化后處理軟件，輕松實現(xiàn)復(fù)雜結(jié)構(gòu)的高保真建模、高精度模擬和靈活的數(shù)據(jù)分析。

采用JAUMIN框架進行區(qū)域分解和多進程并行，通過改造并行流水線掃描算法，實現(xiàn)了ENTER程序的高效并行，通過系列基準(zhǔn)問題和若干大型問題驗證了程序的正確性和并行效果。本文將具體介紹基于有限元和SN方法求解中子輸運方程的算法以及區(qū)域分解和并行掃描的實現(xiàn)機制，最后通過部分計算案例展示ENTER程序的計算精度和并行計算能力。

1 中子輸運方程及其有限元SN求解

定常條件下，多群形式的中子輸運方程[1]可寫為：

g=1,2,…,G

(1)

式中：φg=φ(r,g,Ω)，為能群g在位置r處沿方向角Ω的中子角通量密度；Σt,g=Σt(r,g)，為能群g的宏觀總截面；Sg=S(r,g,Ω)，為能群g的各源項之和，包括散射源項Ss,g、裂變源項Sf,g和獨立外源項Se,g。

Sf,g和Se,g均易于計算和處理，本文不作贅述，重點討論Ss,g在SN方法下的表達(dá)形式。散射截面是散射前后運動方向間夾角的函數(shù)，將散射截面和中子角通量密度分別以有限階的勒讓德函數(shù)和球面諧函數(shù)近似展開后，得到散射源項表達(dá)為：

(2)

(3)

(4)

(5)

從而，選定1組離散方向并確定相應(yīng)的權(quán)重值{Ωm,wm}后，即可將式(3)和(4)的積分轉(zhuǎn)化為求和，進而將式(2)轉(zhuǎn)化為以離散方向處的中子角通量密度φg,m=φ(r,g,Ωm)為待求解量的方程。在直角坐標(biāo)系下，離散方向Ωm=(μm,x,μm,y,μm,z)處的多群中子輸運方程為：

g=1,2,…,G

m=1,2,…,M

(6)

在進行有限元離散和求解時，由于式(6)是非自伴隨的，直接應(yīng)用標(biāo)準(zhǔn)Galerkin加權(quán)方法可能會引起解的振蕩，應(yīng)采用迎風(fēng)Petrov-Galerkin(SUPG)[13]等加強穩(wěn)定性的方法，或采用間斷有限元(Discontinuous Galerkin, DG)[4]、最小二乘有限元(Galerkin Least Square, GLS)[12]等方法求解，在高維問題中GLS方法實現(xiàn)較簡便且計算量較小，以下即對GLS方法簡要說明。

有限元方法是在離散單元內(nèi)將待求變量φg,m(r)近似表達(dá)為1組特別選定的基函數(shù){N(r)}的線性疊加，單元內(nèi)各結(jié)點的變量離散值向量Φg,m即作為基函數(shù)的系數(shù)：

(7)

式中，N為單元的基函數(shù)組成的行向量。

L(φg,m)=Sg,m

(8)

(9)

代入算子L的表達(dá)式，并展開積分式，化簡得到式(6)的矩陣形式離散求解方程：

(μm,kμm,iHk,i+Mt,g+μm,kGg,k)Φg,m=

(10)

式中，以指標(biāo)符號形式簡寫了求和。對于不含獨立外源的問題，式(10)對一維和二/三維幾何是通用的，但一維下所求解變量的物理意義略有不同，以表達(dá)式說明：

(11)

對于獨立外源項，由于其不依賴中子角通量密度φg,m，因而基于式(5)的轉(zhuǎn)換后，計算式在一維下比二/三維下多了因子2π。除此之外，不同幾何維度的計算方法和流程完全一致。

式(10)采用源迭代法求解，在每步迭代中還可采用擴散綜合加速(DSA)算法[14-15]對中子通量密度(標(biāo)通量)和中子流密度(中子角通量密度的1階球諧分量)進行修正，以使迭代更快收斂?；蚧谥凶訑?shù)守恒的思想，采用角度再平衡方法進行加速。

2 區(qū)域分解及并行掃描算法

所謂區(qū)域分解，即是把計算問題的幾何結(jié)構(gòu)劃分成多個相對小的子區(qū)域，將各子區(qū)域分配給不同的CPU核進行求解，一方面減小了求解方程組的規(guī)模，另一方面通過并行增加了單位時間運算量，從而提高計算速度。區(qū)域分解算法的關(guān)鍵是通過子區(qū)域之間的數(shù)據(jù)傳遞，使得分區(qū)并行求解(可能需子區(qū)域間迭代)的結(jié)果與全域整體求解的結(jié)果一致。

區(qū)域分解后會產(chǎn)生新的界面(子區(qū)域邊界)，這些界面原本處于計算域內(nèi)部而不需關(guān)注，但成為子區(qū)域邊界后就需指定合適的邊界條件，以使可在子區(qū)域上求解方程。具體來說，采用SN方法求解中子輸運方程時，需為每個方向的入射邊界指定中子角通量密度。JAUMIN框架基于重疊型區(qū)域分解算法來完成子區(qū)域間數(shù)據(jù)傳遞，即將每個子區(qū)域的劃分邊界向相鄰子區(qū)域中延伸1層網(wǎng)格，形成被稱為影像區(qū)的重疊區(qū)域，用于獲取相鄰子區(qū)域相應(yīng)結(jié)點的數(shù)據(jù)[16]。重疊型區(qū)域分解算法的理論基礎(chǔ)是Schwarz交替法及Lions等對它的并行化拓展[17-18]。子區(qū)域及其影像區(qū)合稱為網(wǎng)格片，輸運方程的求解是在網(wǎng)格片上進行的。以一維幾何下兩個相鄰子區(qū)域的情況作示意說明，如圖1所示。

控制各子區(qū)域的網(wǎng)格單元數(shù)量基本相等的情況下，采用JAUMIN框架對方形區(qū)域的非結(jié)構(gòu)網(wǎng)格進行區(qū)域分解形成的16個子區(qū)域如圖2所示。

應(yīng)用程序在JAUMIN框架上實現(xiàn)并行的原理是框架將所有網(wǎng)格片的初始化、并行計算、數(shù)據(jù)傳遞、歸約同步等操作封裝為積分構(gòu)件，在網(wǎng)格層上通過調(diào)用這些積分構(gòu)件形成完整的并行計算流程，依托C++語言的類繼承和虛函數(shù)這兩個重要機制，用戶只需在繼承類的虛函數(shù)中具體實現(xiàn)特定物理問題在單個網(wǎng)格片上的初邊值條件和數(shù)值求解、在網(wǎng)格層上的歸約同步和迭代操作，則原有并行計算流程中的相應(yīng)模塊被繼承類的用戶函數(shù)替代，框架的并行機制和計算流程保持不變而實現(xiàn)不同問題求解[19]。ENTER程序基于JAUMIN框架的區(qū)域分解并行架構(gòu)如圖3所示。

圖1 網(wǎng)格片影像區(qū)及數(shù)據(jù)傳遞示意圖Fig.1 Sketch of ghost areas and data transfer scheme between mesh patches

圖2 非結(jié)構(gòu)網(wǎng)格區(qū)域分解示意Fig.2 Illustration of domain decomposition under unstructured mesh

由于SN輸運方程(式(6))的求解具有方向性，只有入射邊界(Ωm·n<0，n為邊界外法線向量)上的變量已知時，網(wǎng)格片內(nèi)各離散點才可求解。對于真空邊界，該定解值為0；對于反射邊界，該定解值為對稱出射方向在反射點的中子角通量密度，因此需先求解對稱出射方向方程；對于影像區(qū)邊界，該定解值由該方向上游相鄰網(wǎng)格片傳遞過來，需先求解上游網(wǎng)格片。

因而對幾何空間做區(qū)域分解后，各網(wǎng)格片并非完全獨立可解的，本文借鑒間斷有限元方法中使用的并行流水線掃描算法，將各網(wǎng)格片與所有離散方向進行對應(yīng)組合，對每個網(wǎng)格片獨立判斷未求解的各離散方向是否具備可求解條件，只要有1個方向可求解，則將該網(wǎng)格片加入并行隊列，當(dāng)前計算步只并行求解隊列中的網(wǎng)格片上可求解方向?qū)?yīng)的輸運方程。網(wǎng)格片上的某個方向完成求解后，將更新下游若干個相鄰網(wǎng)格片該方向的可求解條件。每個離散方向最上游的網(wǎng)格片將首先具備求解條件，如果某網(wǎng)格片上沒有離散方向具備求解條件，則將跳過當(dāng)前計算步，等待其邊界上可求解條件被更新。這樣分批次并行求解，直到所有網(wǎng)格片上的所有方向完成求解[8-10]。

圖3 ENTER程序基于JAUMIN框架實現(xiàn)并行的總體架構(gòu)Fig.3 Framework of parallel code ENTER based on JAUMIN infrastructure

在每個并行計算步中，各網(wǎng)格片上求解的離散方向一般不同。以圖2所示的二維16個網(wǎng)格片劃分和4個離散方向(對應(yīng)S2)為例，對并行掃描算法在每個計算步中各網(wǎng)格片求解的方向號說明列于表1，表中“#”表示網(wǎng)格片編號，“×”表示當(dāng)前計算步跳過，“√”表示當(dāng)前網(wǎng)格片上的所有離散方向求解完成。

表1 并行掃描算法在每個計算步中各網(wǎng)格片求解的方向號Table 1 Solving direction number of all patches under each step of parallel sweeping algorithm

每個CPU可處理1個或多個網(wǎng)格片，出于負(fù)載平衡考慮，使用的CPU數(shù)應(yīng)為網(wǎng)格片數(shù)量的約數(shù)。無論是針對結(jié)構(gòu)化網(wǎng)格的KBA掃描算法[20-21]，還是針對非結(jié)構(gòu)網(wǎng)格間斷有限元的并行流水線掃描算法[9]或有向圖掃描算法[22-23]，表達(dá)方式和程序?qū)崿F(xiàn)形式不同但基本思想一致，即對入射邊界上變量值已知的基本求解單位做同步計算，只不過結(jié)構(gòu)化網(wǎng)格下依賴關(guān)系完全明確，有向圖算法根據(jù)網(wǎng)格依賴關(guān)系預(yù)先對掃描計算順序作抽象后進行分解并行。本文的掃描算法也是據(jù)此而來，采用無向區(qū)域分解和對網(wǎng)格片的并行掃描，基本求解單位是網(wǎng)格片與離散方向的組合，并行粒度相對較大，通信開銷較小，但對入射邊界是否滿足求解條件的判斷則更復(fù)雜。

由于區(qū)域分解形成的網(wǎng)格片邊界不規(guī)整，對某些離散方向，網(wǎng)格片之間的依賴關(guān)系復(fù)雜甚至相互依賴，因此會出現(xiàn)某些網(wǎng)格片上的某些方向需多次求解的現(xiàn)象，這樣才能完成可求解條件向下游網(wǎng)格片的傳遞。掃描算法在單步求解中存在的偏差會隨著源迭代求解過程逐漸減小，并最終獲得收斂的結(jié)果。但隨離散方向數(shù)的增多，掃描過程會變得更復(fù)雜，雖然可并行度提高了，但掃描過程中出現(xiàn)網(wǎng)格片互相依賴的情況更頻繁。

3 正確性驗證和并行性能分析

通過對系列基準(zhǔn)問題的計算對比完成了ENTER程序驗證。現(xiàn)通過部分問題的結(jié)果，分析ENTER程序的精度和并行性能。對問題的完整描述及截面參數(shù)參閱相關(guān)文獻[24-28]，源迭代收斂準(zhǔn)則為keff殘差<1×10-5且中子通量密度的均方根偏差<1×10-4。

ISSA基準(zhǔn)題是一維、單群、P0階散射的臨界問題，通過在y軸和z軸方向設(shè)置反射邊界，可將該問題擴展為二維和三維問題。對該問題在不同條件下的計算結(jié)果的對比列于表2，ENTER計算結(jié)果為網(wǎng)格無關(guān)解，即通過多次加密網(wǎng)格獲得一致結(jié)果。與參考值相比，S6條件下keff的最大偏差為4.1×10-4，S16條件下的最大偏差為1.7×10-4。由于一維與二/三維采用的離散求積組不同，在數(shù)值積分精度上不同，因此對于S4及以上條件，二/三維的結(jié)果接近，但與一維略有差異。

圖4展示了兩個外源問題的計算結(jié)果對比，圖4a為一維、2群、P3階散射各向異性問題[25]，圖4b為二維、2群、P1階散射深穿透問題[26]，兩個問題均采用S8條件計算。從圖可見，在中子通量密度跨越7個量級的情況下，計算結(jié)果仍與參考值能保持較好的一致性。

表2 ISSA基準(zhǔn)題keff計算結(jié)果對比Table 2 Comparison of keff results for ISSA benchmark

注：1) 西安交通大學(xué)巨海濤等[24]開發(fā)的最小二乘有限元SN中子輸運程序

2) 采用GAUSS求積組

3) 采用偶階矩條件的全對稱求積組[29]

圖4 外源各向異性和深穿透問題計算結(jié)果Fig.4 Comparison of external source problem and deep penetration problem results under heterogeneous scattering condition

外邊界全部為反射條件的問題對SN方法具有特殊的意義，此時每個離散方向在其他所有象限中均存在1個耦合關(guān)聯(lián)方向，所有離散方向都不能獨立求解，并且中子通量密度沒有唯一解，計算過程中需進行循環(huán)迭代和多次歸一化處理。BWR柵元問題是二維、2群、P0階散射的臨界基準(zhǔn)題，外邊界全部為反射條件[27]。采用ENTER程序S8和S16計算的keff與參考值的相對偏差分別為1.2×10-4和4.9×10-5，獲得的兩群中子通量密度分布如圖5所示。

對Takeda和Ikeda發(fā)表的三維小型LWR堆芯基準(zhǔn)題[28](1/8堆芯，控制棒插入)進行精細(xì)網(wǎng)格劃分，并采用不同離散縱標(biāo)數(shù)求解，通過強擴展并行測試來分析ENTER程序的并行性能(表3)。該問題為2群、P0階散射，采用四面體網(wǎng)格劃分，1/8堆芯單元總數(shù)為1.37×105。定義自由度數(shù)=能群數(shù)×方向數(shù)×變量點數(shù)，則S4條件下最大自由度數(shù)為9.91×106，S8條件下最大自由度數(shù)為3.30×107。

圖5 BWR柵元問題的兩群中子通量密度分布Fig.5 Contours of fast and thermal neutron flux density distribution for BWR cell problem

表3 小型LWR堆芯基準(zhǔn)題計算結(jié)果及并行效率分析(keff參考值為0.962 4±0.000 5)Table 3 Simulation result and parallel efficiency analysis for small LWR benchmark problem(reference value of keff is 0.962 4±0.000 5)

注：1) CPU核數(shù)記為NP，表示使用N個節(jié)點，每個節(jié)點使用P個進程，每個進程處理1個網(wǎng)格片

2) 以最少CPU核數(shù)(CL)時的計算時間(TL)為參考，設(shè)定并行效率為100%，隨著并行CPU核數(shù)(CN)的增加，實際計算時間(TN)對應(yīng)的并行效率為CLTL/CNTN

從表3的結(jié)果可見，該問題在S4條件下保持了很好的并行擴展效率，而在S8條件下并行核數(shù)高于160后效率下降很快。這主要是兩方面原因造成的：一是離散方向數(shù)較多的情況下，當(dāng)存在反射邊界條件時，網(wǎng)格片之間的依賴關(guān)系隨網(wǎng)格片增多而更復(fù)雜，掃描算法更耗時間；二是問題規(guī)模不大時，劃分網(wǎng)格片越多，每個網(wǎng)格片內(nèi)部的單元數(shù)越少、影像區(qū)單元數(shù)越多，計算量增加且通信耗時占比更高。

對包含169個六棱柱組件的FBR三維堆芯基準(zhǔn)題(控制棒半插入)[28]進行全堆芯輸運計算，采用4群、P0階散射截面，離散縱標(biāo)數(shù)S8?；谌庵W(wǎng)格，并采用1次全局網(wǎng)格自動加密(將所有的網(wǎng)格邊二分，三維下1個單元被切分為8個小單元)和不同的CPU核數(shù)，各區(qū)域平均的中子通量密度列于表4，獲得的各群中子通量密度剖面分布如圖6所示，計算得到的keff結(jié)果和耗時列于表5。

從表4、5可見，基于ENTER并行程序?qū)υ摯笠?guī)模問題的模擬結(jié)果是準(zhǔn)確的，區(qū)域平均的中子通量密度與SN文獻值相比，ENTER的計算精度顯著更高，與MCNP參考值符合很好，最大偏差出現(xiàn)在控制棒區(qū)域第4群，這可能是該能群的局部吸收率較大而ENTER求解時局部網(wǎng)格密度不夠造成的。對問題的計算時間是可接受的，計算結(jié)果的精度對網(wǎng)格密度較敏感，這主要是由于控制棒周圍的中子角通量密度梯度較大，目前采用的是一階線性/雙線性單元，高網(wǎng)格密度才能準(zhǔn)確描述通量的空間分布。后續(xù)需進一步考慮引入高階單元和局部自適應(yīng)加密技術(shù)，以實現(xiàn)較少網(wǎng)格量下獲得高精度結(jié)果。

表4 區(qū)域平均中子通量密度與參考值[28]的比較Table 4 Comparison of zone-averaged neutron flux density between simulation results and reference values

注：表中( )內(nèi)的值，對于MCNP參考解為統(tǒng)計誤差，其他則為與MCNP參考解的相對偏差

圖6 六棱柱FBR堆芯的中子通量密度二維剖面分布Fig.6 Contours of neutron flux density distribution for hexagonal FBR core

表5 六棱柱FBR堆芯計算結(jié)果及時間(keff參考值為0.983 3±0.000 4)Table 5 Simulation result and time cost of hexagonal FBR core(reference value of keff is 0.983 3±0.000 4)

注：1) 網(wǎng)格片數(shù)增多會使得影像區(qū)的總結(jié)點數(shù)增加，因而在相同的網(wǎng)格下，劃分網(wǎng)格片越多，自由度數(shù)越大

2) 天河Ⅱ號單個節(jié)點的CPU核數(shù)為24

4 討論及結(jié)論

區(qū)域分解是對輸運方程的空間變量的并行處理，并不排斥對離散角度變量的并行求解，并行掃描算法融合了空間和角度兩個變量的并行度，因而當(dāng)網(wǎng)格片上有多個方向同時具備求解條件時，這些方向之間也是可并行求解的。目前的ENTER程序是基于MPI的空間并行，每個網(wǎng)格片每步只能求解1個方向，因而網(wǎng)格片上的角度并行需采用第2級的OpenMP方式實現(xiàn)。

非結(jié)構(gòu)網(wǎng)格最小二乘有限元方法在區(qū)域分解并行機制下可解決大型復(fù)雜結(jié)構(gòu)堆芯的輸運計算問題，獲得準(zhǔn)確的堆芯參數(shù)分布，在2.81×109自由度規(guī)模下計算時間約7.4 h是可接受的，具有一定的工程應(yīng)用價值。但該方法尚不能獲得大規(guī)模的并行擴展，特別是在方向數(shù)較多(S8以上)的情況下，隨著劃分網(wǎng)格片數(shù)量的增多，網(wǎng)格片之間的依賴關(guān)系變得復(fù)雜，并行效率急劇下降。同樣由于該原因，在并行條件下的計算結(jié)果與串行結(jié)果沒有取得完全一致。

間斷有限元方法是克服以上困難的有效途徑，因為間斷有限元方法以單元與方向的組合為基本求解單位，且單元邊界平直，不存在兩個單元互相依賴的問題，因而有望在這種小的粒度下實現(xiàn)大規(guī)模的并行，借助有向圖掃描算法可能獲得較高的并行效率。網(wǎng)格局部自適應(yīng)加密和高階形函數(shù)的應(yīng)用，也是后續(xù)值得研究的內(nèi)容，可減少初始網(wǎng)格量并快速獲得網(wǎng)格無關(guān)解。

在學(xué)習(xí)和使用JAUMIN框架的過程中，得到了中國工程物理研究院高性能數(shù)值模擬軟件中心的楊章等的指導(dǎo)和幫助，在此對其編程框架團隊表示感謝。