賴劍奇，李樺，張冉，常青

國防科技大學(xué) 空天科學(xué)學(xué)院，長沙 410073

隨著計算流體力學(xué)(Computational Fluid Dynamics, CFD)的日益成熟和計算機技術(shù)的迅猛發(fā)展，CFD在飛行器氣動性能分析與優(yōu)化設(shè)計、復(fù)雜流動機理分析中的作用和地位日益突出[1-2]。對于自然界中普遍存在的高雷諾數(shù)流動如湍流、漩渦、分離、氣動聲學(xué)等多尺度復(fù)雜流動問題，密集的計算網(wǎng)格往往不可或缺。對于真實氣體流動的模擬，包含多種組分的化學(xué)反應(yīng)模型在一定程度上會使得計算網(wǎng)格成倍增加[3]。同時在工程實際中，對于飛行器的全機模擬，一般要求單個計算狀態(tài)的計算網(wǎng)格在千萬量級甚至上億的水平，并且需要計算的狀態(tài)數(shù)目成百上千，如此龐大的計算量迫切需要提高計算效率。因此，快速精確地實現(xiàn)CFD的數(shù)值模擬，是廣大科研工作者共同追求的目標(biāo)。

傳統(tǒng)上,開展基于CPU集群的并行計算是解決大規(guī)模科學(xué)計算問題的有效途徑，而并行計算集群的性能在很大程度上依賴于CPU的更新?lián)Q代。但是隨著CPU單位面積內(nèi)集成的晶體管數(shù)量增多，能量消耗和散熱問題不可忽視，導(dǎo)致CPU的更新速度變緩，發(fā)展陷入瓶頸。圖形處理單元(Graphics Processing Units，GPU)在數(shù)據(jù)并行上具有強大的浮點運算能力和存儲帶寬[4]，其作為一種新型的科學(xué)計算設(shè)備開始進入人們的視野。與Intel CPU相比，GPU具有突出的性能優(yōu)勢。GPU近年來在分子動力學(xué)(MD)、CFD、氣象預(yù)測(WF)、人工智能(AI)等通用計算領(lǐng)域應(yīng)用廣泛[5-9]。GPU的發(fā)展給CFD高性能計算領(lǐng)域帶來了前所未有的機遇和挑戰(zhàn)，在GPU上實現(xiàn)CFD大規(guī)模并行計算一直是GPU在通用計算領(lǐng)域的熱點，將會對CFD的研究與應(yīng)用產(chǎn)生巨大的推動作用。

