徐傳福,車永剛,李大力,王勇獻(xiàn),王正華

(1.國防科技大學(xué)高性能計算國家重點實驗室,湖南 長沙 410073;2.國防科技大學(xué)計算機(jī)學(xué)院,湖南 長沙 410073;3.國防科技大學(xué)氣象海洋學(xué)院,湖南 長沙 410073)

1 引言

計算流體力學(xué)CFD(Computational Fluid Dynamics)通過數(shù)值求解各種流體動力學(xué)控制方程，獲取各種條件下的流動數(shù)據(jù)和作用在繞流物體上的力、力矩和熱量等，從而達(dá)到研究各種流動現(xiàn)象和規(guī)律的目的。CFD是涉及流體力學(xué)、計算機(jī)科學(xué)與技術(shù)、計算數(shù)學(xué)等多個專業(yè)的交叉研究領(lǐng)域。隨著高性能計算機(jī)的飛速發(fā)展，CFD研究和工程實踐都取得了很大進(jìn)步，20世紀(jì)90年代以來，基于雷諾平均Navier-Stokes方程求解的CFD技術(shù)已經(jīng)廣泛應(yīng)用于航空、航天、航海、能源動力、環(huán)境、機(jī)械裝備等諸多國民經(jīng)濟(jì)和國防安全領(lǐng)域，取得了很好的應(yīng)用效果。在航空航天領(lǐng)域，CFD已逐漸成為與理論分析、風(fēng)洞實驗并列的流體力學(xué)3大主要方法之一。美國國家航天局(NASA)預(yù)測，21世紀(jì)，高效能計算機(jī)和CFD技術(shù)的進(jìn)一步結(jié)合將給各類航空航天飛行器的氣動設(shè)計帶來一場革命[1]。

高性能計算技術(shù)的迅猛發(fā)展為大規(guī)模復(fù)雜CFD應(yīng)用提供了重要支撐，當(dāng)前，CFD已經(jīng)成為高性能計算機(jī)上最重要的應(yīng)用之一。隨著應(yīng)用問題復(fù)雜度的增加，CFD要求的幾何外形、數(shù)值方法、物理化學(xué)模型等也日益復(fù)雜、精細(xì)，對大規(guī)模計算提出了更高要求。CFD對大規(guī)模計算的需求可以從能力計算(Capability Computing)和容量計算(Capacity Computing)2方面概括。容量計算指的是利用超級計算機(jī)同時完成大批量生產(chǎn)性業(yè)務(wù)，在CFD中通常用于復(fù)雜工程問題的設(shè)計與優(yōu)化，例如飛行器全包線數(shù)據(jù)庫生產(chǎn)等。有學(xué)者在1997年估計，商業(yè)飛機(jī)巡航一秒鐘的計算，用每秒萬億次計算機(jī)需要數(shù)千年，高保真度全包線氣動數(shù)據(jù)庫的生產(chǎn)被認(rèn)為是CFD一個長期的重大挑戰(zhàn)問題[2,3]。能力計算通常指的是利用超級計算機(jī)全系統(tǒng)計算能力求解單個大型任務(wù)，在CFD中通常用于簡單外形、復(fù)雜流動問題的基礎(chǔ)研究，例如采用直接數(shù)值模擬開展湍流流動機(jī)理研究等。據(jù)美國波音公司Tinoco博士2009年估計，以當(dāng)時高性能計算機(jī)的發(fā)展速度，直到2080年左右才可能進(jìn)行民航客機(jī)全機(jī)的DNS模擬；即便是進(jìn)行大渦模擬也要等到2045年左右[4]。

CFD巨大的計算量對于超級計算機(jī)研制和超大規(guī)模并行計算研究提出了迫切需求，異構(gòu)并行架構(gòu)是當(dāng)前構(gòu)建超大規(guī)模高性能計算機(jī)系統(tǒng)的重要技術(shù)途徑之一[5]。異構(gòu)超級計算機(jī)主要包括異構(gòu)加速器和異構(gòu)眾核2種實現(xiàn)方式。例如，我國的天河系列超級計算機(jī)采用了異構(gòu)加速器，其中天河一號的加速器為通用GPU(Graphics Processing Unit)，而天河二號的加速器為Intel集成眾核MIC(Many Integrated Cores)(升級后的天河二號采用了國產(chǎn)加速器Matrix2000)；神威太湖之光則采用了“申威26010”異構(gòu)眾核處理器，每個處理器包括4個主計算核，每個主計算核配有一個8×8的計算陣列(64個從計算核)。在2020年6月發(fā)布的世界超級計算機(jī)排行榜(TOP500)中，排名前10的超級計算機(jī)有8臺是異構(gòu)超級計算機(jī)。異構(gòu)超級計算機(jī)具有明顯的性能價格比、性能功耗比等優(yōu)勢，但異構(gòu)超級計算機(jī)具有異構(gòu)的計算、存儲和通信能力以及編程環(huán)境，極大增加了高效、大規(guī)模CFD應(yīng)用軟件開發(fā)的難度。在CFD應(yīng)用領(lǐng)域，傳統(tǒng)CPU平臺并行計算主要采用分區(qū)并行方法，每個分區(qū)獨立進(jìn)行求解，利用消息傳遞通信實現(xiàn)分區(qū)之間的任務(wù)并行以及共享存儲實現(xiàn)單個分區(qū)內(nèi)部的線程并行[6 - 9]。異構(gòu)超級計算機(jī)具有多層次、多粒度的異構(gòu)并行性，應(yīng)用并行性開發(fā)需要同時利用消息傳遞任務(wù)級并行、計算結(jié)點內(nèi)CPU與加速器之間的協(xié)同并行、CPU/加速器上的共享存儲線程級并行和CPU/加速器上的向量化指令級并行，需要針對異構(gòu)并行體系結(jié)構(gòu)特征，設(shè)計多層次可擴(kuò)展并行算法，才能實現(xiàn)CFD應(yīng)用、算法與并行體系結(jié)構(gòu)的“最佳適配”，充分挖掘超高性能計算機(jī)潛力。

