一種全堆芯精確到每個(gè)通道的子通道并行模擬方法

王先夢(mèng)，趙民富，呂玉鳳，蔡銀宇，儲(chǔ)根深，盧 旭，王昭順，郭蘇萱，周志鋒，胡長(zhǎng)軍，楊 文

(1.北京科技大學(xué) 計(jì)算機(jī)與通信工程學(xué)院，北京 100083；2.中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413)

子通道分析方法是反應(yīng)堆熱工水力分析中常用的方法。子通道是人為假想的冷卻劑流道，每個(gè)子通道沿軸向分為若干層，每層被稱為1個(gè)控制體，控制體是最小計(jì)算單元，無論控制體的體積大小，同一控制體內(nèi)的冷卻劑壓力、溫度和速度等熱工參數(shù)均被認(rèn)為是相同的。子通道的設(shè)定影響了計(jì)算的精細(xì)程度。為減少計(jì)算量和存儲(chǔ)量，傳統(tǒng)的串行子通道分析程序一般采用簡(jiǎn)化的、粗略的幾何模型，如不計(jì)算全堆芯而只計(jì)算部分子通道，或?qū)?個(gè)組件內(nèi)多個(gè)燃料棒圍成的區(qū)域、1個(gè)組件、甚至多個(gè)組件視為1個(gè)子通道。這種粗略的幾何建模方式引入過多的假設(shè)和近似，已越來越不能滿足研究和應(yīng)用的需求，高精細(xì)的子通道分析受到國內(nèi)外學(xué)者的廣泛重視[1-5]。

Geemert等[6]在F-COBRA-TF的求解壓力方程組的矩陣-向量乘循環(huán)中加入OpenMP[7]編譯指導(dǎo)語句來實(shí)現(xiàn)程序并行化。這種并行策略易于實(shí)現(xiàn)，但因其只適用于共享內(nèi)存的計(jì)算系統(tǒng)，導(dǎo)致其可擴(kuò)展性非常有限。Chaudri等[8]在START程序的多個(gè)循環(huán)中加入OpenMP來并行化，并行后的START能擴(kuò)展至40核，并行效率為78%。在美國CASL數(shù)值堆項(xiàng)目[9]中，Salko等[1,10-13]對(duì)CTF子通道軟件進(jìn)行了優(yōu)化，解決了內(nèi)存使用量大、預(yù)處理過程繁瑣耗時(shí)等問題。文獻(xiàn)[1，11]采用基于區(qū)域分解的并行策略將CTF并行化，每個(gè)組件被劃分到1個(gè)求解域，這種并行策略要求進(jìn)程數(shù)與組件數(shù)相等，限制了程序的可擴(kuò)展性和靈活性。Ramos等[14]開展了針對(duì)CTF軟件的并行化研究，兩個(gè)計(jì)算最耗時(shí)的部分被并行化，其中的線性方程組通過使用數(shù)學(xué)求解庫PETSc[15]來并行化。明平洲等[16-17]對(duì)CORTH軟件進(jìn)行了并行化，實(shí)現(xiàn)了全堆芯計(jì)算，結(jié)果顯示并行化后的CORTH可擴(kuò)展到36核并獲得8倍加速比。Misawa等[18]采用區(qū)域分解的方式對(duì)ACE-3D軟件進(jìn)行并行化，ACE-3D能擴(kuò)展到252核，以63核為比較基準(zhǔn)，獲得83.1%的并行效率。

目前，最精細(xì)的子通道分析是進(jìn)行全堆芯的、精確到每個(gè)通道的子通道計(jì)算。高精細(xì)的子通道分析需大量的計(jì)算和存儲(chǔ)資源，需采用并行計(jì)算才能高效實(shí)現(xiàn)。并行劃分策略是實(shí)現(xiàn)并行系統(tǒng)的關(guān)鍵，與并行系統(tǒng)的效率和可擴(kuò)展性密切相關(guān)。本文提出一種基于子通道的并行任務(wù)劃分和進(jìn)程映射方法。

