非結(jié)構(gòu)CFD軟件MPI+OpenMP混合并行及超大規(guī)模非定常并行計(jì)算的應(yīng)用

王年華，常興華，趙鐘，張來(lái)平,

1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000 2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

隨著計(jì)算機(jī)技術(shù)和數(shù)值方法的飛速發(fā)展，計(jì)算流體力學(xué)(Computational Fluid Dynamics,CFD)數(shù)值模擬在航空航天等領(lǐng)域已得到越來(lái)越廣泛的應(yīng)用。經(jīng)過(guò)幾十年的發(fā)展，基于雷諾平均Navier-Stokes(Reynolds Averaged Navier-Stokes,RANS)方程的常規(guī)狀態(tài)氣動(dòng)力/力矩預(yù)測(cè)已經(jīng)難度不大，但是在遇到旋渦、分離、轉(zhuǎn)捩、湍流噪聲、湍流燃燒等非定常、非線性現(xiàn)象明顯的流動(dòng)時(shí)，在千萬(wàn)量級(jí)網(wǎng)格上求解RANS方程無(wú)法得到足夠精確的數(shù)值解，這時(shí)就需要采用更大規(guī)模的網(wǎng)格，采用更高保真度的數(shù)值方法，如大渦模擬(Large Eddy Simulation,LES)和直接數(shù)值模擬(Direct Numerical Simulation,DNS)等。這些方法的共同特點(diǎn)就是對(duì)網(wǎng)格量要求很高，通常認(rèn)為L(zhǎng)ES方法在黏性底層內(nèi)對(duì)網(wǎng)格量的需求達(dá)到Re1.8量級(jí)，而DNS則要求網(wǎng)格量達(dá)到Re9/4。對(duì)于實(shí)際飛行器外形，雷諾數(shù)Re通常在百萬(wàn)量級(jí)以上，那么網(wǎng)格量則至少要達(dá)到百億量級(jí)以上，才能滿(mǎn)足算法對(duì)多尺度流動(dòng)結(jié)構(gòu)的分辨率要求[1-3]。即便是采用高階精度格式，雖說(shuō)網(wǎng)格量可以適當(dāng)減少，但是所需的網(wǎng)格量仍需以?xún)|計(jì)。

另一方面，E級(jí)計(jì)算時(shí)代即將來(lái)臨，計(jì)算機(jī)集群核心數(shù)已經(jīng)突破千萬(wàn)核，但是目前主流的CFD工程應(yīng)用仍然停留在千萬(wàn)量級(jí)網(wǎng)格數(shù)百核并行計(jì)算規(guī)模，CFD軟件的發(fā)展與計(jì)算機(jī)硬件的發(fā)展速度極不匹配，這也直接制約了CFD精細(xì)數(shù)值模擬在實(shí)際工程問(wèn)題中的應(yīng)用。采用大型計(jì)算機(jī)集群，調(diào)用超大規(guī)模計(jì)算資源進(jìn)行并行計(jì)算，同時(shí)提高并行效率，發(fā)揮出超級(jí)計(jì)算機(jī)的最佳性能，達(dá)到較高的效費(fèi)比，是未來(lái)高分辨率數(shù)值模擬必須解決的重要問(wèn)題，同時(shí)也是E級(jí)計(jì)算時(shí)代到來(lái)之前所必須解決的軟硬件匹配發(fā)展的問(wèn)題[4]。

近年來(lái)，既有分布內(nèi)存又有共享內(nèi)存的多核處理器成為目前高性能計(jì)算機(jī)的主流計(jì)算核心芯片[5]。“天河二號(hào)”(2018年11月世界TOP500排名第4)采用16 000個(gè)計(jì)算節(jié)點(diǎn)，每個(gè)計(jì)算節(jié)點(diǎn)裝備了2顆12核的Ivy Bridge-E Xeon E5 2692處理器和3顆57核的Xeon Phi加速卡，2個(gè) 處理器共享64 GB內(nèi)存，同時(shí)每個(gè)Xeon Phi加速卡板載有8 GB內(nèi)存，計(jì)算節(jié)點(diǎn)間通過(guò)自制的TH Express-2私有高速網(wǎng)絡(luò)互連，MPI通信速率可達(dá)6.36 Gbps[6-7]。“神威太湖之光”(2018年11月 世界TOP500排名第3)采用40 960個(gè)國(guó)產(chǎn)“申威26010”眾核處理器，每個(gè)處理器有260個(gè)核心，每個(gè)處理器固化的板載內(nèi)存為32 GB[8-9]。對(duì)于這種節(jié)點(diǎn)內(nèi)共享內(nèi)存、節(jié)點(diǎn)間通過(guò)高速網(wǎng)絡(luò)互連的架構(gòu)來(lái)說(shuō)，要減少非計(jì)算時(shí)間占比，盡可能提高程序的并行度，就需要充分利用硬件特性，減少通信次數(shù)、通信量，提高并行效率。

傳統(tǒng)的并行模式在節(jié)點(diǎn)間采用MPI消息傳遞機(jī)制進(jìn)行通信，這種模式在規(guī)模較小時(shí)能夠有效擴(kuò)大并行規(guī)模，并且保持較高的并行效率。但是隨著并行規(guī)模增大，分區(qū)數(shù)增多，通信量和通信次數(shù)呈幾何級(jí)數(shù)增大，導(dǎo)致在大規(guī)模并行時(shí)，并行效率急劇下降。為了充分利用現(xiàn)有計(jì)算節(jié)點(diǎn)共享內(nèi)存的特性，可以在節(jié)點(diǎn)內(nèi)部采用共享內(nèi)存模式進(jìn)行并行計(jì)算。