國防科技大學(xué)不僅是我國高性能計算機(jī)系統(tǒng)研制的基地，也是我國高性能計算應(yīng)用軟件研發(fā)的基地。長期以來，國防科技大學(xué)CFD應(yīng)用軟件團(tuán)隊依托天河/銀河系列超級計算機(jī)開展了超大規(guī)模復(fù)雜CFD并行計算和性能優(yōu)化研究，突破了異構(gòu)協(xié)同并行計算等關(guān)鍵技術(shù)，實現(xiàn)了HPC與CFD的深度融合，有力支撐了我國幾套重要的In-house CFD應(yīng)用軟件在天河/銀河系列超級計算機(jī)上的大規(guī)模并行計算。本文歸納總結(jié)了作者團(tuán)隊基于自主高精度CFD軟件，面向航空航天氣動數(shù)值模擬，在天河超級計算機(jī)上開展的超大規(guī)模高精度CFD并行計算研究，并對未來E級超級計算機(jī)上CFD并行應(yīng)用開發(fā)進(jìn)行了分析展望。

2 研究現(xiàn)狀

近年來，隨著P級超級計算機(jī)的研制成功，歐美日等發(fā)達(dá)國家在這些最高端的計算平臺上針對湍流等復(fù)雜流動機(jī)理研究開展了超大規(guī)模CFD并行計算。例如，2013年，德克薩斯州州立大學(xué)研究人員實現(xiàn)了P級高雷諾數(shù)槽道流的直接數(shù)值模擬，并行規(guī)模達(dá)到約78萬處理器核[10]。瑞士蘇黎世聯(lián)邦工學(xué)院、IBM蘇黎世實驗室等單位聯(lián)合完成了基于有限體積法的無粘可壓兩相流模擬，最大網(wǎng)格規(guī)模達(dá)13萬億網(wǎng)格點，獲得了11PFLOPS的持續(xù)性能，達(dá)到系統(tǒng)峰值性能的55%，該項工作獲得了2013年度戈登·貝爾獎[11]。

高性能異構(gòu)并行體系結(jié)構(gòu)發(fā)展至今，很多學(xué)者在GPU、MIC等異構(gòu)超級計算機(jī)上開展了大規(guī)模CFD異構(gòu)并行計算研究，取得了不錯的加速效果。GPU出現(xiàn)的早期，研究人員通常僅移植一些簡化的CFD代碼，以二階精度和一些簡單的流動問題模擬為主，計算平臺通常也僅包含1塊或幾塊GPU卡。通過早期實踐驗證GPU計算的加速效果后，大規(guī)模GPU異構(gòu)并行逐漸成為CFD并行計算研究的熱點。例如，Jacobsen等[12]實現(xiàn)了一個支持GPU集群的不可壓CFD求解器，在美國國家超級計算應(yīng)用中心的Tesla集群上利用64個結(jié)點(共128塊Tesla C1060 GPU卡)進(jìn)行了測試，相對于8 CPU核獲得約130倍的加速比。他們同時在256塊GPU上開展了頂蓋方腔管道湍流的大渦模擬[13]，采用了1維區(qū)域分解，對比了MPI+OpenMP+CUDA 并行和MPI+CUDA并行實現(xiàn)，但實際上他們的工作并未實現(xiàn)CPU-GPU協(xié)同，OpenMP僅用于代替結(jié)點內(nèi)MPI通信。作者團(tuán)隊在天河一號超級計算機(jī)上設(shè)計了基于MPI+CUDA+OpenMP的CPU-GPU異構(gòu)協(xié)同并行算法，成功實現(xiàn)了當(dāng)時世界上最大規(guī)模CPU-GPU協(xié)同并行高精度計算，模擬了三段翼構(gòu)型高精度氣動聲學(xué)問題和大型客機(jī)C919的高精度氣動力預(yù)測問題，問題規(guī)模達(dá)8億網(wǎng)格單元，并行規(guī)模達(dá)1 024個計算結(jié)點[14 - 17]。MIC架構(gòu)產(chǎn)品的出現(xiàn)晚于GPU，相關(guān)CFD問題應(yīng)用優(yōu)化與并行算法研究工作相對較少。德國宇航中心于2011年啟動了名為“面向眾核架構(gòu)的CFD代碼高效實現(xiàn)HICFD(Highly Efficient Implementation of CFD codes for HPC many-core architectures)”的研究項目[18]，面向眾核高性能計算機(jī)研究新的方法與工具，最優(yōu)地利用系統(tǒng)的全部并行層級，包括在最高層使用MPI，在眾核加速卡級使用高度可擴(kuò)展的MPI/OpenMP混合并行，在處理器核級高效利用SIMD部件。隨著天河二號超級計算機(jī)的發(fā)布，作者團(tuán)隊又開展了MIC平臺上CFD并行計算和性能優(yōu)化研究，設(shè)計了基于MPI+Offload+OpenMP+SIMD的CPU-MIC異構(gòu)協(xié)同并行算法，實現(xiàn)了數(shù)十億網(wǎng)格規(guī)模的可壓縮拐角直接數(shù)值模擬，并行規(guī)模擴(kuò)展到百萬異構(gòu)計算核心[19 - 24]。

