鄒 航，梁 亮，張 乾,*，宋佩濤，趙 強(qiáng)

(1.哈爾濱工程大學(xué) 核安全與仿真技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001；2.西安核創(chuàng)能源科技有限公司，陜西 西安 710077；3.中國(guó)輻射防護(hù)研究院，山西 太原 030006)

為解決中子在不均勻介質(zhì)中的輸運(yùn)問(wèn)題，特征線方法(MOC)以其幾何普適性的優(yōu)點(diǎn)，成為近年來(lái)反應(yīng)堆物理設(shè)計(jì)分析的研究熱點(diǎn)。為提高三維非均勻全堆芯輸運(yùn)計(jì)算效率，國(guó)際上提出2D/1D耦合方法來(lái)加速計(jì)算[1-3]，其中高效率的2D-MOC求解器在全堆輸運(yùn)計(jì)算中占有重要地位。由于MOC方法本身具有計(jì)算量大、收斂較慢的不足之處，且在中子動(dòng)力學(xué)計(jì)算中變得尤為明顯，因此國(guó)內(nèi)外提出了一系列加速收斂方法在數(shù)值層面上提高M(jìn)OC的動(dòng)力學(xué)計(jì)算效率。其中，粗網(wǎng)有限差分方法(CMFD)能顯著降低每個(gè)時(shí)間步輸運(yùn)計(jì)算的迭代次數(shù)，但由于動(dòng)力學(xué)問(wèn)題本身收斂較慢，美國(guó)密歇根大學(xué)提出多級(jí)瞬態(tài)方法(TML)并應(yīng)用在MPACT上[2,4]，TML采用預(yù)估矯正準(zhǔn)靜態(tài)(PCQS)的思想，將點(diǎn)堆中子動(dòng)力學(xué)、CMFD以及輸運(yùn)計(jì)算相耦合，使得在大時(shí)間步長(zhǎng)條件下能獲得較高計(jì)算精度，且在CMFD計(jì)算中使用1G-CMFD加速M(fèi)G-CMFD收斂[5]。西安交通大學(xué)使用簡(jiǎn)化源縮放法(SSSM)解決CMFD固定源計(jì)算中通量幅值與形狀收斂速度不匹配的問(wèn)題，并在PCQS中點(diǎn)堆中子動(dòng)力學(xué)部分使用在線高階拉格朗日插值來(lái)獲得更高精度的反應(yīng)性計(jì)算值[6]。隨著圖形處理器(GPU)在算力、功耗方面的快速發(fā)展，針對(duì)GPU自身硬件特性而開(kāi)發(fā)的并行算法已經(jīng)被大量應(yīng)用于科學(xué)計(jì)算領(lǐng)域。MOC方法具有天然并行特性，國(guó)內(nèi)外已開(kāi)展了基于GPU加速M(fèi)OC計(jì)算的相關(guān)研究[7-14]。本文針對(duì)全隱式差分方法，同時(shí)采用CUDA編程和CMFD加速技術(shù)，實(shí)現(xiàn)基于GPU并行的二維MOC瞬態(tài)計(jì)算。

1 中子輸運(yùn)動(dòng)力學(xué)計(jì)算理論

1.1 中子輸運(yùn)動(dòng)力學(xué)方程

求解中子輸運(yùn)動(dòng)力學(xué)方程的核心思想是對(duì)其進(jìn)行時(shí)間差分[3]，并利用先驅(qū)核方程積分方法將先驅(qū)核方程與中子輸運(yùn)方程進(jìn)行耦合，推導(dǎo)得到一系列與穩(wěn)態(tài)方程形式相當(dāng)?shù)乃矐B(tài)固定源方程(TFSP)，并利用穩(wěn)態(tài)中子輸運(yùn)方程求解方法對(duì)其進(jìn)行逐時(shí)間步求解。中子輸運(yùn)方程和先驅(qū)核方程如式(1)與式(2)所示：

χp,g(r,t)(1-β(r,t))SF(r,t)+Sd,g(r,t))

(1)

k=1,2,…,6

(2)