本文在原有MPI并行模式的基礎(chǔ)上，在課題組自研的基于非結(jié)構(gòu)網(wǎng)格的二階精度有限體積CFD軟件HyperFLOW中，引入共享內(nèi)存并行模式，并通過(guò)OpenMP改造，實(shí)現(xiàn)了MPI+OpenMP兩級(jí)混合并行。首先，根據(jù)節(jié)點(diǎn)內(nèi)并行數(shù)據(jù)粒度大小，將混合并行分為粗粒度并行和細(xì)粒度并行，簡(jiǎn)單介紹了兩種并行的實(shí)現(xiàn)原理和基于C++編程語(yǔ)言的實(shí)現(xiàn)細(xì)節(jié)。其次，在國(guó)產(chǎn)in-house集群上，通過(guò)CRM(Common Research Model)標(biāo)模定常湍流算例[10]對(duì)兩種混合并行模式進(jìn)行測(cè)試和比較。隨后，為了驗(yàn)證兩種混合并行模式在非定常計(jì)算中的可擴(kuò)展性，在機(jī)翼外掛物投放標(biāo)模算例的3.6億非結(jié)構(gòu)重疊網(wǎng)格上進(jìn)行效率測(cè)試，并采用12 288核完成了基于混合并行模式投放過(guò)程的非定常計(jì)算，得到了較高的并行計(jì)算效率。最后，在in-house集群上采用4.9萬(wàn)核，在28.8億非結(jié)構(gòu)重疊網(wǎng)格上進(jìn)行效率測(cè)試，驗(yàn)證了混合并行改造的效果。

1 HyperFLOW軟件簡(jiǎn)介

HyperFLOW軟件平臺(tái)[11-12]是作者團(tuán)隊(duì)自主發(fā)展的結(jié)構(gòu)和非結(jié)構(gòu)混合CFD解算器。該解算器采用C++面向?qū)ο蟮脑O(shè)計(jì)思路，建立了“運(yùn)行數(shù)據(jù)庫(kù)”“模塊化求解”等先進(jìn)的軟件架構(gòu)，具有良好的通用性和可擴(kuò)展性。其可以進(jìn)行結(jié)構(gòu)解算器和非結(jié)構(gòu)解算器的獨(dú)立求解，也可實(shí)現(xiàn)兩種解算器在空間上的耦合求解。解算器的空間離散采用了二階精度的格心型有限體積方法，集成了Roe、vanLeer、AUSM、Steger-Warming等通量計(jì)算格式以及Barth、Venkatakrishnan等多種限制器函數(shù)。湍流模擬方法有SA一方程湍流模型、SST兩方程湍流模型以及DES方法等。時(shí)間離散格式有Runge-Kutta、隱式BLU-SGS等方法，非定常計(jì)算采用雙時(shí)間步方法進(jìn)行物理時(shí)間推進(jìn)。非定常計(jì)算的動(dòng)網(wǎng)格生成采用基于運(yùn)動(dòng)分解的耦合變形和重疊的動(dòng)態(tài)網(wǎng)格生成技術(shù)。為了適應(yīng)大規(guī)模工程計(jì)算的需求，早前已經(jīng)發(fā)展了基于網(wǎng)格分區(qū)的非阻塞MPI通信的大規(guī)模并行計(jì)算技術(shù)，并且在中等規(guī)模時(shí)達(dá)到了較高的MPI并行效率。本文僅對(duì)HyperFLOW中的非結(jié)構(gòu)解算器進(jìn)行MPI+OpenMP混合并行改造。

2 MPI+OpenMP混合并行計(jì)算方法

HyperFLOW軟件的并行通信數(shù)據(jù)結(jié)構(gòu)是基于網(wǎng)格分區(qū)構(gòu)建的，1個(gè)進(jìn)程可以處理1個(gè)或者串行處理幾個(gè)網(wǎng)格分區(qū)，進(jìn)程之間采用MPI消息傳遞機(jī)制進(jìn)行數(shù)據(jù)傳遞，如圖1(a)所示。每個(gè)分區(qū)內(nèi)在進(jìn)行迭代計(jì)算時(shí)，采用基于面(或單元)的數(shù)據(jù)結(jié)構(gòu)。

圖1 幾種并行模式的比較Fig.1 Comparison of different parallelization modes

這種基于網(wǎng)格分區(qū)的MPI通信模式，其通信過(guò)程可以完全人為指定，適合具有分布式內(nèi)存的計(jì)算機(jī)系統(tǒng)，能夠有效地提高并行規(guī)模。但是，隨著分區(qū)數(shù)的增大，分區(qū)交界面單元數(shù)量、通信量、通信次數(shù)呈幾何級(jí)數(shù)增大，導(dǎo)致通信時(shí)間占計(jì)算總時(shí)間也明顯增大，從而使得并行效率存在瓶頸。

為了減少通信量，提高并行效率，可以利用共享內(nèi)存的特點(diǎn)減少通信次數(shù)和通信量。結(jié)合HyperFLOW軟件的數(shù)據(jù)結(jié)構(gòu)特點(diǎn)，主要可以采用兩種共享內(nèi)存方式。

1) 粗粒度OpenMP模式,如圖1(b)所示。由于進(jìn)程之間需要通信，導(dǎo)致效率下降，因而可以保持分區(qū)數(shù)不變，減少進(jìn)程數(shù)，同時(shí)引入線程級(jí)并行，將一部分原來(lái)的進(jìn)程用線程替代，進(jìn)程間采用MPI通信，但是線程間采用共享內(nèi)存，從而減少通信量。由于線程對(duì)網(wǎng)格分區(qū)進(jìn)行并行處理，數(shù)據(jù)粒度大，因而也稱(chēng)這種并行模式為“粗粒度”O(jiān)penMP并行模式，這與文獻(xiàn)[13]提及的基于兩級(jí)分區(qū)的粗粒度混合并行模式是一致的。

2) 細(xì)粒度OpenMP模式，如圖1(c)所示。除上述模式外，也可以減少分區(qū)數(shù)，同時(shí)減少進(jìn)程數(shù)，進(jìn)程間仍然采用MPI通信，但是在網(wǎng)格分區(qū)內(nèi)部進(jìn)行線程級(jí)并行，達(dá)到減少通信量的目的。具體來(lái)說(shuō)，在基于面(或單元)的數(shù)據(jù)結(jié)構(gòu)基礎(chǔ)上引入線程級(jí)并行，也即線程對(duì)基于面(或單元)的數(shù)據(jù)進(jìn)行并行處理，數(shù)據(jù)粒度小，因而這種并行模式也被稱(chēng)為“細(xì)粒度”O(jiān)penMP并行模式。

2.1 純MPI并行

