多GPU并行可壓縮流求解器及其性能分析

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

國防科技大學 空天科學學院，長沙 410073

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

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

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

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

1 控制方程和數值計算方法

1.1 控制方程

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

(1)

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

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

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

(2)

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

(3)

1.2 數值計算方法

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

(4)

(5)

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

對于黏性通量的計算，需要用到控制體單元界面上的物理量及其一階導數值。考慮到黏性項的橢圓型數學性質，界面上的物理量由兩側單元中心值的平均值得到,即

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

(6)

一階導數值由高斯積分公式進行求解，即

(7)

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

(8)

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

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

μ=μL+μT

(9)

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

2 多GPU并行計算求解器

2.1 CUDA編程模型

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

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

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

2.2 GPU并行模式

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

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

2.3 MPI+CUDA多GPU并行算法

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

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

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

受到計算資源的限制，本文僅針對單節點并行系統進行測試，每個節點包含1顆CPU和4塊GPU。測試使用的CPU為Intel Xeon E5-2670，八核，主頻為2.6 GHz，主機內存容量為64 GB；GPU為NVIDIA GTX 1070，其性能參數見表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]對混合壓縮進氣道進行了大量的數值與實驗分析，獲得了清晰的流場結構和壁面壓力實驗數據。超聲速進氣道的流場結構復雜，包含邊界層、流動分離、膨脹區、反射激波和邊界層誘導激波等復雜的物理現象。圖6是進氣道構型，進氣道總長400 mm，高52 mm，隔離段總長l為79.3 mm，隔離段高度h為15 mm，上楔面傾角δ1為21.5°，下楔面傾角δ2為9.5°，擴張段擴張角δ3為5°。計算條件：來流馬赫數為2.41，單位雷諾數為5.07×107，總壓為540 kPa，總溫為305 K，攻角為-10°。

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

序號網格大小網格量/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 計算結果

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

3.3 性能分析

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

SP=tCPU/tGPU

(10)

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

CE=ts/tp

(11)

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

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

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

表3 不同規模網格計算時間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是單節點并行系統弱擴展性分析示意圖，整個節點上每塊GPU上負載的網格規模保持一致。從圖中可以看出隨著參與計算GPU數量的增加，并行效率有所下降，但依舊維持在一個較高的值，證明本文建立的并行算法具有良好的弱擴展性。例如對于grid4，兩塊GPU并行效率達到90.18%，4塊GPU并行效率達到77.83%。與兩塊GPU并行計算相比，4塊GPU并行計算雖然能提高計算效率，但其并行效率有所降低。參與并行計算的GPU設備越多，通信所花費的時間也隨之增加，導致并行效率下降。根據計算規模的不同，合理調用計算節點對于充分發揮算法在大規模并行計算上的優勢具有重要的意義。

4 結 論

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

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

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