式中：vg為能群g下的中子速度；φg為位置r處方向?yàn)棣傅闹凶咏峭棵芏?；Ck為緩發(fā)中子先驅(qū)核濃度；Σt,g為宏觀總截面；βk為緩發(fā)中子份額；λk為衰變常量；Ss,g為散射源項(xiàng)；SF為裂變?cè)错?xiàng)；Sd,g為緩發(fā)中子源項(xiàng)；χp,g為瞬發(fā)中子裂變譜。各源項(xiàng)定義如式(3)～(6)所示。

(3)

(4)

(5)

(6)

對(duì)式(1)的時(shí)間項(xiàng)進(jìn)行全隱式差分，可得式(7)：

(7)

(8)

(9)

(10)

(11)

(12)

(13)

1.2 CMFD加速

多群CMFD TFSP方程可由中子擴(kuò)散動(dòng)力學(xué)方程[2,4-5]推導(dǎo)得出，格式如下：

(14)

(15)

對(duì)式(14)使用與推導(dǎo)輸運(yùn)動(dòng)力學(xué)方程相同的時(shí)間差分方法，可得式(16)：

(16)

(17)

(18)

(19)

(20)

Str=AΦ+BFΦ+C

(21)

將各算子代入式(16)可得CMFD TFSP方程：

(M-S)Φ=χFΦ+AΦ+BFΦ+C

(22)

1.3 MOC動(dòng)力學(xué)計(jì)算

由1.1節(jié)可知，2D MOC動(dòng)力學(xué)方程可寫(xiě)為如下形式[2,4,6]：

(23)

其中瞬態(tài)固定源項(xiàng)定義與式(9)相同：

(24)

在原有源項(xiàng)上加上瞬態(tài)固定源項(xiàng)后，式(24)可由穩(wěn)態(tài)MOC求解器逐時(shí)間步求解，且現(xiàn)有的并行加速算法不需作較大改動(dòng)即可套用至動(dòng)力學(xué)計(jì)算上。

2 并行算法設(shè)計(jì)及程序開(kāi)發(fā)

本文基于上述理論，以現(xiàn)有的柵格物理程序ALPHA穩(wěn)態(tài)版本為基礎(chǔ)，調(diào)用其MOC求解器及CMFD求解器來(lái)開(kāi)發(fā)其動(dòng)力學(xué)計(jì)算功能。ALPHA中采用GPU加速CMFD及MOC計(jì)算，MOC計(jì)算中采用特征線-單能群并行策略，CMFD計(jì)算采用紅-黑排序并行求解策略[6-12]，具體計(jì)算流程如圖1所示。

圖1 GPU加速的MOC耦合CMFD求解瞬態(tài)固定源問(wèn)題流程圖Fig.1 TFSP solution procedure for MOC with CMFD acceleration based on GPU

整體計(jì)算流程如下：1) 在CPU端進(jìn)行初始穩(wěn)態(tài)和動(dòng)力學(xué)的初始化計(jì)算，包括幾何預(yù)處理和截面初始化等計(jì)算所需信息；2) 將幾何、截面信息從系統(tǒng)內(nèi)存一次性拷貝到顯存；3) 執(zhí)行CMFD計(jì)算，為MOC迭代提供初值；4) 通過(guò)本征值和通量偏差判斷整體迭代是否收斂。

CMFD計(jì)算主要包括：1) 粗網(wǎng)多群參數(shù)歸并及線性系統(tǒng)的構(gòu)造，由于各粗網(wǎng)參數(shù)和CMFD方程組中各方程相互獨(dú)立，可將粗網(wǎng)參數(shù)歸并與方程構(gòu)造分配給各GPU線程并行執(zhí)行；2) CMFD源項(xiàng)更新也可在GPU端并行執(zhí)行，每個(gè)線程負(fù)責(zé)單個(gè)粗網(wǎng)單個(gè)能群的源項(xiàng)計(jì)算；3) CMFD線性系統(tǒng)求解，本文所采用的線性方程組解法為逐次超松弛迭代(SOR)，為保證在GPU并行求解，本文采用紅-黑排序，將相互獨(dú)立的方程求解分配給不同的GPU線程并行執(zhí)行，因此求解過(guò)程中有1/2的方程可并行求解，隨后提供更新后的結(jié)果，用于另一半方程的并行求解；4) CMFD收斂后，利用CMFD結(jié)果更新MOC細(xì)網(wǎng)通量，為下一次MOC計(jì)算提供初值。采用GPU線程并行后，1個(gè)線程執(zhí)行1個(gè)粗網(wǎng)內(nèi)所有細(xì)網(wǎng)在單個(gè)能群下的通量更新。