1 子通道模型

精確到每個(gè)通道的子通道幾何建模，分為冷卻劑為中心和燃料棒為中心兩種方式[19]。如圖1的藍(lán)色部分所示，冷卻劑為中心的子通道是由4根燃料棒圍成的流動(dòng)區(qū)域。圖1的黃色部分為以燃料棒為中心的子通道。由于冷卻劑為中心的子通道使用更廣泛，本文采用冷卻劑為中心的子通道建模方式。

圖1 冷卻劑為中心和燃料棒為中心的子通道Fig.1 Coolant-centered and rod-centered subchannels

全堆芯、精確到每個(gè)通道的壓水堆子通道建模如圖2所示。壓水堆堆芯一般由上百個(gè)燃料組件組成，每個(gè)燃料組件由幾百根燃料棒組成(圖2中橙色部分)。燃料棒被冷卻劑(圖2中藍(lán)色部分)圍繞。每個(gè)子通道從軸向被劃分為幾十到幾百不等若干層。每一層被稱為1個(gè)控制體或網(wǎng)格，每個(gè)控制體內(nèi)質(zhì)量、動(dòng)量、能量守恒。

圖2 子通道幾何建模Fig.2 Geometrical model of subchannel

2 基于組件均分的并行劃分方法

并行任務(wù)劃分將整體求解問題分解為若干個(gè)子問題，子問題被多個(gè)進(jìn)程同時(shí)處理。本文提出基于區(qū)域分解的組件均分法進(jìn)行并行任務(wù)劃分，按進(jìn)程數(shù)少于組件數(shù)和進(jìn)程數(shù)大于組件數(shù)兩種情況分別考慮。

2.1 進(jìn)程數(shù)少于組件數(shù)的劃分策略

用戶的計(jì)算資源為單機(jī)或小型服務(wù)器時(shí)，目前市場(chǎng)上的小型服務(wù)器的核數(shù)一般少于64，此時(shí)用戶設(shè)定的進(jìn)程數(shù)少于組件數(shù)。以將全堆芯12組件劃分為4求解域?yàn)槔?，圖3示出組件均分法在進(jìn)程數(shù)少于組件數(shù)時(shí)的劃分過程。該堆芯中每個(gè)組件包含3×3棒束。劃分過程如下：1) 引入假想燃料棒和假想子通道，將整個(gè)求解域擴(kuò)展為1個(gè)正方形(稱作堆芯補(bǔ)全圖)，圖3中淺灰色圓形表示假想燃料棒，淺灰色正方形表示假想子通道；2) 將真實(shí)存在的棒/子通道、假想的棒/子通道用不同的標(biāo)志位分別標(biāo)志；3) 將用戶設(shè)置的進(jìn)程數(shù)拆解成兩個(gè)最接近的整數(shù)相乘，較小的作為行數(shù)，較大的作為列數(shù)；4) 將整個(gè)堆芯按行列均分，得到進(jìn)程數(shù)個(gè)子區(qū)域，如果該區(qū)域內(nèi)存在真實(shí)的棒/子通道，該子區(qū)域有效，否則該子區(qū)域無效，如果存在無效區(qū)域，最終得到的求解域數(shù)將少于用戶設(shè)置的進(jìn)程數(shù)；5) 按行優(yōu)先的順序，將有效子區(qū)域依次編號(hào)，如圖3，有效子區(qū)域依次編號(hào)為1～4；6) 使用有效子區(qū)域編號(hào)標(biāo)志燃料棒和子通道，如子區(qū)域1內(nèi)的棒和子通道均標(biāo)志為1，至此，求解域劃分完成。

圖3 將全堆芯劃分為4個(gè)求解域的過程Fig.3 Process of partitioning full-core into 4 domains

2.2 進(jìn)程數(shù)大于組件數(shù)的劃分策略

圖4 全堆芯121組件劃分為484進(jìn)程的過程Fig.4 Process of partitioning 121 assemblies of full-core into 484 processes