整體而言，國內(nèi)多數(shù)CFD并行計算規(guī)模為數(shù)十到數(shù)百核，與國外相比仍然有較大差距。長期以來國內(nèi)多數(shù)CFD軟件更加關(guān)注CFD自身的專業(yè)性，與高性能計算機(jī)系統(tǒng)的研制相互脫節(jié)，CFD軟件適應(yīng)新型高性能并行體系結(jié)構(gòu)的能力較弱，迫切需要開展CFD算法、應(yīng)用和系統(tǒng)的深度融合研究，充分發(fā)揮新一代超級計算機(jī)系統(tǒng)的潛能。

3 CFD方法、軟件和計算平臺

3.1 數(shù)值方法和軟件

這里以一個In-house多區(qū)結(jié)構(gòu)網(wǎng)格高精度CFD軟件為例，簡要介紹高精度CFD數(shù)值方法和軟件實現(xiàn)。直角坐標(biāo)系下強(qiáng)守恒形式控制方程表示為：

上述方程通過坐標(biāo)變換(x,y,z,t)→(ξ,η,ζ,τ)轉(zhuǎn)換為曲線坐標(biāo)下方程：

該結(jié)構(gòu)網(wǎng)格高精度CFD軟件中實現(xiàn)了WCNS(Weighted Compact Non-linear Scheme)[25]、HDCS(Hybrid cell-edge and cell-node Dissipative Compact Scheme)[26]等我國自主發(fā)展的有限差分高精度計算格式，這里以5階顯式WCNS格式WCNS-E-5沿方向無粘通量導(dǎo)數(shù)離散為例，其內(nèi)點格式可以表示為：

高精度CFD軟件計算流程如圖1所示。迭代(定常時間步迭代或非定常子迭代)開始施加邊界條件，之后交換虛網(wǎng)格單元和奇異網(wǎng)格單元原始變量值，接著在計算和交換原始變量梯度之前計算譜半徑增量和時間步。WCNS、HDCS等高精度格式在右端項(粘性項和無粘項)計算中實現(xiàn)，右端項計算結(jié)果也需要進(jìn)行交換。軟件中實現(xiàn)了常見的LU-SGS、PR-SGS等隱式方法以及顯式Runge-Kutta方法，求得守恒變量增量后更新原始變量及殘差，循環(huán)結(jié)束。

Figure 1 Flowchart of the structured high-order CFD圖1 高精度結(jié)構(gòu)網(wǎng)格CFD軟件計算流程

3.2 異構(gòu)計算平臺

圖2給出了采用加速器的異構(gòu)計算平臺。圖2中每個計算結(jié)點包含若干共享內(nèi)存的多核CPU，計算結(jié)點間通過高速互連網(wǎng)絡(luò)進(jìn)行通信，每個計算結(jié)點包含若干加速器ACC(ACCelerator)，加速器通常具有片上存儲，通過PCI-e與CPU進(jìn)行數(shù)據(jù)交互。以天河一號為例，每個計算結(jié)點包含雙路Intel Xeon X5670 CPU和1個NVIDIA Tesla M2050 GPU；天河二號每個計算結(jié)點則包含雙路Intel Xeon E5-2692 v2 CPU和3塊MIC協(xié)處理器(Intel Xeon Phi 31S1P)。圖2同時給出了異構(gòu)計算平臺各層次對應(yīng)的編程模型。計算結(jié)點間通常采用消息傳遞接口MPI(Message Passing Interface)實現(xiàn)分布式并行，計算結(jié)點內(nèi)各CPU核上通常采用OpenMP實現(xiàn)共享存儲并行。不同加速器通常需要采用不同的編程模型，既有GPU專用的CUDA和MIC專用的Intel OffLoad編程模型，也有同時支持GPU、MIC等多種計算平臺的OpenACC、OpenMP4.X等編程模型。CFD應(yīng)用只有綜合利用上述并行編程模型，才能實現(xiàn)多層次并行算法。

Figure 2 Accelerator-based heterogeneous computing platform and its programming models圖2 加速器異構(gòu)計算平臺及其編程模型

4 超大規(guī)模異構(gòu)協(xié)同并行計算

4.1 異構(gòu)并行區(qū)域分解

區(qū)域分解是開展并行算法設(shè)計的第1步。圖3以三維多區(qū)網(wǎng)格CFD計算為例，給出了支持多層次異構(gòu)協(xié)同并行算法的區(qū)域分解示意圖。圖3中CFD流場區(qū)域首先根據(jù)負(fù)載均衡策略劃分為多個網(wǎng)格塊(為了滿足大規(guī)模并行計算及其負(fù)載均衡需求，通常需要對原始生成的單塊或多塊網(wǎng)格進(jìn)行二次剖分)，每個MPI進(jìn)程負(fù)責(zé)1個包含若干網(wǎng)格塊的分組，為了簡化編程，通常1個計算結(jié)點分配1個MPI進(jìn)程。結(jié)點內(nèi)網(wǎng)格塊分組需要根據(jù)CPU和加速器的計算、存儲能力在兩者之間進(jìn)行均衡分配，考慮到編程復(fù)雜度和PCI-e通信開銷等，通常以整個網(wǎng)格塊作為分配單位。對于CPU或MIC，每個網(wǎng)格塊內(nèi)沿著特定維度劃分為數(shù)據(jù)片(data chunk)分配給計算核實現(xiàn)OpenMP線程并行。MIC計算核較多，網(wǎng)格塊規(guī)模較小或所劃分的維度網(wǎng)格單元均較小時，可以采用嵌套OpenMP對網(wǎng)格塊的第2個維度進(jìn)行進(jìn)一步劃分。在CPU或MIC上，針對每個線程處理的數(shù)據(jù)片內(nèi)的網(wǎng)格線可以實現(xiàn)向量化并行。