傳統(tǒng)上，GPU僅限于圖形渲染領(lǐng)域。NVIDIA在2007年推出了統(tǒng)一計算設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)，降低了編譯程序的復(fù)雜性。CUDA作為編程模型的出現(xiàn)標(biāo)志著GPU開始在通用計算領(lǐng)域廣泛應(yīng)用。在CFD領(lǐng)域，國內(nèi)外學(xué)者開始將GPU應(yīng)用于Euler方程和Navier-Stokes方程的求解，并廣泛應(yīng)用于層流和湍流物理流動規(guī)律的分析，取得了一系列重要成果。Brandvik和Pullan[10]首次在NVIDIA 8800GTX GPU上實現(xiàn)三維Euler方程的顯式求解，計算加速比達到16倍。Jacobsen等[11]在多GPU并行計算集群上建立了基于消息傳遞接口+統(tǒng)一計算設(shè)備架構(gòu)(MPI+CUDA)的不可壓流求解器，并通過方腔頂蓋驅(qū)動流動算例對多GPU并行系統(tǒng)的性能進行了分析。Castonguay等[12]在Telsa C2050 GPU上實現(xiàn)了基于非結(jié)構(gòu)網(wǎng)格的可壓縮黏性流動求解，計算問題復(fù)雜性進一步提高。Emelyanov等[13]在GPU上求解超聲速平板流動問題，計算加速比達到45倍，同時分析了影響GPU計算效率的因素，提出了優(yōu)化方式。Watkins等[14]在多GPU上實現(xiàn)了高階格式的隱式求解，進一步拓寬了GPU的應(yīng)用范圍。Aissa等[15]對在GPU上實現(xiàn)顯式和隱式推進求解定常流動問題進行了對比，并對算法的收斂性進行了分析。國內(nèi)宋慎義等[16]在GPU上實現(xiàn)了基于非結(jié)構(gòu)網(wǎng)格CFD求解器的設(shè)計，討論了科學(xué)計算程序的優(yōu)化方法。徐傳福等[17-19]在天河-1A巨型機上基于WCNS格式并行求解了復(fù)雜外形流動問題。對于C919大飛機外形，計算網(wǎng)格達到了1.5億。李大力等[20]在天河2超級計算機上實現(xiàn)了GPU的高階格式隱式求解。馬文鵬等[21]在GPU上實現(xiàn)了基于混合網(wǎng)格的非定常流動問題的求解，并針對NACA0012等簡單外形算例進行了分析。劉楓等[22]建立了基于MPI+CUDA的異構(gòu)并行可壓縮流求解器，并對求解器的性能進行了測試。數(shù)值結(jié)果顯示求解器魯棒性好，計算效率較CPU同構(gòu)并行提高10倍以上。曹文斌等[23]應(yīng)用多GPU對可壓縮湍流開展并行計算，取得了良好的加速效果。國內(nèi)外在GPU上實現(xiàn)CFD的大規(guī)模并行計算方面開展了大量的工作，但是對于算法的性能缺乏系統(tǒng)的評價。

本文從GPU的架構(gòu)特點和CUDA編程模型出發(fā)，在NVIDIA GTX 1070 GPU上建立了基于MPI+CUDA的多GPU并行可壓縮流求解器。針對超聲速進氣道算例，對求解器的單GPU并行性能和多GPU可擴展性能進行分析。對算法的性能進行系統(tǒng)的評價，對于在多GPU上開展CFD大規(guī)模并行計算具有重要意義。

1 控制方程和數(shù)值計算方法

1.1 控制方程

在不考慮外部加熱和徹體力等源項的影響時，積分形式的可壓縮三維流動Navier-Stokes方程為[24]

(1)

式中：W為守恒變量矢量；Fc為無黏通量矢量；Fv為黏性通量矢量；Ω為控制體單元；?Ω為控制體單元的邊界；dΩ為體積微元；dS為面積微元。

如圖1所示，采用基于結(jié)構(gòu)網(wǎng)格的格心式有限體積法對式(1)進行離散。圖中：n為控制體表面外法線單位向量；(I,J)為控制體單元的中心；(i,j)為控制體單元的頂點。

對于控制體ΩI,J,K，對式(1)進行離散可得

(2)

式中：NF為控制體表面數(shù)量；ΔSm為控制體表面面積。定義式(2)中的右端項為殘值，記為RI,J,K，則有

(3)

1.2 數(shù)值計算方法

對于無黏通量的計算，本文采用AUSM+UP格式[25]。AUSM+UP格式是一種應(yīng)用較為廣泛的上風(fēng)格式，具有較高的間斷分辨率和計算效率。基于質(zhì)量流量函數(shù)將無黏通量分裂為對流項和壓力項：

(4)

(5)

采用MUSCL (Monotone Upstream-centered Schemes for Conservation Law)插值方法[26]獲得二階精度迎風(fēng)格式開展數(shù)值計算，通過Van Albada限制器[27]抑制間斷處振蕩和非物理解的產(chǎn)生。

對于黏性通量的計算，需要用到控制體單元界面上的物理量及其一階導(dǎo)數(shù)值?？紤]到黏性項的橢圓型數(shù)學(xué)性質(zhì)，界面上的物理量由兩側(cè)單元中心值的平均值得到,即

UI+1/2=(UI+UI+1)/2

(6)

一階導(dǎo)數(shù)值由高斯積分公式進行求解，即

