有限元SN中子輸運(yùn)模擬的區(qū)域分解并行

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

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

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

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

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

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

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

1 中子輸運(yùn)方程及其有限元SN求解

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

g=1,2,…,G

(1)

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

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

(2)

(3)

(4)

(5)

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

g=1,2,…,G

m=1,2,…,M

(6)

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

有限元方法是在離散單元內(nèi)將待求變量φg,m(r)近似表達(dá)為1組特別選定的基函數(shù){N(r)}的線性疊加，單元內(nèi)各結(jié)點(diǎn)的變量離散值向量Φ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)符號(hào)形式簡寫了求和。對(duì)于不含獨(dú)立外源的問題，式(10)對(duì)一維和二/三維幾何是通用的，但一維下所求解變量的物理意義略有不同，以表達(dá)式說明：

(11)

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

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

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

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

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

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

應(yīng)用程序在JAUMIN框架上實(shí)現(xiàn)并行的原理是框架將所有網(wǎng)格片的初始化、并行計(jì)算、數(shù)據(jù)傳遞、歸約同步等操作封裝為積分構(gòu)件，在網(wǎng)格層上通過調(diào)用這些積分構(gòu)件形成完整的并行計(jì)算流程，依托C++語言的類繼承和虛函數(shù)這兩個(gè)重要機(jī)制，用戶只需在繼承類的虛函數(shù)中具體實(shí)現(xiàn)特定物理問題在單個(gè)網(wǎng)格片上的初邊值條件和數(shù)值求解、在網(wǎng)格層上的歸約同步和迭代操作，則原有并行計(jì)算流程中的相應(yīng)模塊被繼承類的用戶函數(shù)替代，框架的并行機(jī)制和計(jì)算流程保持不變而實(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輸運(yùn)方程(式(6))的求解具有方向性，只有入射邊界(Ωm·n<0，n為邊界外法線向量)上的變量已知時(shí)，網(wǎng)格片內(nèi)各離散點(diǎn)才可求解。對(duì)于真空邊界，該定解值為0；對(duì)于反射邊界，該定解值為對(duì)稱出射方向在反射點(diǎn)的中子角通量密度，因此需先求解對(duì)稱出射方向方程；對(duì)于影像區(qū)邊界，該定解值由該方向上游相鄰網(wǎng)格片傳遞過來，需先求解上游網(wǎng)格片。

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

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

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

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

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

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

3 正確性驗(yàn)證和并行性能分析

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

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

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

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

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

2) 采用GAUSS求積組

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

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

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

對(duì)Takeda和Ikeda發(fā)表的三維小型LWR堆芯基準(zhǔn)題[28](1/8堆芯，控制棒插入)進(jìn)行精細(xì)網(wǎng)格劃分，并采用不同離散縱標(biāo)數(shù)求解，通過強(qiáng)擴(kuò)展并行測試來分析ENTER程序的并行性能(表3)。該問題為2群、P0階散射，采用四面體網(wǎng)格劃分，1/8堆芯單元總數(shù)為1.37×105。定義自由度數(shù)=能群數(shù)×方向數(shù)×變量點(diǎn)數(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)題計(jì)算結(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個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)使用P個(gè)進(jìn)程，每個(gè)進(jìn)程處理1個(gè)網(wǎng)格片

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

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

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

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

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

表5 六棱柱FBR堆芯計(jì)算結(jié)果及時(shí)間(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ù)增多會(huì)使得影像區(qū)的總結(jié)點(diǎn)數(shù)增加，因而在相同的網(wǎng)格下，劃分網(wǎng)格片越多，自由度數(shù)越大

2) 天河Ⅱ號(hào)單個(gè)節(jié)點(diǎn)的CPU核數(shù)為24

4 討論及結(jié)論

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

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

間斷有限元方法是克服以上困難的有效途徑，因?yàn)殚g斷有限元方法以單元與方向的組合為基本求解單位，且單元邊界平直，不存在兩個(gè)單元互相依賴的問題，因而有望在這種小的粒度下實(shí)現(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)和幫助，在此對(duì)其編程框架團(tuán)隊(duì)表示感謝。