當(dāng)用戶的計(jì)算資源是大型集群，甚至超級(jí)計(jì)算機(jī)時(shí)，可將整個(gè)堆芯劃分為幾百、甚至幾千個(gè)求解域，以充分利用計(jì)算資源。參考Salko等[10-13]將每個(gè)組件劃分為1個(gè)求解域的劃分方式，當(dāng)進(jìn)程數(shù)大于組件數(shù)時(shí)，劃分以組件為基準(zhǔn)進(jìn)行。整個(gè)堆芯可劃分為組件數(shù)的整數(shù)倍，如圖4所示。劃分過程如下：1) 讀入用戶設(shè)置的組件倍數(shù)，圖4中，組件數(shù)為121，組件倍數(shù)為4，進(jìn)程數(shù)為484；2) 將組件倍數(shù)拆成最小的兩個(gè)正整數(shù)相乘，較小的作為行數(shù)，較大的作為列數(shù)；3) 每個(gè)組件按行列均分，整個(gè)堆芯被劃分成組件數(shù)乘以倍數(shù)個(gè)子區(qū)域；4) 以行優(yōu)先的方式，依次對(duì)子區(qū)域進(jìn)行編號(hào)；5) 使用子區(qū)域的編號(hào)標(biāo)志燃料棒和子通道，如子區(qū)域1內(nèi)的燃料棒和子通道都標(biāo)志為1；6) 子區(qū)域編號(hào)即為求解域編號(hào)，至此求解域劃分完成。

2.3 劃分策略的實(shí)現(xiàn)

1) 自動(dòng)編號(hào)及數(shù)據(jù)結(jié)構(gòu)

劃分策略的實(shí)現(xiàn)，需建立子通道/棒編號(hào)、子通道/棒與組件對(duì)應(yīng)關(guān)系以及子通道與棒映射關(guān)系。根據(jù)用戶輸入的組件分布信息，按行優(yōu)先對(duì)組件連續(xù)編號(hào)，無組件的位置編號(hào)為0。圖5示出全堆芯57組件的組件編號(hào)圖，其中綠色位置表示存在組件，灰色位置表示無組件。

圖5 組件編號(hào)圖Fig.5 Numbered layout of assembly

全堆芯57組件，每組件3×3燃料棒圖6 棒的所屬組件標(biāo)志圖Fig.6 Identified layout of rod

根據(jù)組件編號(hào)圖和用戶輸入的棒分布信息，建立棒的所屬組件標(biāo)志圖，如圖6所示。棒的所屬組件標(biāo)志圖給出每個(gè)棒所屬組件索引，主要用于進(jìn)程數(shù)大于組件數(shù)時(shí)區(qū)分不同的組件。

建立通道的所屬關(guān)系圖，采取行優(yōu)先遍歷方式，設(shè)置每個(gè)棒周圍4個(gè)通道的所屬組件值為棒的所屬組件值(如某通道已設(shè)置則跳過該通道)。棒/通道的所屬組件圖是進(jìn)行任務(wù)劃分的依據(jù)。圖7示出子通道的所屬組件標(biāo)志圖。

圖7 子通道的所屬組件標(biāo)志圖Fig.7 Identified layout of subchannel

建立子通道全局編號(hào)圖，按行優(yōu)先對(duì)子通道所屬組件編號(hào)圖中的子通道連續(xù)編號(hào)。圖8示出子通道的全局編號(hào)圖。棒的全局編號(hào)圖建立方式與子通道的全局編號(hào)圖類似。

所有編號(hào)采用廣義表存儲(chǔ)。定義LS為廣義表，LS=(No_assem，IdAssem_rod，IdAssem_chan，No_rod，No_chan)。其中，子表No_assem存儲(chǔ)No_assem，子表IdAssem_rod和IdAssem_chan分別存儲(chǔ)棒的所屬組件編號(hào)圖和子通道的所屬組件編號(hào)圖，子表No_rod和No_chan分別存儲(chǔ)棒和子通道的全局編號(hào)圖。

2) 區(qū)域分解