(7)

時間推進采用顯式多步Runge-Kutta格式[28]，該格式只存儲最初的守恒變量和最近的殘值，可以減少存儲開銷，并具有良好的數(shù)據(jù)并行性。表達式為

(8)

式中：Δt為時間步長；αk為k階系數(shù)。

湍流模型采用k-ω剪切應(yīng)力輸運(SST)兩方程模型，與式(1)相比，方程的形式保持不變，僅增加了湍流源項。此時黏性系數(shù)μ由層流黏性系數(shù)μL和渦黏性系數(shù)μT兩部分組成，即

μ=μL+μT

(9)

式中：μL通過Sutherland公式得到，μT由湍流模型給出，當(dāng)不考慮湍流模型時，μT的值為0。

2 多GPU并行計算求解器

2.1 CUDA編程模型

GPU最初用于圖形處理，專門執(zhí)行復(fù)雜的數(shù)學(xué)和幾何計算。隨著CUDA、OpenCL、OpenACC等通用編程模型的推廣，GPU被用來進行通用數(shù)值計算，稱為GPGPU(General Purpose Graphics Processing Units)[3, 29-30]。本文測試所用的GPU為NVIDIA GTX 1070，其主要性能參數(shù)如表1所示。從表中可以看出GTX 1070具有強大的浮點運算能力和存儲帶寬，在并行求解可壓縮Navier-Stokes方程時具有較大的優(yōu)勢。由于GTX 1070的單精度浮點計算能力遠(yuǎn)大于雙精度浮點計算能力，GPU計算使用單精度。

CUDA支持C/C++、Fortran等語言的擴展，具有較強的通用性。完整的CUDA程序包含主機端和設(shè)備端兩部分，主機端代碼在CPU上執(zhí)行，設(shè)備端代碼調(diào)用內(nèi)核在GPU上執(zhí)行[31]。內(nèi)核是GPU上的函數(shù)，對應(yīng)于線程網(wǎng)格，一個線程網(wǎng)格由若干個線程塊組成。對于GTX 1070，一個線程塊至多包含1 024個線程。在CUDA中，線程塊數(shù)目和線程塊中的線程數(shù)目對GPU計算性能具有重要影響。通過測試，當(dāng)線程塊中的線程數(shù)目設(shè)置為256時，GTX 1070的計算性能可以達到最優(yōu)。而線程塊的數(shù)目需要根據(jù)計算規(guī)模來確定，在計算過程中，確保每個線程負(fù)載一個網(wǎng)格單元。

表1 GTX 1070的主要性能參數(shù)Table 1 Main performance parameters of GTX 1070

2.2 GPU并行模式

傳統(tǒng)上基于CPU的并行計算采用物理網(wǎng)格分塊的方法，將計算域劃分為若干個網(wǎng)格量相當(dāng)?shù)淖訁^(qū)域，然后每個CPU核心負(fù)載一個子區(qū)域的計算量。在編程上采用單控制流多數(shù)據(jù)流(SPMD)模型，采用消息傳遞接口(MPI)實現(xiàn)數(shù)據(jù)的交換及同步。多CPU并行計算效率很大程度上依賴于中央處理器單元的性能，考慮到CPU發(fā)展速度變緩及CPU之間的通信帶寬，不斷增加CPU的數(shù)目并不能增加并行計算效率[32]。