MOC將采用CMFD更新后的通量進(jìn)行計(jì)算。本文采用Jacobi格式的輸運(yùn)掃描算法進(jìn)行輸運(yùn)求解，其中源項(xiàng)計(jì)算采用GPU多線程并行求解，輸運(yùn)掃描計(jì)算采用的特征線-能群并行的策略。

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

本文采用TWIGL基準(zhǔn)題與2D MINI-CORE基準(zhǔn)題測(cè)試程序精度及效率。計(jì)算平臺(tái)采用多核工作站，工作站所用CPU為Intel Core i7-8700K，所用GPU為NVIDIA Geforce 2080Ti。

3.1 TWIGL基準(zhǔn)題

TWIGL基準(zhǔn)題[15]是一個(gè)兩群三區(qū)問(wèn)題，其幾何結(jié)構(gòu)及堆芯網(wǎng)格劃分如圖2所示。TWIGL堆芯內(nèi)部由3種不同材料區(qū)域構(gòu)成，通過(guò)改變區(qū)域1材料的熱吸收截面來(lái)驅(qū)動(dòng)堆芯功率發(fā)生瞬態(tài)變化。ALPHA的計(jì)算參數(shù)如下：柵元尺寸2.0 cm×2.0 cm，每個(gè)柵元使用5×5矩形網(wǎng)格劃分，特征線間距0.02 cm，角度離散采用32個(gè)輻角和最佳三極角，初始穩(wěn)態(tài)和動(dòng)力學(xué)計(jì)算通量收斂準(zhǔn)則均為1.0×10-7，時(shí)間步長(zhǎng)分別選用2、5和10 ms，CPU和GPU計(jì)算均使用雙精度。

圖2 TWIGL基準(zhǔn)題幾何結(jié)構(gòu)及堆芯網(wǎng)格劃分Fig.2 TWIGL benchmark geometry and core division strategy

堆芯功率時(shí)變計(jì)算結(jié)果如圖3、4所示，參考解由VARIANT和DeCART給出[15]，可看出ALPHA計(jì)算結(jié)果與參考解符合良好，堆芯功率最大相對(duì)誤差為0.43%。各區(qū)域平均棒功率計(jì)算結(jié)果比較列于表1。區(qū)域3功率相對(duì)誤差較區(qū)域1、2大，最大誤差出現(xiàn)在0.2 s，恰好是功率開(kāi)始瞬變的時(shí)刻。帶有CMFD加速的不同硬件條件下計(jì)算效率與精度比較列于表2。由于TWIGL問(wèn)題規(guī)模較小，GPU算例中對(duì)于每個(gè)時(shí)間步的計(jì)算時(shí)間并不能達(dá)到所用GPU的峰值算力，且調(diào)用GPU過(guò)程中所消耗計(jì)算資源占比較多，所以GPU算例計(jì)算時(shí)間并未隨著步長(zhǎng)增大而呈比例地減小，但相比CPU算例而言，加速效果依然明顯。對(duì)于2 ms步長(zhǎng)算例，加速比達(dá)到6.6。GPU并行條件下的CMFD加速效率對(duì)比列于表3。CMFD加速性能十分可觀，加速比基本在40左右。對(duì)比表3中3種不同的時(shí)間步長(zhǎng)計(jì)算結(jié)果可發(fā)現(xiàn)，CMFD時(shí)間占比隨時(shí)間步長(zhǎng)的減小而減小，從而CMFD加速比也隨之減小。對(duì)比表2與表3結(jié)果可發(fā)現(xiàn)，GPU并行條件下的加速比隨著MOC計(jì)算時(shí)間占比增大而增大。

圖3 TWIGL功率變化及相對(duì)誤差Fig.3 Core power behavior of TWIGL and relative error

圖4 不同時(shí)刻TWIGL精細(xì)功率分布 Fig.4 Detailed power distribution of TWIGL at different time

表1 TWIGL各區(qū)域平均棒功率計(jì)算結(jié)果比較Table 1 Region wise pin power comparison of TWIGL

