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

王年華，常興華，趙鐘，張來平,

1.中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室，綿陽 621000 2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

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

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

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

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

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

1 HyperFLOW軟件簡介

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

2 MPI+OpenMP混合并行計算方法

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

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

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

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

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

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

2.1 純MPI并行

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

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

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

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

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

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

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

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

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

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

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

2.2 混合并行改造與實例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 CRM標模并行效率測試

3.1 OpenMP效率

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

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

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

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

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

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

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

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

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

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

3.2 混合并行的MPI效率測試

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

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

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

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

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

4 大規(guī)模并行計算測試

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

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

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

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

致 謝

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