GPU并行模式與多CPU并行模式存在較大的不同。以GTX 1070為例，其包含15個多流處理器(SM)，同一個線程塊內(nèi)的線程保證在同一個SM上同時執(zhí)行，并通過共享存儲器實現(xiàn)數(shù)據(jù)的快速訪存。SM的架構(gòu)稱之為單指令多線程(SIMT)，與單指令多數(shù)據(jù)(SIMD)架構(gòu)類似。GPU具有強大的計算能力，但其存儲能力較弱，而CPU卻相反，充分發(fā)揮CPU與GPU的優(yōu)勢，實現(xiàn)CPU與GPU的合理分工具有重要意義。圖2為GPU并行模式示意圖。CPU負(fù)責(zé)流場的初始化和后處理，主機與設(shè)備之間的數(shù)據(jù)傳輸，內(nèi)核程序的啟動、完成與同步等管理性事物，而GPU內(nèi)核程序則負(fù)責(zé)處理與網(wǎng)格相關(guān)的數(shù)值計算。對于CPU與GPU之間的數(shù)據(jù)傳輸，需要用到全局內(nèi)存，其位于顯存中，訪問延遲很大，因此只針對原始變量進行數(shù)據(jù)傳輸，本文選取的原始變量為U=[ρuvwT]T。

2.3 MPI+CUDA多GPU并行算法

對于在多GPU上實現(xiàn)并行計算，目前應(yīng)用最廣泛的兩種應(yīng)用程序接口(API)是OpenMP和MPI，形成了許多OpenMP+CUDA、MPI+CUDA和MPI+OpenMP+CUDA并行程序包[33-35]。OpenMP基于共享式存儲模型實現(xiàn)單節(jié)點內(nèi)部的消息傳遞和同步，其傳遞效率較高。而MPI基于分布式存儲模型或共享式存儲模型實現(xiàn)并行編程的消息傳遞和同步，其可用于節(jié)點內(nèi)部和節(jié)點之間。雖然MPI的傳遞效率略低于OpenMP，但是其編程性能良好，對于實現(xiàn)程序的可擴展性具有重要意義，因此MPI是目前應(yīng)用最廣泛的并行編程消息傳遞庫。本文建立了基于MPI+CUDA編程模型的多GPU并行可壓縮流求解器，實現(xiàn)了GPU并行計算程序的擴展。

圖3為MPI+CUDA編程模型示意圖，通過MPI實現(xiàn)多節(jié)點多GPU并行計算。在進行設(shè)備間的數(shù)據(jù)交換時，GPU間的數(shù)據(jù)傳遞不能直接完成，需要通過CPU作為“媒介”來完成。對于位于不同設(shè)備上的相鄰網(wǎng)格區(qū)域，首先將其邊界數(shù)據(jù)通過PCI-e總線傳遞至主機端，然后調(diào)用MPI，在交換池中實現(xiàn)數(shù)據(jù)交換，并通過PCI-e總線傳遞至設(shè)備端。GPU間的數(shù)據(jù)交換示意圖如圖4所示。目前一些計算設(shè)備架構(gòu)采用NVLink技術(shù)，能夠在CPU與GPU間和GPU間實現(xiàn)超高速的數(shù)據(jù)傳輸，傳輸速度是傳統(tǒng)PCI-e總線速度的5～12倍，但現(xiàn)有的計算設(shè)備并不支持該項技術(shù)。

在多GPU的CFD并行計算中，考慮多個計算節(jié)點間及節(jié)點內(nèi)各GPU之間的負(fù)載平衡問題，需要使各計算設(shè)備上負(fù)載的網(wǎng)格量一致。常用的網(wǎng)格分區(qū)方法有一維區(qū)域分解法和三維區(qū)域分解法[11, 29]，與三維區(qū)域分解技術(shù)相比，一維區(qū)域分解技術(shù)雖然增加了通信量，但是由于數(shù)據(jù)傳遞的連續(xù)性使得其傳遞效率高，因此本文網(wǎng)格分區(qū)使用一維區(qū)域分解方法，如圖5所示。

受到計算資源的限制，本文僅針對單節(jié)點并行系統(tǒng)進行測試，每個節(jié)點包含1顆CPU和4塊GPU。測試使用的CPU為Intel Xeon E5-2670，八核，主頻為2.6 GHz，主機內(nèi)存容量為64 GB；GPU為NVIDIA GTX 1070，其性能參數(shù)見表1。本文編程采用CUDA 7.5，C語言編譯采用Intel icc 11.1編譯器，MPI采用MPICH2 1.4.1。

3 算例分析