如圖1(a)所示，對(duì)于純MPI并行模式，1個(gè)節(jié)點(diǎn)內(nèi)部包含多個(gè)CPU，每個(gè)CPU上運(yùn)行1個(gè)進(jìn)程，處理1個(gè)或幾個(gè)網(wǎng)格分區(qū)。進(jìn)程間獨(dú)立計(jì)算，通過(guò)MPI通信分區(qū)界面數(shù)據(jù)。HyperFLOW中采用對(duì)所有網(wǎng)格分區(qū)進(jìn)行遍歷的模式進(jìn)行通信，具體分為以下兩個(gè)步驟：

1) 通信數(shù)據(jù)準(zhǔn)備。遍歷每個(gè)網(wǎng)格分區(qū)，計(jì)算出該分區(qū)需要發(fā)送的數(shù)據(jù)內(nèi)容和長(zhǎng)度。

2) 非阻塞式數(shù)據(jù)發(fā)送和接收。每個(gè)進(jìn)程均遍歷所有網(wǎng)格塊，當(dāng)所遍歷的網(wǎng)格塊屬于當(dāng)前進(jìn)程時(shí)，則向鄰居網(wǎng)格塊所在進(jìn)程發(fā)送數(shù)據(jù)；當(dāng)所遍歷的網(wǎng)格塊屬于其鄰居網(wǎng)格塊時(shí)，則從鄰居網(wǎng)格塊所在進(jìn)程接收數(shù)據(jù)。

圖2是MPI通信模式示意，圖中給出了4個(gè)進(jìn)程、4個(gè)網(wǎng)格塊，將每個(gè)網(wǎng)格塊分配至具有相同編號(hào)的進(jìn)程。以1號(hào)進(jìn)程為例，遍歷4個(gè)網(wǎng)格塊，該進(jìn)程上只有1號(hào)網(wǎng)格塊(與3號(hào)進(jìn)程中的3號(hào)網(wǎng)格塊為鄰居關(guān)系)，當(dāng)遍歷到1號(hào)網(wǎng)格塊時(shí)向3號(hào)進(jìn)程發(fā)送消息，當(dāng)遍歷到3號(hào)網(wǎng)格塊時(shí)從1號(hào)網(wǎng)格塊接收消息，而對(duì)于0號(hào)、2號(hào)網(wǎng)格塊則跳過(guò)[13]。

圖2 MPI通信模式[13]Fig.2 MPI parallelization mode[13]

以下采用CRM標(biāo)模定常湍流算例[10]對(duì)MPI效率進(jìn)行測(cè)試，標(biāo)模網(wǎng)格為非結(jié)構(gòu)混合網(wǎng)格，網(wǎng)格量約為4 000萬(wàn)單元，具體網(wǎng)格量如表1所示。

表1 CRM標(biāo)模測(cè)試算例網(wǎng)格信息Table 1 Grid information of CRM test case

定常湍流數(shù)值模擬采用二階精度有限體積法對(duì)RANS方程進(jìn)行離散求解，無(wú)黏通量格式選擇Roe格式，梯度重構(gòu)為GG-Cell方法，黏性通量選擇法向?qū)?shù)法[14]，湍流模型采用SA一方程模型，時(shí)間推進(jìn)采用隱式LU-SGS方法。具體方法參見(jiàn)文獻(xiàn)[3]，這里不再贅述。計(jì)算馬赫數(shù)Ma=0.85，雷諾數(shù)Re=5×106，攻角α=2.5°，下文中CRM算例測(cè)試條件均與本節(jié)相同。

測(cè)試集群為中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所的國(guó)產(chǎn)in-house集群(簡(jiǎn)稱(chēng)國(guó)產(chǎn)in-house集群)，CPU采用Phytium FT1500A芯片，單節(jié)點(diǎn)16核，共享32 GB內(nèi)存，采用國(guó)產(chǎn)高速互連網(wǎng)絡(luò)，雙向鏈路速率為200 Gbps，運(yùn)行國(guó)產(chǎn)麒麟操作系統(tǒng)，C++編譯器版本為Phytium 1.0.0(Based on GNU GCC 4.9.3)，MPI庫(kù)版本為MPICH version 3.2，OpenMP庫(kù)版本為OpenMP API 3.1。

圖3給出了加速比和并行效率的測(cè)試結(jié)果。最小測(cè)試規(guī)模為64核，最大測(cè)試規(guī)模為8 192核并行，不同并行規(guī)模的網(wǎng)格均采用8 192分區(qū)。結(jié)果顯示在1 024核并行時(shí)，相對(duì)于64核的MPI并行效率為99.8%，加速比為15.97，接近理想加速比。但是在并行規(guī)模進(jìn)一步增大時(shí)，并行效率急劇下降，當(dāng)并行核數(shù)為8 192核時(shí)，程序的并行效率只有37.9%，加速比僅達(dá)48左右，與理想加速比128存在較大差距。這是因?yàn)殡S著并行規(guī)模增大，單核處理的網(wǎng)格量減少，在8 192核時(shí)，單核處理的物理網(wǎng)格量只有不到5 000個(gè)單元，此時(shí)單核的計(jì)算量很小，而通信量隨著核數(shù)增加而急劇增大，從而使得并行效率嚴(yán)重下降。這體現(xiàn)出MPI并行模式在超大規(guī)模并行計(jì)算時(shí)存在的效率瓶頸問(wèn)題，必須通過(guò)減少通信時(shí)間占比來(lái)提高并行效率。利用多核處理器節(jié)點(diǎn)內(nèi)共享內(nèi)存的特性，將程序改造成節(jié)點(diǎn)間采用MPI通信、節(jié)點(diǎn)內(nèi)采用OpenMP共享內(nèi)存的兩級(jí)混合并行模式是一種減少通信量的可行辦法。

圖3 純MPI模式的并行效率及加速比Fig.3 Parallel efficiency and speedup of MPI parallelization mode

2.2 混合并行改造與實(shí)例

粗粒度混合并行是對(duì)網(wǎng)格分區(qū)進(jìn)行線程級(jí)并行，因此在不改變MPI通信模式的情況下，將對(duì)網(wǎng)格分區(qū)(zone)的循環(huán)進(jìn)行OpenMP并行改造，即可實(shí)現(xiàn)粗粒度混合并行，如圖4所示。