GPU上的任務(wù)分配則較為復(fù)雜。圖3中給出了一種2層分配策略，由于當(dāng)前GPU均支持流處理，首先將網(wǎng)格塊分配給GPU流實現(xiàn)GPU上的任務(wù)級并行；進(jìn)一步，在每個網(wǎng)格塊內(nèi)，若CFD算法模型在計算中各網(wǎng)格單元沒有依賴關(guān)系，則可以設(shè)置三維的GPU線程塊，每個GPU線程處理1個網(wǎng)格單元。如果某一維度網(wǎng)格單元之間存在依賴關(guān)系，則可以考慮采用二維GPU線程塊。作者團(tuán)隊在天河一號上針對多區(qū)結(jié)構(gòu)網(wǎng)格高精度CFD軟件實現(xiàn)了這種2層策略，這種方法不僅可以充分挖掘GPU的多層次并行，同時可進(jìn)一步利用流處理對GPU計算的區(qū)塊進(jìn)行分組，克服了GPU存儲空間小對計算規(guī)模的限制，實現(xiàn)了CPU與GPU之間負(fù)載的靈活控制(CPU計算能力弱但存儲容量大，GPU計算能力強(qiáng)而存儲空間小)。

Figure 3 Domain decomposition for multi-level heterogeneous collaborative parallel computing圖3 多層次異構(gòu)協(xié)同并行區(qū)域分解示意圖

在區(qū)域分解過程中，網(wǎng)格剖分通常采用獨立的工具實現(xiàn)，CPU或加速器上的計算任務(wù)分配需要在CFD求解器代碼中實現(xiàn)。與同構(gòu)CPU平臺相比，異構(gòu)并行平臺對區(qū)域分解和負(fù)載均衡提出了更新的要求。例如，為了更好的負(fù)載均衡，異構(gòu)計算結(jié)點通常要求剖分的網(wǎng)格塊更多，考慮到加速器豐富的并行能力，分配給加速器的網(wǎng)格區(qū)塊又不宜太小。

4.2 異構(gòu)協(xié)同并行

異構(gòu)加速器極大地提升了異構(gòu)超級計算機(jī)整體性能。以天河一號為例，每個計算結(jié)點CPU雙精度浮點性能約140 GFLOPS，GPU雙精度浮點性能約500 GFLOPS，在異構(gòu)超級計算機(jī)上開展大規(guī)模CFD計算不僅需要用好CPU，更需要用好加速器，使得兩者能夠?qū)崿F(xiàn)高效協(xié)同并行。本文歸納總結(jié)了2種協(xié)同并行模式[14]：基于嵌套OpenMP的協(xié)同和基于ACC異步執(zhí)行的協(xié)同，如圖4所示。以每個計算結(jié)點配置1個加速器、網(wǎng)格包括NBLKCOMS個分塊為例：基于嵌套OpenMP的協(xié)同首先在CPU上創(chuàng)建第1層2個OpenMP線程，其中第1個線程控制編號為[1，ACC_BLOCK_NUM]的分塊在ACC上的計算，在第2個線程內(nèi)創(chuàng)建嵌套OpenMP線程，啟動其他CPU核計算編號為[ACC_BLOCK_NUM+1,NBLKCOMS]的分塊?；贏CC異步執(zhí)行的協(xié)同則利用了加速器異步編程模型和異步調(diào)度執(zhí)行機(jī)制，CPU主線程異步啟動ACC計算任務(wù)后，控制權(quán)立刻返回CPU主線程，此時啟動CPU上的多線程計算。在上述2種模式中CPU與ACC上都能夠同時運行計算任務(wù)，從而實現(xiàn)協(xié)同并行。CFD計算過程通常涉及邊界處理、邊界數(shù)據(jù)交換等操作，為了保證CPU上具有最新的計算結(jié)果，協(xié)同并行結(jié)束時需要在CPU與ACC之間進(jìn)行數(shù)據(jù)傳輸、同步等。

Figure 4 Two heterogeneous collaborative parallel modes圖4 2種異構(gòu)協(xié)同并行模式

上述2種協(xié)同并行模式編程實現(xiàn)均較為簡單，目前GPU、MIC等異構(gòu)加速器及CUDA、OpenMP4.X、OpenACC異構(gòu)編程模型都支持異步執(zhí)行和數(shù)據(jù)傳輸，開發(fā)人員可以在實際CFD程序中實現(xiàn)2種并行模式后針對不同的算例測試異構(gòu)協(xié)同并行效果。盡管理論上基于ACC異步執(zhí)行的協(xié)同并行似乎較為高效(無需專門留出CPU核控制ACC，所有CPU核均可參與計算)，但在作者過去的實踐中，基于嵌套OpenMP的協(xié)同并行效果更佳。

4.3 GPU并行