基于本文所建立的多GPU并行計算求解器，針對CFD中可壓縮流動Navier-Stokes方程的求解，對超聲速進氣道算例進行了計算，分析了求解器的單GPU并行性能和多GPU可擴展性能。

3.1 計算模型與條件

Reinartz等[36]對混合壓縮進氣道進行了大量的數(shù)值與實驗分析，獲得了清晰的流場結(jié)構(gòu)和壁面壓力實驗數(shù)據(jù)。超聲速進氣道的流場結(jié)構(gòu)復(fù)雜，包含邊界層、流動分離、膨脹區(qū)、反射激波和邊界層誘導(dǎo)激波等復(fù)雜的物理現(xiàn)象。圖6是進氣道構(gòu)型，進氣道總長400 mm，高52 mm，隔離段總長l為79.3 mm，隔離段高度h為15 mm，上楔面傾角δ1為21.5°，下楔面傾角δ2為9.5°，擴張段擴張角δ3為5°。計算條件：來流馬赫數(shù)為2.41，單位雷諾數(shù)為5.07×107，總壓為540 kPa，總溫為305 K，攻角為-10°。

為了對求解器的并行性能和可擴展性能進行分析，設(shè)計了6套網(wǎng)格，如表2所示。用Pointwise軟件生成結(jié)構(gòu)網(wǎng)格，各套網(wǎng)格在壁面和拐角處進行了加密，第1層網(wǎng)格高度為1×10-6m，局部網(wǎng)格如圖7所示。

序號網(wǎng)格大小網(wǎng)格量/104grid1 601×129×11 76.8grid21 201×129×11 153.6grid31 201×257×11 307.2grid42 401×257×11 614.4grid52 401×513×111 228.8grid64 801×513×112 457.6

3.2 計算結(jié)果

圖8是4套網(wǎng)格條件下的壁面壓力分布，p/pt,0為壓力與來流總壓的比值。從圖中可以看出，不同規(guī)模網(wǎng)格對計算結(jié)果影響不大，本文選取grid3作為計算網(wǎng)格。圖9是壁面壓力計算結(jié)果與實驗數(shù)據(jù)的對比，從圖中可以看出，計算結(jié)果和實驗數(shù)據(jù)吻合較好。圖10是馬赫數(shù)等值線圖，數(shù)值計算清晰地捕捉了流場內(nèi)的復(fù)雜波系結(jié)構(gòu)。

3.3 性能分析

加速比(SP)和并行效率(CE)是衡量并行算法性能的重要參數(shù)，也是研究并行算法可擴展性能的重要指標(biāo)。加速比定義為單顆CPU平均計算時間與相應(yīng)GPU時間的比值,即

SP=tCPU/tGPU

(10)

并行效率定義為單塊GPU處理相應(yīng)負(fù)載所需的平均計算時間與多GPU處理相應(yīng)倍數(shù)負(fù)載時間的比值,即

CE=ts/tp

(11)

在計算過程中，通過模擬10 000時間步得到單個時間步的平均計算時間。

表3給出了不同網(wǎng)格規(guī)模條件下CPU和單GPU計算時間。圖11給出了單GPU并行計算加速比隨網(wǎng)格規(guī)模的變化曲線，從圖中可以看出GPU并行計算加速比達到了37倍以上。隨著網(wǎng)格規(guī)模的增大，單GPU加速比從最初的37倍逐漸增加到46倍。這是因為數(shù)據(jù)交換等分支指令任務(wù)與迭代計算等單一指令任務(wù)相比，占計算任務(wù)的比例隨網(wǎng)格規(guī)模的增大而減小，說明GPU非常適合求解大網(wǎng)格規(guī)模的復(fù)雜運算空間。當(dāng)網(wǎng)格規(guī)模較小時，加速比隨網(wǎng)格規(guī)模的增加迅速增大，并逐漸趨于平緩。當(dāng)網(wǎng)格規(guī)模增加到飽和值時，加速比基本保持不變并達到加速比極限。飽和值和加速比極限與GPU架構(gòu)特點有關(guān)，對于本文所用的NVIDIA GTX 1070，飽和值約為1 200萬，加速比極限約為46。