細(xì)粒度混合并行是在網(wǎng)格分區(qū)內(nèi)部，在基于面(或單元)的數(shù)據(jù)結(jié)構(gòu)上進(jìn)行線程級(jí)并行。因此，在不改變MPI通信模式的情況下，將基于面(或單元)的循環(huán)進(jìn)行OpenMP并行改造，即可實(shí)現(xiàn)細(xì)粒度混合并行，如圖5所示。

圖5 細(xì)粒度混合并行改造原理Fig.5 Principle of fine-grain hybrid parallelization

改造中需要注意的兩點(diǎn)：① 避免不同線程之間的數(shù)據(jù)競(jìng)爭(zhēng)，主要是數(shù)據(jù)寫(xiě)競(jìng)爭(zhēng)，如多個(gè)線程同時(shí)對(duì)同一地址變量進(jìn)行賦值操作；② 對(duì)計(jì)算耗時(shí)較大的熱點(diǎn)函數(shù)和模塊逐個(gè)進(jìn)行改造，如無(wú)黏通量、黏性通量、時(shí)間步長(zhǎng)計(jì)算、時(shí)間推進(jìn)等模塊。

借助性能分析工具，可以分析得到程序中的熱點(diǎn)類(lèi)/函數(shù)如表2所示，將類(lèi)和模塊對(duì)應(yīng)到各個(gè)程序功能。結(jié)果顯示，最耗時(shí)的類(lèi)是層流方程及湍流方程的LU-SGS推進(jìn)，其次依次是黏性通量計(jì)算、梯度計(jì)算、無(wú)黏通量、時(shí)間步長(zhǎng)及流場(chǎng)更新、湍流源項(xiàng)計(jì)算等。

表2 程序熱點(diǎn)函數(shù)及功能模塊對(duì)照表Table 2 Table of hot spots class/function and modules

如圖6所示，在非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)結(jié)構(gòu)中，由于不同面可能對(duì)應(yīng)相同的左(右)單元，因而可能會(huì)出現(xiàn)不同線程對(duì)同一個(gè)左(右)單元的右端項(xiàng)值rhsField進(jìn)行更新的情況，如圖6中face1和face2的左單元均為cell1，若face1和face2被分別分配到不同線程進(jìn)行并行計(jì)算，就存在對(duì)cell1的右端項(xiàng)值進(jìn)行同時(shí)更新的可能性。此時(shí)就必須采用atomic原子更新來(lái)避免不同線程對(duì)同一個(gè)地址進(jìn)行累加/減操作，避免造成數(shù)據(jù)沖突，導(dǎo)致不可預(yù)測(cè)的錯(cuò)誤，以基于C++語(yǔ)言的HyperFLOW程序?yàn)槔a如圖7所示。

圖6 非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)結(jié)構(gòu)示意圖Fig.6 Sketch of unstructured mesh data structure

圖7 避免數(shù)據(jù)競(jìng)爭(zhēng)的原子更新操作Fig.7 Avoidance of data races via atomic manipulation

如果不采用原子更新，該部分程序或其他類(lèi)似情況就無(wú)法實(shí)現(xiàn)OpenMP并行，會(huì)導(dǎo)致程序的并行度降低，但同時(shí)也應(yīng)指出atomic原子更新是在并行程序中強(qiáng)制只允許某單一線程對(duì)數(shù)據(jù)進(jìn)行更新，在保證結(jié)果穩(wěn)定正確的同時(shí)也會(huì)帶來(lái)額外的并行開(kāi)銷(xiāo)，是否采用原子更新需要對(duì)程序并行度和并行開(kāi)銷(xiāo)進(jìn)行權(quán)衡。

以下采用CRM標(biāo)模算例考察原子更新對(duì)并行效率和運(yùn)行時(shí)間的影響，分別采用4個(gè)線程不同進(jìn)程(16進(jìn)程依次增加到2 048進(jìn)程，核數(shù)從64核依次增大到8 192核)進(jìn)行并行測(cè)試，結(jié)果如圖8 所示，atomic原子更新可以提高程序并行度，提高大規(guī)模并行時(shí)的并行效率，但是也會(huì)導(dǎo)致程序在小規(guī)模并行時(shí)存在額外的并行開(kāi)銷(xiāo)，耗時(shí)增大。如圖8所示，圖8(a)的并行效率比較顯示atomic 原子更新效率相對(duì)更高。圖8(b)給出了采用atomic更新和不采用atomic更新的墻上時(shí)間差值，差值大于0，采用atomic更新耗時(shí)更長(zhǎng)，反之，采用atomic更新耗時(shí)更短。顯然，atomic 原子更新在大規(guī)模時(shí)效率高，耗時(shí)短，優(yōu)勢(shì)更加明顯，本文測(cè)試規(guī)模達(dá)到萬(wàn)核級(jí)別，因此測(cè)試中均采用原子更新提高程序并行度，減少程序耗時(shí)。

圖8 原子更新對(duì)并行效率和程序耗時(shí)的影響Fig.8 Effects of atomic manipulation on parallel efficiency and elapsed time

除采用atomic原子更新避免數(shù)據(jù)競(jìng)爭(zhēng)之外，由于C++中采用大量基于類(lèi)(Class)的數(shù)據(jù)結(jié)構(gòu)，如InviscidFluxSolverClass為計(jì)算無(wú)黏通量的類(lèi)，直接在類(lèi)對(duì)象中完成計(jì)算，可以大量減少以形參和實(shí)參形式的數(shù)據(jù)傳值，具有較好的封裝性和編程效率，但是在進(jìn)行線程級(jí)并行時(shí)，類(lèi)中的成員數(shù)據(jù)需要進(jìn)行并行化以避免數(shù)據(jù)競(jìng)爭(zhēng)。HyperFLOW中采用直接對(duì)類(lèi)對(duì)象進(jìn)行線程并行化，如圖9所示，每個(gè)線程中都有一個(gè)計(jì)算無(wú)黏通量的inviscidFluxSolverObject對(duì)象，線程間互相獨(dú)立，各自進(jìn)行對(duì)象初始化和無(wú)黏通量計(jì)算，徹底避免數(shù)據(jù)競(jìng)爭(zhēng)的問(wèn)題。

圖9 避免數(shù)據(jù)競(jìng)爭(zhēng)的對(duì)象并行Fig.9 Avoidance of data races via object parallelization