當(dāng)進(jìn)程數(shù)小于組件數(shù)時(shí)，區(qū)域分解流程圖如圖9a所示。當(dāng)進(jìn)程數(shù)大于組件數(shù)時(shí)，區(qū)域分解流程圖如圖9b所示。

3) 求解域邊界處理

每個(gè)子通道的求解都需要來自其鄰居子通道的信息。進(jìn)行任務(wù)劃分后，處于求解域交界處的子通道，其鄰居子通道可能在其他求解域。為處理求解域交界處子通道從鄰居子通道獲取信息的問題，引入了ghost子通道。ghost子通道不進(jìn)行計(jì)算，只用來與鄰居子通道進(jìn)行通信。需要與鄰居子通道交換信息時(shí)，每個(gè)求解域先遍歷本地?fù)碛術(shù)host區(qū)域的子通道，使用MPI發(fā)送函數(shù)將本地信息發(fā)送到ghost區(qū)域。發(fā)送完成后，每個(gè)求解域遍歷自己的ghost區(qū)域，接收其他進(jìn)程發(fā)送來的信息。建立ghost區(qū)域的方式如下。

圖8 子通道的全局編號(hào)圖Fig.8 Global numbered layout of subchannel

假設(shè)當(dāng)前子通道0所在求解域的值為D0，其周圍4個(gè)相鄰位置的子通道1～4所在的求解域編號(hào)分別為D1～D4。如果Di≠D0，i=1，2，3，4，則稱該通道為求解域D0的ghost子通道。

圖10示出求解域交界處子通道的示意圖。假設(shè)將1個(gè)堆芯(內(nèi)含4盒組件)劃分為16個(gè)求解域(圖10a)，用藍(lán)色框突出顯示的子通道在第13號(hào)求解域，其ghost子通道分別在第7和第10號(hào)求解域(圖10b)。第13號(hào)求解域的ghost區(qū)域如圖11所示，淺藍(lán)色表示本地子通道，深藍(lán)色表示ghost子通道。

4) 混合編程實(shí)現(xiàn)

為充分利用集群的層次性存儲(chǔ)結(jié)構(gòu)特點(diǎn)，獲得良好的并行性能，采用消息傳遞模型和共享內(nèi)存模型混合的并行編程模型實(shí)現(xiàn)[20-21]。進(jìn)行并行任務(wù)劃分后，采用消息傳遞編程模型MPI，每個(gè)進(jìn)程處理1個(gè)求解域。在每個(gè)進(jìn)程內(nèi)，在求解傳熱系數(shù)的部分，遍歷控制體的循環(huán)無相關(guān)性；在求解燃料棒的熱傳導(dǎo)方程部分，遍歷燃料棒的循環(huán)無相關(guān)性，在這幾部分應(yīng)用共享內(nèi)存編程模型OpenMP對(duì)程序進(jìn)行線程級(jí)的并行化處理。圖12示出混合編程示意圖。

a——進(jìn)程數(shù)小于組件數(shù)；b——進(jìn)程數(shù)大于組件數(shù)圖9 區(qū)域分解流程圖Fig.9 Flow chart of domain decomposition algorithm

圖10 求解域交界處子通道Fig.10 Subchannel in boundary of domain

圖11 求解域13及其ghost子通道Fig.11 Domain thirteen and its ghost channel

3 并行劃分策略測(cè)試

3.1 算例描述

共采用3個(gè)測(cè)試算例。算例1為全堆芯157組件穩(wěn)態(tài)工況模擬，每組件內(nèi)有17×17棒束，輸入設(shè)置來自文獻(xiàn)[22-24]。算例1的幾何和熱工參數(shù)列于表1。測(cè)試分為3個(gè)規(guī)模：小規(guī)模測(cè)試，軸向分21層，控制體數(shù)目約為100萬；中等規(guī)模測(cè)試，軸向分為125層，控制體數(shù)目約為600萬；大規(guī)模測(cè)試，軸向分為249層，控制體數(shù)目約為1 200萬。

圖12 混合編程示意圖Fig.12 Diagram of hybrid programming

