王先夢，趙民富，呂玉鳳，蔡銀宇，儲根深，盧 旭，王昭順，郭蘇萱，周志鋒，胡長軍，楊 文

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

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

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

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

1 子通道模型

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

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

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

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

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

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

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

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

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

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

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

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

圖5 組件編號圖Fig.5 Numbered layout of assembly

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

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

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

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

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

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

2) 區(qū)域分解

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

3) 求解域邊界處理

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

圖8 子通道的全局編號圖Fig.8 Global numbered layout of subchannel

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

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

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

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

3.1 算例描述

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

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

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

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

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

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

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

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

表3 混合編程與純MPI版本的總運行時間和通信時間對比Table 3 Comparison of total running time and communication time of hybrid programming and pure MPI version

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

3) 程序初步驗證

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

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

4 小結(jié)

為實現(xiàn)全堆芯、精確到每個真實通道的子通道分析，提出了組件均分法進(jìn)行并行任務(wù)劃分。實驗結(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倍時，程序運行時間最短。當(dāng)進(jìn)程數(shù)逐漸增大時，運行時間反而上升。引入消息傳遞模型和共享內(nèi)存模型混合的編程模型實現(xiàn)后發(fā)現(xiàn)，當(dāng)進(jìn)程數(shù)較少時，純MPI版本性能更優(yōu)，隨進(jìn)程數(shù)增多，MPI和OpenMP混合編程性能優(yōu)于純MPI版本。目前的并行策略根據(jù)壓水堆的幾何結(jié)構(gòu)設(shè)計，下一步計劃根據(jù)快堆的幾何結(jié)構(gòu)特點設(shè)計相應(yīng)并行模擬方法。