此外，程序中還引入了C++智能指針auto_ptr和shared_ptr，利用智能指針自動(dòng)刪除的特性，避免多個(gè)線程在刪除指針時(shí)，出現(xiàn)內(nèi)存錯(cuò)誤，如圖10所示。

圖10 智能指針與傳統(tǒng)指針比較Fig.10 Comparison of shared pointer and traditional pointer

同時(shí)，要注意改造前后計(jì)算結(jié)果應(yīng)該保持一致。具體來(lái)說(shuō)，改造后單線程混合并行的計(jì)算結(jié)果(包括收斂歷程和最終結(jié)果)應(yīng)當(dāng)與改造前完全一致，也即混合并行在單線程時(shí)可以退化為純MPI并行。而在多線程時(shí)，由于非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)存儲(chǔ)無(wú)序，導(dǎo)致自動(dòng)分配到各個(gè)線程間的數(shù)據(jù)依賴(lài)關(guān)系不明確，在隱式LUSGS改造時(shí)，難以實(shí)現(xiàn)像結(jié)構(gòu)網(wǎng)格那樣具有明確依賴(lài)關(guān)系的流水線式(即前一線程執(zhí)行完、數(shù)據(jù)更新后，再通知依賴(lài)前一線程數(shù)據(jù)的下一線程繼續(xù)執(zhí)行)的并行求解[15]。

由于粗粒度混合并行是在網(wǎng)格分區(qū)上進(jìn)行的并行，各個(gè)分區(qū)數(shù)據(jù)相互獨(dú)立，發(fā)生數(shù)據(jù)競(jìng)爭(zhēng)的可能性很小，只需要將全局變量可能存在的數(shù)據(jù)競(jìng)爭(zhēng)消除，因此代碼的改動(dòng)量很少。相反，細(xì)粒度混合并行要對(duì)多個(gè)熱點(diǎn)函數(shù)/類(lèi)進(jìn)行改造來(lái)實(shí)現(xiàn)線程級(jí)并行，同時(shí)避免數(shù)據(jù)競(jìng)爭(zhēng)，代碼改動(dòng)量和難度是粗粒度混合并行的數(shù)倍。

3 CRM標(biāo)模并行效率測(cè)試

3.1 OpenMP效率

在2.1節(jié)中，采用CRM標(biāo)模算例進(jìn)行了純MPI的效率測(cè)試，本節(jié)將采用CRM標(biāo)模算例分別對(duì)粗粒度和細(xì)粒度混合并行方式進(jìn)行測(cè)試。

將測(cè)試網(wǎng)格分別劃分為1 024和8 192個(gè)網(wǎng)格分區(qū)，分別采用64進(jìn)程和512進(jìn)程進(jìn)行粗粒度和細(xì)粒度混合并行計(jì)算，最小線程數(shù)為1線程，最大線程數(shù)為16線程。同時(shí)，對(duì)于細(xì)粒度混合并行，也可將測(cè)試網(wǎng)格分別劃分為64和512個(gè)分區(qū)，采用對(duì)應(yīng)進(jìn)程數(shù)(64進(jìn)程和512進(jìn)程，即分區(qū)數(shù)與進(jìn)程數(shù)對(duì)應(yīng))進(jìn)行并行效率和計(jì)算耗時(shí)測(cè)試，測(cè)試算例設(shè)置如表3所示。

表3 OpenMP并行效率測(cè)試算例Table 3 Test cases for OpenMP parallel efficiency

測(cè)試結(jié)果如圖11所示，圖中p代表進(jìn)程，z代表分區(qū)，結(jié)果顯示：

圖11 混合并行OpenMP效率測(cè)試結(jié)果Fig.11 Hybrid OpenMP parallel efficiency test results

1) 隨著線程數(shù)增大，幾種并行方式的OpenMP并行效率均下降，且分區(qū)更多時(shí)效率隨著線程數(shù)增加下降更快。

以512進(jìn)程為例，在采用不同線程時(shí)，粗粒度混合并行(圖11(a)中黑色虛線)各線程的OpenMP并行效率分別為81%、56%、35%和17%，細(xì)粒度混合并行(大分區(qū)數(shù)，圖11(a)中紅色虛線)的并行效率分別為71%、40%、19%和14%；而細(xì)粒度混合并行(對(duì)應(yīng)分區(qū)數(shù)，圖11(a)中藍(lán)色虛線)的并行效率分別為86%、68%、45%和25%。可見(jiàn)，OpenMP的線程數(shù)增加會(huì)導(dǎo)致額外的并行開(kāi)銷(xiāo)，減少線程數(shù)有助于提高OpenMP并行效率。

2) 采用對(duì)應(yīng)分區(qū)數(shù)的細(xì)粒度混合并行可以大大減少大規(guī)模并行時(shí)網(wǎng)格分區(qū)數(shù)量，提高OpenMP并行效率。采用對(duì)應(yīng)分區(qū)數(shù)的細(xì)粒度混合并行(512進(jìn)程，512分區(qū)，圖11(a)中藍(lán)色虛線)效率明顯高于大分區(qū)數(shù)的粗粒度混合并行和細(xì)粒度混合并行(512進(jìn)程，8 192分區(qū)，圖11(a)中黑色和紅色虛線)。

3) 粗粒度混合并行在小規(guī)模時(shí)更有優(yōu)勢(shì)，但是隨著并行規(guī)模增大，并行效率下降較快，而對(duì)應(yīng)分區(qū)數(shù)的細(xì)粒度混合并行在大規(guī)模并行時(shí)更有優(yōu)勢(shì)。

比如，在進(jìn)程規(guī)模較小時(shí)(64進(jìn)程)，粗粒度16線程并行效率能夠達(dá)到55%，而細(xì)粒度(對(duì)應(yīng)分區(qū)數(shù))并行效率只能達(dá)到40%左右，此時(shí)粗粒度混合并行具有一定優(yōu)勢(shì)，但是在進(jìn)程規(guī)模較大時(shí)(512進(jìn)程)，粗粒度并行效率迅速下降到17%，而采用對(duì)應(yīng)分區(qū)數(shù)(512分區(qū))的細(xì)粒度并行效率也下降到了25%，但是比粗粒度混合并行效率高。同時(shí)，從圖11(b)中的結(jié)果也可看到，在64進(jìn)程時(shí)，粗粒度混合并行耗時(shí)最少，而在512進(jìn)程，16線程時(shí)，采用對(duì)應(yīng)分區(qū)數(shù)的細(xì)粒度混合并行耗時(shí)最少。