表1 熱工和幾何參數(shù)Table 1 Thermal parameter and geometric parameter

算例2為全堆芯49組件，每組件9×9棒模擬。算例3為全堆芯49組件，每組件11×11棒模擬。算例2、3除組件和燃料棒數(shù)目外，其余輸入與算例1的一致。

3.2 測(cè)試結(jié)果及分析

1) 可擴(kuò)展性測(cè)試

使用算例1的中等規(guī)模測(cè)試在天河二號(hào)超級(jí)計(jì)算機(jī)上進(jìn)行可擴(kuò)展性測(cè)試。表2列出使用不同進(jìn)程時(shí)的總運(yùn)行時(shí)間、通信時(shí)間和通信時(shí)間占總運(yùn)行時(shí)間的比例(通信占比)。進(jìn)程數(shù)較少時(shí)，每個(gè)進(jìn)程上的計(jì)算量都較大，所以進(jìn)程數(shù)的增加會(huì)使單個(gè)進(jìn)程的計(jì)算量顯著下降，程序性能提升較為明顯。使用4進(jìn)程進(jìn)行全堆芯模擬需約22 h，使用1 884進(jìn)程，時(shí)間縮短至417 s。隨著進(jìn)程數(shù)的進(jìn)一步增加，程序性能反而降低，使用6 280進(jìn)程時(shí)，模擬時(shí)間增加到470 s。由于子通道方法的計(jì)算量有限，隨進(jìn)程數(shù)的增多，通信占比不斷增大。使用4進(jìn)程時(shí)，通信占比僅2.37%，使用6 280進(jìn)程時(shí)，通信占比已達(dá)到34.84%。

表2 全堆芯157組件中等規(guī)模測(cè)試模擬時(shí)間Table 2 Simulation time of full-core with 157 assemblies

進(jìn)一步使用算例1在另一集群進(jìn)行測(cè)試，如圖13a所示。在小規(guī)模測(cè)試中，使用942進(jìn)程時(shí)運(yùn)行時(shí)間最短(小于1 min)，與4進(jìn)程的基準(zhǔn)(運(yùn)行約1 h)相比，運(yùn)行時(shí)間快約59倍，處理器數(shù)目超過942時(shí)，運(yùn)行時(shí)間增加。在中等規(guī)模測(cè)試中，使用942進(jìn)程時(shí)運(yùn)行時(shí)間也最短(約10 min)，與4進(jìn)程的基準(zhǔn)相比(約23 h)，運(yùn)行時(shí)間加速約143倍。在大規(guī)模測(cè)試中，使用1 884進(jìn)程時(shí)運(yùn)行時(shí)間最短(約32 min)，與4進(jìn)程的基準(zhǔn)相比(約84 h)，運(yùn)行時(shí)間大約快158倍。算例2、3的測(cè)試結(jié)果也呈現(xiàn)出類似規(guī)律，如圖13b、c所示。

2) 混合編程結(jié)果

使用算例1的中等規(guī)模測(cè)試在天河二號(hào)超級(jí)計(jì)算機(jī)上進(jìn)行混合編程結(jié)果測(cè)試。在純MPI版本中，MPI的數(shù)量等于求解域的數(shù)量，每個(gè)MPI進(jìn)程在1個(gè)物理核上運(yùn)行并處理1個(gè)求解域，測(cè)試時(shí)進(jìn)程數(shù)從1 256擴(kuò)展到6 280。在混合編程版本中，每個(gè)MPI進(jìn)程創(chuàng)建2個(gè)線程，每個(gè)MPI進(jìn)程運(yùn)行在2個(gè)物理核上。每個(gè)進(jìn)程處理1個(gè)求解域，因此每個(gè)求解域由2個(gè)線程處理，進(jìn)程數(shù)從628擴(kuò)展到3 140，相應(yīng)的線程數(shù)則從1 256擴(kuò)展到6 280。混合編程與純MPI版本在使用不同進(jìn)程數(shù)/線程數(shù)測(cè)試時(shí)的總運(yùn)行時(shí)間和通信時(shí)間列于表3。