表2 TWIGL基準(zhǔn)題不同硬件條件下計(jì)算效率及精度比較Table 2 Computational efficiency and accuracy comparison for TWIGL in different hardware conditions

表3 TWIGL基準(zhǔn)題CMFD加速效率比較Table 3 Effectiveness of CMFD for TWIGL

3.2 MINI-CORE 2D基準(zhǔn)題

MINI-CORE 2D基準(zhǔn)題[15]幾何結(jié)構(gòu)以及邊界條件如圖5所示，該堆芯由5個(gè)鈾富集度為4.2%的UOX組件和4個(gè)钚富集度為4.0%的MOX組件組成。控制棒僅位于中心的UOX組件內(nèi)，通過(guò)勻速提棒使得堆芯產(chǎn)生瞬態(tài)功率變化。ALPHA的計(jì)算參數(shù)如下：每個(gè)組件為17×17柵元排布，每個(gè)柵元使用5×5矩形網(wǎng)格劃分，選用特征線間距為0.02 cm，角度離散選用32輻角和最佳三極角，初始穩(wěn)態(tài)和動(dòng)力學(xué)計(jì)算通量收斂準(zhǔn)則均為1.0×10-6，CPU和GPU計(jì)算均使用雙精度。

圖5 MINI-CORE 2D基準(zhǔn)題幾何結(jié)構(gòu)Fig.5 MINI-CORE 2D benchmark geometry

堆芯功率時(shí)變計(jì)算結(jié)果如圖6、7所示，參考解由DeCART與VARIANT給出[15]，ALPHA計(jì)算結(jié)果與參考解符合良好，2 ms步長(zhǎng)算例中堆芯功率最大相對(duì)誤差為3.3%，出現(xiàn)位置在0.01 s。不同時(shí)刻下組件功率計(jì)算結(jié)果列于表4，各時(shí)刻下組件功率相對(duì)誤差均在0.05%以內(nèi)。帶有CMFD加速的不同硬件條件下計(jì)算效率與精度比較列于表5，可看出GPU計(jì)算加速效果良好，對(duì)于2 ms算例加速比為3.8倍。GPU條件下的CMFD加速效率對(duì)比列于表6，CMFD加速效果明顯，10 ms算例加速比可達(dá)到51.6。MINI-CORE基準(zhǔn)題中GPU加速比與CMFD時(shí)間占比的關(guān)系與TWIGL結(jié)果一致，即GPU加速比隨CMFD時(shí)間占比增大而減小。

圖6 MINI-CORE 2D功率變化及相對(duì)誤差Fig.6 Core power behavior of MINI-CORE 2D and relative error

圖7 不同時(shí)刻MINI-CORE 2D精細(xì)功率分布 Fig.7 Detailed power distribution of MINI-CORE 2D at different time

表4 MINI-CORE 2D各區(qū)域平均棒功率計(jì)算結(jié)果比較Table 4 Region wise pin power comparison of MINI-CORE 2D

表5 MINI-CORE 2D基準(zhǔn)題keff及計(jì)算效率比較Table 5 Steady-state keff and calculation efficiency of MINI-CORE 2D

表6 MINI-CORE基準(zhǔn)題CMFD加速效率比較Table 6 Effectiveness of CMFD for MINI-CORE

4 結(jié)論

本文介紹了采用GPU加速的TFSP并行算法，實(shí)現(xiàn)了GPU并行加速動(dòng)力學(xué)MOC及CMFD計(jì)算。TWIGL基準(zhǔn)題和MINI-CORE 2D基準(zhǔn)題數(shù)值結(jié)果表明：程序滿足瞬態(tài)輸運(yùn)求解的精度要求；同時(shí)，采用單GPU進(jìn)行小規(guī)模堆芯瞬態(tài)輸運(yùn)計(jì)算時(shí)，與單CPU相比加速倍數(shù)約在2～6倍，且加速比隨CMFD計(jì)算時(shí)間占比的增大而減小。本文僅初步實(shí)現(xiàn)了CMFD TFSP加速計(jì)算，所涉及的并行算法可推廣至預(yù)估-校正準(zhǔn)靜態(tài)方法以及TML方法，提升全堆瞬態(tài)計(jì)算的效率。