3.2 混合并行的MPI效率測(cè)試

由2.1節(jié)和3.1節(jié)結(jié)論可知，當(dāng)進(jìn)程數(shù)增大時(shí)，MPI效率下降，而當(dāng)線程數(shù)增大時(shí)，OpenMP效率也會(huì)降低。因此，對(duì)于某一固定的并行規(guī)模，線程數(shù)增加，OpenMP效率降低，而同時(shí)相應(yīng)進(jìn)程數(shù)減少，MPI效率提升，因而會(huì)存在一個(gè)最佳的線程數(shù)，此時(shí)MPI效率提升效果大于OpenMP效率下降幅度，整體能夠達(dá)到最佳的混合并行效率。

如圖12所示，分別給出了純MPI并行、不同線程數(shù)時(shí)粗粒度混合并行效率和細(xì)粒度混合并行效率，其中純MPI測(cè)試和粗粒度測(cè)試網(wǎng)格均分為8 192個(gè)分區(qū)，細(xì)粒度測(cè)試網(wǎng)格分區(qū)數(shù)對(duì)應(yīng)于最大進(jìn)程數(shù)。

圖12 混合并行MPI效率測(cè)試結(jié)果Fig.12 MPI hybrid parallel efficiency test results

由圖12(a)可見(jiàn)，當(dāng)線程數(shù)增大時(shí)，粗粒度混合并行效率先下降，在8線程時(shí)，并行效率出現(xiàn)小幅回升，并在16線程時(shí)最高，達(dá)到32.5%(8 192核相對(duì)64核)，原因是在線程數(shù)較少時(shí)，OpenMP效率隨著線程數(shù)增加下降幅度大，如圖11(a)所示，在引入2線程時(shí)，粗粒度OpenMP效率下降到80%，而MPI效率上升幅度小，如圖3所示，8 192 核并行效率從37.9%上升到4 096核的52.5%，導(dǎo)致混合并行效率整體下降，而隨著線程數(shù)增大、進(jìn)程數(shù)減少，MPI并行效率上升幅度變大，混合并行效率逐漸恢復(fù)。

圖12(b)結(jié)果顯示，當(dāng)線程數(shù)增大到8線程時(shí)，細(xì)粒度混合并行效率達(dá)到最高，相對(duì)于64核，8 192 核并行效率達(dá)到56.5%，高于純MPI并行效率和粗粒度混合并行效率。

4 大規(guī)模并行計(jì)算測(cè)試

采用雙時(shí)間步方法進(jìn)行非定常多體分離數(shù)值模擬通常對(duì)內(nèi)迭代的計(jì)算精度要求較高，這樣才能保證計(jì)算誤差的累積不會(huì)對(duì)計(jì)算造成致命后果。采用大規(guī)模網(wǎng)格進(jìn)行多體分離數(shù)值模擬是減少內(nèi)迭代誤差的手段之一。

針對(duì)機(jī)翼外掛物投放標(biāo)模算例，本文在機(jī)翼和掛彈原始非結(jié)構(gòu)網(wǎng)格(表4中original，4 500萬(wàn)網(wǎng)格單元)的基礎(chǔ)上，通過(guò)一次自適應(yīng)加密，生成了3.6億非結(jié)構(gòu)混合網(wǎng)格(表4中adapt1)，再次進(jìn)行自適應(yīng)加密，生成了28.8億非結(jié)構(gòu)混合網(wǎng)格(表4中adapt2)，自適應(yīng)加密的具體過(guò)程參見(jiàn)文獻(xiàn)[16]。分離過(guò)程動(dòng)態(tài)網(wǎng)格采用非結(jié)構(gòu)重疊網(wǎng)格技術(shù)進(jìn)行裝配，方法細(xì)節(jié)可參考文獻(xiàn)[17-19]，此處不再贅述。3套網(wǎng)格裝配結(jié)果如圖13所示。

圖13 機(jī)翼外掛物投放大規(guī)模測(cè)試網(wǎng)格Fig.13 Large scale grids for wing store separation tests

如表4所示，由于網(wǎng)格量巨大，網(wǎng)格文件達(dá)到幾十甚至幾百GB量級(jí)，為減少文件存儲(chǔ)和IO耗時(shí)，保證程序整體的運(yùn)行效率，采用基于分組的文件存儲(chǔ)模式和對(duì)等(Peer to Peer,P2P)模式的文件IO機(jī)制[13]。

將adapt1網(wǎng)格分別存儲(chǔ)在1 024個(gè)文件中，adapt2網(wǎng)格分別存儲(chǔ)在8 192個(gè)文件中，單個(gè)文件大小為15～30 MB左右，在國(guó)產(chǎn)in-house集群上，分別采用1 024進(jìn)程讀入網(wǎng)格3.6億網(wǎng)格，8 192進(jìn)程讀入28.8億網(wǎng)格，P2P模式耗時(shí)分別為13.72 s和37.28 s，達(dá)到了較高的IO效率。而相比之下，服務(wù)器(Client/Server mode,CS)模式讀入網(wǎng)格的時(shí)間要長(zhǎng)得多，如表4所示。

表4 機(jī)翼外掛物標(biāo)模網(wǎng)格大小及網(wǎng)格讀入耗時(shí)Table 4 Grid size of wing-store standard model and time spent on reading grid files

在國(guó)產(chǎn)in-house集群上分別采用不同進(jìn)程數(shù)對(duì)3套不同規(guī)模網(wǎng)格進(jìn)行重疊網(wǎng)格隱式裝配和非定常計(jì)算。結(jié)果如表5所示，對(duì)4 500萬(wàn)網(wǎng)格采用128進(jìn)程進(jìn)行重疊網(wǎng)格裝配和非定常計(jì)算，重疊網(wǎng)格裝配耗時(shí)42.65 s，內(nèi)迭代50步耗時(shí)757.5 s，重疊網(wǎng)格裝配耗時(shí)與非定常內(nèi)迭代計(jì)算耗時(shí)比值在5%左右。對(duì)3.6億網(wǎng)格采用1 536進(jìn)程×8線程進(jìn)行重疊網(wǎng)格裝配和非定常計(jì)算，重疊網(wǎng)格裝配耗時(shí)17.56 s，內(nèi)迭代50步耗時(shí)277.4 s，重疊網(wǎng)格裝配耗時(shí)與非定常內(nèi)迭代計(jì)算耗時(shí)比值也在5%左右。對(duì)28.8億網(wǎng)格采用12 288進(jìn)程進(jìn)行重疊網(wǎng)格裝配，耗時(shí)僅需52.52 s。