MIC上的OpenMP并行與CPU上的類似，這里重點介紹作者團(tuán)隊針對多塊網(wǎng)格CFD計算提出的2層GPU并行算法[14]：網(wǎng)格區(qū)塊內(nèi)基于CUDA的細(xì)粒度數(shù)據(jù)并行算法和網(wǎng)格區(qū)塊間基于CUDA流處理機(jī)制的粗粒度任務(wù)并行算法，如圖5所示。圖5左邊給出的多塊網(wǎng)格CFD計算過程包含nblk個網(wǎng)格分塊的塊循環(huán)(block-loop)，以及(K,J,I)三維大小為(NK,NJ,NI)網(wǎng)格分塊(虛邊界擴(kuò)充ngn個網(wǎng)格單元)內(nèi)的單元循環(huán)(cell-loop)；右邊給出了2層并行算法實現(xiàn)的偽代碼，其中對網(wǎng)格區(qū)塊的循環(huán)(block-loop)和對網(wǎng)格單元的循環(huán)(cell-loop)分別映射到CUDA中的流處理循環(huán)(stream-loop)和CUDA計算核函數(shù)(CUDA kernel)，這里假設(shè)GPU流的數(shù)量為num_stream。

Figure 5 Schematic diagram of two-level GPU parallel algorithm for multi-block structured CFD圖5 多區(qū)結(jié)構(gòu)網(wǎng)格CFD計算2層GPU并行算法示意圖

細(xì)粒度并行算法根據(jù)CFD不同計算過程所包括的循環(huán)內(nèi)部對網(wǎng)格單元的處理是否具有數(shù)據(jù)依賴關(guān)系而設(shè)計。對于網(wǎng)格單元間完全獨立的計算過程(例如Jacobi迭代、譜半徑計算等)，網(wǎng)格區(qū)塊采用三維分解，計算過程實現(xiàn)為三維的CUDA kernel，根據(jù)索引每個GPU線程可獨立計算1個網(wǎng)格單元。對于網(wǎng)格單元間存在依賴關(guān)系的計算過程(例如各方向的無粘、粘性通量計算)，由于CUDA不支持全局線程同步，因此可以采用二維分解，實現(xiàn)為二維CUDA kernel，每個GPU線程計算1條網(wǎng)格線。

在多區(qū)網(wǎng)格CFD中，任意2次流場邊界信息交換之間的不同網(wǎng)格分區(qū)之間的計算是獨立的，可以并行處理。CUDA通過流處理(Streaming)機(jī)制支持任務(wù)級并行，允許用戶將應(yīng)用問題分為多個相互獨立的任務(wù)，每個任務(wù)或者流定義了一個操作序列，同一流內(nèi)的操作需要滿足一定的順序，而不同流則可以在GPU上亂序執(zhí)行。粗粒度并行算法將每個GPU流綁定到一個網(wǎng)格區(qū)塊，同時在GPU上執(zhí)行多個GPU流，實現(xiàn)多個網(wǎng)格區(qū)塊的并發(fā)處理。流機(jī)制的引入一方面滿足了應(yīng)用問題多層次并行性開發(fā)的需求，另一方面很好地適應(yīng)了GPU的硬件資源特點，提高了資源利用效率。在多GPU流處理機(jī)制上進(jìn)一步設(shè)計了分組多流機(jī)制GBMS(Group-Based Multiple Streams)[14,15]，圖6a給出了多GPU流并發(fā)處理多個區(qū)塊時的狀態(tài)圖，可以看出多GPU流能夠重疊多個分區(qū)的拷入、計算和拷出，隱藏CPU和GPU間數(shù)據(jù)傳輸開銷；圖6b GBMS將4個流/分區(qū)分為2組，這種方式可有效克服天河一號Tesla M2050存儲容量小的限制，允許GPU計算更多網(wǎng)格區(qū)塊。

Figure 6 Schematic diagram of multi-stream parallel execution and GBMS for 4 grid blocks圖6 4個分區(qū)時多流并發(fā)執(zhí)行和分組多流并行示意圖

4.4 OpenMP并行

若計算過程中網(wǎng)格單元之間沒有依賴關(guān)系，則網(wǎng)格分區(qū)內(nèi)的OpenMP并行可以通過沿著網(wǎng)格單元循環(huán)(通常為最外層循環(huán))添加OpenMP編譯指導(dǎo)語句實現(xiàn)，較為簡單，本節(jié)重點介紹作者團(tuán)隊針對CFD中常用的具有強(qiáng)數(shù)據(jù)依賴關(guān)系的Gauss-Seidel迭代類算法(例如LU-SGS)改進(jìn)的共享存儲并行算法。以三維結(jié)構(gòu)網(wǎng)格為例(如圖7所示)，在LU-SGS算法向前(下三角矩陣)掃描時，網(wǎng)格點(i,j,k)的更新計算需要依賴小號鄰居點(i-1,j,k)、(i,j-1,k)和(i,j,k-1)的更新值，反之向后(上三角矩陣)掃描過程中則需要依賴大號鄰居點(i+1,j,k)、(i,j+1,k)和(i,j,k+1)的更新值。由于這一數(shù)據(jù)依賴特點，無法直接添加OpenMP指導(dǎo)語句實現(xiàn)LU-SGS的共享存儲OpenMP并行。

Figure 7 Data dependence in forward computing of LU-SGS algorithm 圖7 LU-SGS算法向前掃描時的數(shù)據(jù)依賴關(guān)系