a——全堆芯157組件，每組件17×17棒束；b——全堆芯49組件，每組件9×9燃料棒；c——全堆芯49組件，每組件11×11 棒束圖13 可擴(kuò)展性測(cè)試Fig.13 Scalability test

表3 混合編程與純MPI版本的總運(yùn)行時(shí)間和通信時(shí)間對(duì)比Table 3 Comparison of total running time and communication time of hybrid programming and pure MPI version

由表4可見，混合編程帶來的性能提升并不明顯，特別是當(dāng)進(jìn)程數(shù)較少時(shí)(如進(jìn)程數(shù)少于2 512時(shí))，多線程版本比純MPI版本更耗時(shí)。這是因?yàn)榧尤隣penMP的部分(求解傳熱系數(shù)和導(dǎo)熱方程兩部分)在外迭代的內(nèi)部，每個(gè)時(shí)間步都會(huì)產(chǎn)生產(chǎn)生和撤銷線程的開銷。另外在求解傳熱系數(shù)的部分，需要先求解橫向質(zhì)量流量，而求解每個(gè)子通道的質(zhì)量流量需考慮與周圍所有相鄰子通道的質(zhì)量交換，因此需與鄰居子通道通信。求解傳熱系數(shù)部分需使用多個(gè)OpenMP并行區(qū)域，導(dǎo)致OpenMP帶來的額外開銷較大。隨進(jìn)程數(shù)的增加，混合編程性能稍好于純MPI版本，這是通過減少通信時(shí)間獲得的。使用同樣核數(shù)時(shí)，混合編程的MPI數(shù)量是純MPI版本數(shù)量的一半，混合編程中使用的求解域的數(shù)量是純MPI版本的一半，通信只發(fā)生在進(jìn)程間，混合編程的通信時(shí)間基本都少于純MPI版本的通信時(shí)間。

3) 程序初步驗(yàn)證

采用全堆芯157組件穩(wěn)態(tài)模擬算例，通過比較計(jì)算輸出的堆芯壓降來分析并行程序的正確性，結(jié)果列于表4。由表4可見，串行計(jì)算與使用18、68、157、942進(jìn)程并行計(jì)算所得到的堆芯壓降基本一致，與文獻(xiàn)[23]描述的900 MW壓水堆公布的堆芯壓降0.156 6 MPa吻合較好。

表4 不同進(jìn)程數(shù)計(jì)算輸出的堆芯壓降Table 4 Calculation pressure loss using different processes

4 小結(jié)

為實(shí)現(xiàn)全堆芯、精確到每個(gè)真實(shí)通道的子通道分析，提出了組件均分法進(jìn)行并行任務(wù)劃分。實(shí)驗(yàn)結(jié)果顯示，該方法可擴(kuò)展至上千核。全堆芯157組件，軸向分125層的穩(wěn)態(tài)工況模擬僅需不到10 min即可完成。采用全堆芯49組件，每組件9×9、11×11棒，軸向分層分別為21、125、249的算例進(jìn)行分析，發(fā)現(xiàn)當(dāng)進(jìn)程數(shù)為組件數(shù)6～20倍時(shí)，程序運(yùn)行時(shí)間最短。當(dāng)進(jìn)程數(shù)逐漸增大時(shí)，運(yùn)行時(shí)間反而上升。引入消息傳遞模型和共享內(nèi)存模型混合的編程模型實(shí)現(xiàn)后發(fā)現(xiàn)，當(dāng)進(jìn)程數(shù)較少時(shí)，純MPI版本性能更優(yōu)，隨進(jìn)程數(shù)增多，MPI和OpenMP混合編程性能優(yōu)于純MPI版本。目前的并行策略根據(jù)壓水堆的幾何結(jié)構(gòu)設(shè)計(jì)，下一步計(jì)劃根據(jù)快堆的幾何結(jié)構(gòu)特點(diǎn)設(shè)計(jì)相應(yīng)并行模擬方法。