表5 機(jī)翼外掛物標(biāo)模網(wǎng)格裝配信息Table 5 Grid assemble information of wing store model

4.1 3.6億非結(jié)構(gòu)重疊網(wǎng)格萬(wàn)核并行測(cè)試

將3.6億adapt1網(wǎng)格分為12 288個(gè)分區(qū)(機(jī)翼網(wǎng)格8 192分區(qū)+外掛物網(wǎng)格4 096分區(qū))，分別在國(guó)產(chǎn)in-house集群和天河2號(hào)上采用8線程細(xì)粒度和8線程粗粒度混合并行進(jìn)行非定常多體分離測(cè)試，最小測(cè)試規(guī)模為384核(48進(jìn)程×8線程)，最大規(guī)模為12 288核(1 536進(jìn)程×8線程)。

湍流數(shù)值模擬仍采用前述二階精度有限體積離散的RANS方程求解，測(cè)試算法為雙時(shí)間步非定常方法，內(nèi)迭代步采用LUSGS隱式時(shí)間推進(jìn)，內(nèi)迭代步數(shù)取50步，無(wú)黏通量選擇Roe格式，黏性通量選擇法向?qū)?shù)法，湍流模型選擇SA模型，氣動(dòng)/運(yùn)動(dòng)采用松耦合方式進(jìn)行求解[20]。

圖14給出了在國(guó)產(chǎn)in-house集群(圖中標(biāo)記為CARDC)和天河2號(hào)(圖中標(biāo)記為T(mén)H2)上的并行效率和墻上時(shí)間測(cè)試結(jié)果。結(jié)果顯示，在CARDC的in-house集群上，細(xì)粒度12 288核并行效率能夠達(dá)到90%(以768核為基準(zhǔn))。在天河2號(hào)上，并行效率能保持在70%(以384核為基準(zhǔn))，而此時(shí)純MPI效率和粗粒度混合并行效率已經(jīng)下降至50%(以384核為基準(zhǔn))。這再次證實(shí)了MPI效率在大規(guī)模并行時(shí)存在瓶頸，同時(shí)在大規(guī)模并行時(shí)，粗粒度混合并行效率也低于細(xì)粒度混合并行效率。

圖14(b)給出了不同并行方式墻上時(shí)間的對(duì)比。結(jié)果顯示，在小規(guī)模時(shí)純MPI并行的墻上時(shí)間最小，而由于并行開(kāi)銷(xiāo)大，細(xì)粒度混合并行的墻上時(shí)間最大。隨著并行規(guī)模增大，由于純MPI效率逐漸低于混合并行，導(dǎo)致墻上時(shí)間的優(yōu)勢(shì)逐漸消失。在12 288核時(shí)，純MPI和混合并行墻上時(shí)間接近，如果進(jìn)一步增大并行規(guī)模，隨著混合并行效率優(yōu)勢(shì)擴(kuò)大，混合并行墻上時(shí)間有望小于純MPI。

圖14 機(jī)翼外掛物投放大規(guī)模并行測(cè)試Fig.14 Large-scale parallel tests for wing store separation

采用12 288核混合并行(1 536進(jìn)程×8線程)進(jìn)行多體分離計(jì)算，并與試驗(yàn)結(jié)果進(jìn)行對(duì)比，如圖15所示，圖中ux、uy、uz、dx、dy、dz、Cfx、Cfy、Cfz和ψ、θ、φ分別表示慣性系中掛彈在3個(gè)方向的速度、位移、氣動(dòng)力系數(shù)和歐拉角；ωx、ωy、ωz和Cmx、Cmy、Cmz分別表示掛彈在體軸系3個(gè)方向上的旋轉(zhuǎn)角速度和氣動(dòng)力矩系數(shù)，其中“EXP”表示試驗(yàn)結(jié)果，不帶“EXP”為本文數(shù)值模擬結(jié)果。結(jié)果顯示計(jì)算結(jié)果與試驗(yàn)結(jié)果[21]符合良好，驗(yàn)證了混合并行計(jì)算結(jié)果的正確性，圖16給出了分離過(guò)程中典型時(shí)刻的重疊網(wǎng)格。

圖15 混合并行計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.15 Comparison of hybrid parallel numerical and experimental results

圖16 機(jī)翼外掛物投放典型時(shí)刻重疊網(wǎng)格Fig.16 Overset grids during wing store separation process

4.2 28.8億非結(jié)構(gòu)重疊網(wǎng)格4.9萬(wàn)核并行測(cè)試

通過(guò)機(jī)翼外掛物28.8億超大規(guī)模網(wǎng)格進(jìn)行非定常多體分離效率測(cè)試，驗(yàn)證本文發(fā)展的混合并行方法的可擴(kuò)展性。

如4.1節(jié)中所述，將adapt2網(wǎng)格劃分為98 304個(gè) 分區(qū)，分別存儲(chǔ)在8 192個(gè)文件中，最少采用512進(jìn)程×8線程，即4 096核進(jìn)行測(cè)試。最大采用6 144進(jìn)程×8線程，即49 152核進(jìn)行測(cè)試。測(cè)試算法與4.1節(jié)相同。

圖17給出了在國(guó)產(chǎn)in-house集群上的并行效率和墻上時(shí)間測(cè)試結(jié)果。結(jié)果顯示在4.9 萬(wàn)核時(shí)，以4 096核為基準(zhǔn)并行效率達(dá)到55.3%，加速比達(dá)到6.6。

圖17 機(jī)翼外掛物投放超大規(guī)模并行測(cè)試Fig.17 Massive parallel test for wing store separation