為了實現(xiàn)LU-SGS的OpenMP并行計算，NASA等機(jī)構(gòu)的學(xué)者提出了對角-超平面和流水線2種并行算法[32,33]，主要思想是開發(fā)LU-SGS算法中各個不同網(wǎng)格面和網(wǎng)格線上計算沒有依賴關(guān)系的網(wǎng)格點進(jìn)行并行計算。在早期的多核處理器上，流水線LU-SGS并行算法取得了很好的并行加速比，然而我們的測試分析表明，在MIC新型眾核加速器上，228線程時并行效率急劇下降到25%以下，對于較小規(guī)模的算例甚至低至1%以下。其原因在于流水線并行效率受流水線建立和排空過程的限制。隨著流水線深度(線程數(shù))增加，流水線建立和排空的開銷也隨之增加。對于外層循環(huán)長度小于流水線深度的網(wǎng)格塊而言，甚至沒有足夠的計算任務(wù)充分填充所有的流水線級。另一方面，隨著線程數(shù)的增加，各個線程上的計算負(fù)載均衡性也會越來越差。

為了破解傳統(tǒng)流水線并行LU-SGS在MIC加速器上的并行擴(kuò)展性瓶頸，作者團(tuán)隊針對三維網(wǎng)格(3個維度記為K、J、I，KmJn指的是K、J維的第m、n個網(wǎng)格單元)提出了一種改進(jìn)的基于線程嵌套的2層流水線并行LU-SGS算法(簡稱TL_Pipeline[22 - 24])，如圖8所示。其基本思想就是通過將原來的1條長流水線(深度為dp)轉(zhuǎn)換成相對較短的2條嵌套的流水線，從而進(jìn)一步開發(fā)原先每個線程上二維子任務(wù)(子平面)內(nèi)蘊含的并行性。TL_Pipeline外層(第1層)流水線的深度為dp1，在每個流水線階段內(nèi)又包含1個嵌套深度為dp2的內(nèi)層(第2層)流水線，并且dp1×dp2=dp。其中內(nèi)層流水線是在Jsub子平面上構(gòu)造的，子任務(wù)中Jsub維各網(wǎng)格線組成流水線任務(wù)隊列，每一條網(wǎng)格線沿著I維靜態(tài)地剖分成dp2個子網(wǎng)格線。圖8給出的是算法的負(fù)載剖分和任務(wù)調(diào)度示意圖，其中dp1和dp2均為4，即共有16個并行線程。以第1層流水線的K4J1子任務(wù)為例，該子任務(wù)將進(jìn)一步被劃分為若干更細(xì)粒度的子任務(wù)JmIi(i=1,2,3,4)，然后在內(nèi)層流水線上進(jìn)行調(diào)度執(zhí)行。

Figure 8 Illustration of task decomposition and execution timelines for the pipeline approach (left) with 4 threads/pipeline-stages,and the two-level pipeline approach (right) with 4 threads/pipeline-stages in the sub-pipeline圖8 4個線程/流水段時傳統(tǒng)流水線方法(左)和兩層流水線方法(右，子流水線也包含4個流水段)的任務(wù)分解和執(zhí)行過程

4.5 向量化并行

當(dāng)前很多CPU和加速器(例如MIC)均采用了寬向量設(shè)計(例如天河二號CPU的雙精度向量寬度為4，MIC 的雙精度向量寬度為8)，如果沒有充分利用向量化并行，則應(yīng)用的實際浮點性能將下降為峰值性能的1/v(設(shè)v為向量寬度)。對于一些復(fù)雜的CFD計算過程，直接添加向量化指導(dǎo)語句難以實現(xiàn)編譯器自動向量化或者向量化效率較低，需要對CFD計算及訪存特征等進(jìn)行深入分析，必要時對數(shù)據(jù)結(jié)構(gòu)進(jìn)行重構(gòu)并采用intrinsic向量化指令實現(xiàn)高效向量化并行。這里簡單介紹一下作者團(tuán)隊針對高階精度有限差分格式WCNS在MIC架構(gòu)上開展的向量化并行化研究[34]。

圖9給出了5階顯式WCNS格式(WCNS-E-5) 半節(jié)點重構(gòu)模板計算特點。測試表明，250萬網(wǎng)格規(guī)模時半節(jié)點重構(gòu)計算約占了總計算時間的1/3，說明這部分是整個計算的性能熱點。作者團(tuán)隊采用MIC平臺特有的intrinsic向量化指令對WCNS-E-5代碼進(jìn)行了重寫，使用的 intrinsic 語句主要包括_mm512_load_pd(對齊取數(shù)據(jù)，可一次訪存取8個雙精度浮點，有效降低訪存次數(shù))、_mm512_fmadd_pd/_mm512_fmsub_pd(乘加/乘減，可通過運算指令融合提升性能)、_mm512_storenrngo_pd(對齊寫數(shù)據(jù)，將向量寄存器的8個雙精度浮點數(shù)寫入內(nèi)存不緩存，對于只寫訪問有很好的Cache優(yōu)化效果)等。此外，為了使WCNS-E-5更好地適應(yīng)向量化計算，還對相關(guān)數(shù)據(jù)結(jié)構(gòu)做了相應(yīng)的調(diào)整。如圖10所示，首先將數(shù)據(jù)結(jié)構(gòu)由結(jié)構(gòu)體數(shù)組AOS(Array Of Structure)調(diào)整為數(shù)組結(jié)構(gòu)SOA(Structure Of Array)形式，使得數(shù)組能夠大跨度地連續(xù)訪問；其次為了訪存的對齊，對原數(shù)組做擴(kuò)充填補(bǔ)，保證數(shù)組的對齊(只寫)訪問。

Figure 9 WCNS-E-5 interpolation template圖9 WCNS-E-5插值模板

Figure 10 Data structure adjustment圖10 數(shù)據(jù)結(jié)構(gòu)調(diào)整