并行算法具有良好的可擴展性能是開展大規(guī)模并行計算的基礎(chǔ)，本文從強擴展性和弱擴展性兩方面對并行算法的可擴展性能進行評價。表4給出了不同網(wǎng)格規(guī)模條件下多GPU并行計算的時間。圖12是單節(jié)點并行算法強擴展性分析示意圖，整個節(jié)點上負(fù)載的網(wǎng)格規(guī)模保持一致。從圖中可以看出在網(wǎng)格規(guī)模較小時，多GPU并行計算并沒有顯示出明顯的優(yōu)勢，隨著網(wǎng)格規(guī)模的增大，多GPU并行計算能夠極大地提高計算效率，證明本文建立的并行算法具有良好的強擴展性。例如對于grid6，單GPU并行計算加速比為46，兩塊GPU為84，而4塊GPU則達到143。在多GPU并行計算中，隨著網(wǎng)格規(guī)模的增大，與計算時間相比，通信時間所占的比例逐漸減小，因此更能發(fā)揮多GPU并行計算的優(yōu)勢。

表3 不同規(guī)模網(wǎng)格計算時間Table 3 Computing time for different number of grids

序號2 GPUs/ms4 GPUs/msgrid1 13.5 12.7grid2 21.8 16.6grid3 39.1 26.0grid4 75.6 44.5grid5146.3 85.5grid6283.9166.4

圖13是單節(jié)點并行系統(tǒng)弱擴展性分析示意圖，整個節(jié)點上每塊GPU上負(fù)載的網(wǎng)格規(guī)模保持一致。從圖中可以看出隨著參與計算GPU數(shù)量的增加，并行效率有所下降，但依舊維持在一個較高的值，證明本文建立的并行算法具有良好的弱擴展性。例如對于grid4，兩塊GPU并行效率達到90.18%，4塊GPU并行效率達到77.83%。與兩塊GPU并行計算相比，4塊GPU并行計算雖然能提高計算效率，但其并行效率有所降低。參與并行計算的GPU設(shè)備越多，通信所花費的時間也隨之增加，導(dǎo)致并行效率下降。根據(jù)計算規(guī)模的不同，合理調(diào)用計算節(jié)點對于充分發(fā)揮算法在大規(guī)模并行計算上的優(yōu)勢具有重要的意義。

4 結(jié) 論

1) 在NVIDIA GTX 1070 GPU上建立基于MPI+CUDA的多GPU并行可壓縮流求解器，應(yīng)用本文所建立的求解器求解超聲速進氣道算例，計算結(jié)果與實驗數(shù)據(jù)吻合良好，求解器具有較好的準(zhǔn)確性。

2) 本文建立的求解器可以充分發(fā)揮GPU在數(shù)據(jù)并行上的強大浮點運算能力和存儲帶寬，對于超聲速進氣道算例，GPU并行計算取得了37倍以上的加速比。隨著網(wǎng)格規(guī)模的增加，單GPU加速比從最初的37倍逐漸增加到46倍。同時單GPU并行計算加速比隨著網(wǎng)格規(guī)模的增加而逐漸趨于平緩，當(dāng)網(wǎng)格規(guī)模增加到“飽和值”時，加速比保持不變。

3) 對于多GPU并行算法，對其可擴展性能進行分析。隨著網(wǎng)格規(guī)模的增大，多GPU并行計算能夠極大地提高計算效率，證明本文建立的并行算法具有良好的強擴展性；同時隨著參與計算的GPU數(shù)量的增加，并行效率有所下降，但依舊維持在一個較高的值，證明本文建立的并行算法具有良好的弱擴展性。本文建立的多GPU可壓縮流求解器可擴展性能良好，可以推廣至CFD的大規(guī)模并行計算中。