表6給出了不同并行規(guī)模時(shí)，28.8億網(wǎng)格計(jì)算單元壁面距離、重疊網(wǎng)格裝配及內(nèi)迭代等幾個(gè)主要計(jì)算步驟的墻上時(shí)間。由表中數(shù)據(jù)可見(jiàn)，重疊網(wǎng)格裝配耗時(shí)與50步內(nèi)迭代耗時(shí)之比約為7%，在整個(gè)計(jì)算耗時(shí)中僅占很小一部分。而在并行規(guī)模較小時(shí)，由于單核網(wǎng)格量較大，如512進(jìn)程×8線程，即4 096核，此時(shí)單核網(wǎng)格量達(dá)到70.3萬(wàn)，單元最近壁面的查詢(xún)量較大，導(dǎo)致壁面距離計(jì)算耗時(shí)較大。

表6 超大規(guī)模網(wǎng)格并行測(cè)試各功能耗時(shí)Table 6 Time spent on main operations of massive scale parallel test

在采用4.9萬(wàn)核時(shí)，計(jì)算壁面距離耗時(shí)約3 200 s，重疊網(wǎng)格裝配時(shí)間約為59 s，內(nèi)迭代計(jì)算時(shí)間約為700 s，總的來(lái)說(shuō)單個(gè)外迭代步耗時(shí)約為1 h。因此，在該集群上，對(duì)于一個(gè)28億網(wǎng)格量級(jí)的多體分離非定常數(shù)值模擬任務(wù)，采用5萬(wàn) 核進(jìn)行數(shù)值模擬，預(yù)計(jì)在7～10天左右能夠得到200步外迭代的數(shù)值模擬結(jié)果。該結(jié)果比目前采用百核進(jìn)行千萬(wàn)量級(jí)網(wǎng)格多體分離數(shù)值模擬周期(大概2～3天)仍長(zhǎng)出不少，主要是因?yàn)樽钚”诿嬗?jì)算耗時(shí)太多，在后續(xù)工作中，需要專(zhuān)門(mén)針對(duì)壁面距離算法進(jìn)行優(yōu)化和改進(jìn)[21-22]。為進(jìn)一步提高加速比，縮短模擬時(shí)間，還可以引入GPU并行計(jì)算，比如基于MPI+OpenACC的CPU/GPU混合并行[23-24]等。

5 結(jié) 論

本文首先實(shí)現(xiàn)了兩種方式的MPI+OpenMP混合并行：粗粒度混合并行和細(xì)粒度混合并行。其次，通過(guò)CRM標(biāo)模(4 000萬(wàn)非結(jié)構(gòu)網(wǎng)格)算例，對(duì)OpenMP并行效率和混合并行效率進(jìn)行了測(cè)試。最后，為驗(yàn)證混合并行在超大規(guī)模非定常算例中的可擴(kuò)展性，分別采用3.6億和28.8億重疊網(wǎng)格進(jìn)行了多體分離非定常數(shù)值模擬和并行效率測(cè)試，得到如下主要結(jié)論：

1) 隨著進(jìn)程數(shù)增大，純MPI并行效率存在瓶頸。如在本文的測(cè)試中，并行效率由1 024核時(shí)的99%下降到8 192核時(shí)的38%。

2) 由于OpenMP并行效率隨著線程數(shù)增加而單調(diào)下降，因此對(duì)于固定的并行核數(shù)規(guī)模，存在一個(gè)最優(yōu)線程數(shù)。在本文的測(cè)試環(huán)境中，粗粒度混合并行最優(yōu)線程數(shù)為16線程，而細(xì)粒度混合并行最優(yōu)線程數(shù)為8。此時(shí)細(xì)粒度混合并行效率高于純MPI并行效率，具有一定的效率優(yōu)勢(shì)。

3) 分別采用1 536進(jìn)程和12 288進(jìn)程對(duì)3.6億網(wǎng)格和28.8億網(wǎng)格進(jìn)行文件并行讀入和重疊網(wǎng)格并行裝配，耗時(shí)均在幾十秒內(nèi)，對(duì)大規(guī)模非定常計(jì)算而言，耗時(shí)基本可以忽略，已經(jīng)達(dá)到了較高的效率。

4) 采用細(xì)粒度混合并行完成了對(duì)3.6億網(wǎng)格的非定常多體分離數(shù)值模擬，計(jì)算結(jié)果與試驗(yàn)結(jié)果符合較好，驗(yàn)證了混合并行計(jì)算的精度。

5) 對(duì)于3.6億網(wǎng)格，在國(guó)產(chǎn)in-house集群上，12 888核并行效率達(dá)到90%(以768核為基準(zhǔn))，在天河二號(hào)上，并行效率達(dá)到70%(以384核為基準(zhǔn))；對(duì)于28.8億網(wǎng)格，在國(guó)產(chǎn)in-house集群上，4.9萬(wàn)核并行效率達(dá)到55.3%(以4 096核為基準(zhǔn))，通過(guò)上述算例驗(yàn)證了本文混合并行改造方法的可行性。

當(dāng)然，需要特別指出的是，目前僅在國(guó)產(chǎn)in-house集群和“天河二號(hào)”上進(jìn)行了測(cè)試。這兩臺(tái)機(jī)器均屬于“胖節(jié)點(diǎn)”型的集群，即每個(gè)計(jì)算節(jié)點(diǎn)內(nèi)有若干CPU，每個(gè)CPU內(nèi)有若干核，而計(jì)算節(jié)點(diǎn)間利用高速互聯(lián)網(wǎng)絡(luò)相連。對(duì)于其他架構(gòu)的并行計(jì)算機(jī)，本文的結(jié)論是否同樣適用尚待進(jìn)一步驗(yàn)證。

未來(lái)將進(jìn)行更大規(guī)模并行測(cè)試，同時(shí)根據(jù)計(jì)算機(jī)體系結(jié)構(gòu)的特點(diǎn)，開(kāi)展大規(guī)模異構(gòu)并行計(jì)算方法研究，比如開(kāi)展MPI+OpenMP+CUDA或者M(jìn)PI+OpenACC的混合并行研究，以使程序適應(yīng)基于CPU+GPU等環(huán)境的異構(gòu)體系計(jì)算機(jī)。此外，將進(jìn)一步改進(jìn)隱式方法的細(xì)粒度混合并行，提高隱式格式收斂效率；改進(jìn)壁面距離計(jì)算方法，減少壁面距離計(jì)算耗時(shí)。

致 謝

感謝中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所鄧亮博士的指導(dǎo)。