4.6 實驗結(jié)果

本節(jié)將給出在天河系列超級計算機(jī)上的部分測試結(jié)果。首先定義協(xié)同效率CE(Collaborative Efficiency)以評估CPU-GPU協(xié)同并行中的效率損失:

其中,SPCPU和SPGPU分別是僅實現(xiàn)CPU和GPU并行時的加速比，SPCPU+GPU是CPU-GPU協(xié)同并行獲得的加速比(以SPCPU作為基準(zhǔn))。例如，SPCPU+GPU=1.8,SPGPU=1.3時，協(xié)同效率CE為1.8/(1.0+1.3)×100%≈78.3%，意味著協(xié)同并行中的效率損失約為22%。

表1給出了天河一號超級計算機(jī)的單塊Tesla M2050 GPU上不同流實現(xiàn)策略時的加速比。網(wǎng)格規(guī)模固定為128×128×128，網(wǎng)格分區(qū)數(shù)由2增加到8，以單流實現(xiàn)的性能作為基準(zhǔn)?？梢钥闯?，采用CUDA多流在GPU上同時執(zhí)行多個網(wǎng)格分區(qū)可提升25%～30%的性能。由于CPU-GPU同步以及一些變量的PCI-e傳輸，GBMS有一定的額外開銷，相對于多流有一定的性能損失(最多28%左右)，但作者團(tuán)隊設(shè)計的GBMS策略可以將高精度CFD軟件運行在單個M2050 GPU上的最大模擬容量從2百萬網(wǎng)格單元提升到4百萬網(wǎng)格單元，這為后續(xù)CPU-GPU協(xié)同并行的負(fù)載均衡奠定了基礎(chǔ)。

Table 1 Performance comparison of different CUDA stream implementations表1 GPU上不同流實現(xiàn)策略時的加速比

圖11給出了天河一號超級計算機(jī)單個計算結(jié)點內(nèi)GPU并行、CPU-GPU協(xié)同并行以及實現(xiàn)GBMS策略時的加速比和協(xié)同并行效率。由于GPU存儲容量限制，僅針對前4個相對較小規(guī)模的網(wǎng)格給出了GPU并行加速比，對于網(wǎng)格規(guī)模256×128×128(約4百萬網(wǎng)格單元)，必須采用GBMS策略，SPGPU由1.3降到1.05，對于前面4個小規(guī)模網(wǎng)格問題，協(xié)同并行不需要GBMS，SPCPU+GPU和CE分別可達(dá)到1.8和79%，協(xié)同并行相對于純GPU并行(GPU-only)提升了約45%的性能。對于更大規(guī)模的網(wǎng)格(大于4百萬網(wǎng)格單元)，GBMS對于提升協(xié)同并行加速比和協(xié)同效率非常重要。以256×256×18(約8百萬網(wǎng)格單元)網(wǎng)格規(guī)模為例，如果不采用GBMS策略，則GPU只能模擬其中的2百萬網(wǎng)格單元，其他的6百萬網(wǎng)格單元只能在CPU上進(jìn)行模擬，由于嚴(yán)重的負(fù)載不均衡，SPCPU+GPU僅為1.3，CE僅為57%。采用GBMS策略后，CPU和GPU均可以模擬4百萬網(wǎng)格單元，SPCPU+GPU和CE分別提高到1.79和89%。采用GBMS負(fù)載均衡情況下，高精度CFD軟件在一個天河一號計算結(jié)點上的模擬容量由3.5百萬網(wǎng)格單元提升到8百萬網(wǎng)格單元，提升了約2.3倍。進(jìn)一步增加問題規(guī)模會導(dǎo)致協(xié)同并行加速比和效率明顯下降，原因是即使采用GBMS，GPU模擬負(fù)載最多仍然為4百萬網(wǎng)格單元。

Figure 11 Intra-node collaborative speedup and efficiency圖11 計算結(jié)點內(nèi)協(xié)同并行加速比及效率

圖12給出了天河二號MIC加速器(Intel Xeon Phi 31S1P)上針對LU-SGS設(shè)計的2層流水線OpenMP并行算法TL-Pipeline與傳統(tǒng)流水線OpenMP并行算法的加速比。可以看出，傳統(tǒng)流水線并行的最大加速比僅為15.5(此時對應(yīng)的問題規(guī)模為64×64×64，線程數(shù)為57)，當(dāng)使用MIC全部的228個線程時，加速比降為5.0(32×64×128),7.1(64×64×64)和6.2(128×64×32)，表明在眾核加速器上，隨著線程數(shù)量增加，傳統(tǒng)流水線并行存在嚴(yán)重的可擴(kuò)展性瓶頸。與此同時，對于上面3個問題規(guī)模，即使采用全部228個線程，TL-Pipeline加速比也可分別達(dá)到69.1,82.3和98.3。注意到TL-Pipeline加速比的提升對于不同dp1×dp2的組合變化較大，這是因為流水線開銷以及線程/流水段負(fù)載均衡隨不同模擬變化較大?？傊琓L-Pipeline為在MIC這樣的眾核加速器上實現(xiàn)LU-SGS等高數(shù)據(jù)依賴求解方法的可擴(kuò)展OpenMP并行提供了新的方法。

Figure 12 Improvements of TL-Pipeline over the traditional pipeline approach for three grids dimensions on Intel Xeon Phi圖12 3種網(wǎng)格規(guī)模下MIC上2層流水線并行相對于傳統(tǒng)流水線并行的性能提升

本文采用規(guī)模為256×256×256的算例，在天河二號的MIC協(xié)處理器上測試對比了WCNS半節(jié)點重構(gòu)計算的基準(zhǔn)版本和intrinsics向量化優(yōu)化版本的性能改進(jìn)(如圖13所示)。基準(zhǔn)版本在MIC端當(dāng)線程數(shù)為224的時候性能達(dá)到最優(yōu)，相對MIC上單線程計算的加速比約為83.8倍。采用深度intrinsics向量化優(yōu)化實現(xiàn)之后，在線程數(shù)規(guī)模不大的情況下，intrinsics優(yōu)化實現(xiàn)同樣具有很好的線程可擴(kuò)展性。得益于MIC更寬的向量位寬，深度intrinsics優(yōu)化將最優(yōu)性能提升到基準(zhǔn)實現(xiàn)的4.5倍左右，在112線程時取得最大加速比，說明向量優(yōu)化之后MIC上的寬向量部件性能得到了更充分的開發(fā)。

Figure 13 Performance comparison of baseline and intrinsic implementation for WCNS-E-5 interpolation on MIC圖13 WCNS半節(jié)點重構(gòu)基準(zhǔn)版本和intrinsic版本在MIC上的性能對比

5 結(jié)束語

當(dāng)前，天河系列、神威太湖之光等超級計算機(jī)已多次排名世界第一，標(biāo)志著我國超級計算機(jī)系統(tǒng)研制能力已進(jìn)入世界前列。與此同時，國內(nèi)很多單位在這些國產(chǎn)超級計算機(jī)上開展了大量CFD并行應(yīng)用軟件開發(fā)和性能優(yōu)化工作，取得了不錯的成果。長期以來，國防科技大學(xué)CFD應(yīng)用軟件團(tuán)隊以天河/銀河系列超級計算機(jī)為依托，開展了超大規(guī)模復(fù)雜CFD并行應(yīng)用開發(fā)和性能優(yōu)化研究，突破了新型異構(gòu)協(xié)同并行計算等一系列關(guān)鍵技術(shù)，初步實現(xiàn)了HPC與CFD的深度融合，有力支撐了國產(chǎn)CFD軟件在天河/銀河系列超級計算機(jī)上的高效超大規(guī)模并行應(yīng)用。

當(dāng)前高性能計算發(fā)展的下一個里程碑是E級計算(Exascale Computing，百億億次浮點運算/每秒)，美、日、歐洲、俄羅斯等都制定了E級超級計算機(jī)研制計劃。例如，2016年美國能源部正式啟動了E級計算計劃ECP(Exascale Computing Project)[35]。ECP特別強(qiáng)調(diào)應(yīng)用軟件開發(fā)，將其作為第1個重點關(guān)注領(lǐng)域。我國也將E級超級計算機(jī)研制納入了國家“十三五”規(guī)劃，目前國防科技大學(xué)、江南計算技術(shù)研究所、中科曙光3家單位牽頭開展的E級超算原型系統(tǒng)研制項目已通過驗收，實際的國產(chǎn)E級超級計算機(jī)系統(tǒng)預(yù)期在2021年左右研制成功。我國E級計算計劃中同樣非常關(guān)注與E級計算機(jī)配套的高性能計算應(yīng)用軟件研制，國家重點研發(fā)計劃中已部署了一批應(yīng)用軟件研發(fā)項目，其中包括數(shù)值飛行器原型、數(shù)值發(fā)動機(jī)原型等CFD相關(guān)項目。可以預(yù)見，基于CFD與E級計算機(jī)融合的“數(shù)值風(fēng)洞試驗”“數(shù)值優(yōu)化設(shè)計”“數(shù)值虛擬飛行”將給航空航天飛行器設(shè)計帶來革命性的變化,并將推動流體力學(xué)和空氣動力學(xué)等學(xué)科的創(chuàng)新發(fā)展[36]。當(dāng)前，E級系統(tǒng)CFD應(yīng)用開發(fā)面臨巨大挑戰(zhàn)。一方面超級計算機(jī)體系結(jié)構(gòu)異構(gòu)、眾核、寬向量趨勢明顯，目前大多數(shù)CFD應(yīng)用軟件只能在純CPU系統(tǒng)上運行，通常僅支持MPI并行計算，尚不具備利用新型異構(gòu)眾核寬向量并行體系結(jié)構(gòu)的能力，難以充分發(fā)揮超級計算機(jī)的多層次異構(gòu)并行性能潛力。另一方面，CFD數(shù)值模擬實際浮點計算性能不高、并行可擴(kuò)展性差的情況仍然普遍存在。對真實復(fù)雜CFD應(yīng)用而言，機(jī)器浮點效率常常低于10%甚至5%，擴(kuò)展到千核以上并行效率嚴(yán)重下降。如何高效地利用大量寬向量計算核心和異構(gòu)體系結(jié)構(gòu)，獲得實際計算的高性能，是一個嚴(yán)峻的挑戰(zhàn)，“應(yīng)用墻”問題依然突出。因此，迫切需要從大規(guī)模CFD數(shù)值模擬應(yīng)用的數(shù)值模型和算法特點出發(fā)，緊密結(jié)合新型異構(gòu)眾核體系結(jié)構(gòu)特征，針對性地開展并行算法研究工作，使應(yīng)用程序充分發(fā)掘大規(guī)模新一代并行計算機(jī)性能，支撐實際CFD應(yīng)用的高效能計算，滿足國家航空航天飛行器等重大型號工程氣動設(shè)計需求。