基于異構平臺的通量分裂格式性能研究

梁正虹，黃 俊，劉志勤，陳 波，楊 茂

(西南科技大學 計算機科學與技術學院，四川 綿陽 621010)

0 引言

計算流體力學(CFD,computational fluid dynamics)通過數(shù)值求解流體動力學控制方程，獲取各種條件下的流動數(shù)據(jù)和作用在繞流物體上的力、力矩和熱量等,從而達到研究各種流動現(xiàn)象和規(guī)律的目的。CFD是涉及流體力學、計算機科學與技術、計算數(shù)學等多個專業(yè)的交叉研究領域[1]。在航空航天領域，CFD已成為獲取高超聲速飛行器氣動力和氣動熱數(shù)據(jù)、開展高超聲速流動基礎科學問題研究的三大手段之一[2]。

對于高超聲速流動，保持格式的迎風特性對于計算結果十分有益。迎風格式主要沿著矢通量分裂(FVS, flux vector splitting)和通量差分分裂(FDS,flux difference splitting)兩個路線發(fā)展。FVS的典型格式有Steger-Warming[2]和Van Leer[2]；FDS最具代表性的是Roe[2]格式。FVS將對流項通量按其特征值的正負進行分裂，穩(wěn)定性好但數(shù)值粘性較大；FDS通過近似求解Riemann問題獲得控制單元通量，計算量較大且需要引入熵修正，增大了數(shù)值粘性[3]。20世紀90年代，Liou和Steffen提出兼具FVS和FDS優(yōu)點的AUSM[4](advection upstream splitting method)格式，經(jīng)過20多年的發(fā)展，其衍生格式已成為目前航空航天領域應用最廣泛的格式之一。

CPU/GPU異構并行架構是當今構建大規(guī)模計算集群的主要架構之一，計算機技術的快速發(fā)展給CFD高性能計算領域帶來前所未有的機遇與挑戰(zhàn)。美國國家航天局(NASA)預測，21世紀，高效能計算機和CFD技術的進一步結合將給各類航空航天飛行器的氣動設計帶來一場革命[5]。近十年來，國內外專家學者通過CPU/GPU異構平臺加速CFD應用，取得一系列重要成果。D.A.Jacobsen[6]等成功地將CFD應用從單GPU擴展至多GPU和集群；M.Aissa[7]等在GPU平臺上開展了基于結構網(wǎng)格的顯格式隱格式加速性能研究；Y.Xia等[8-9]研究結構網(wǎng)格、非結構網(wǎng)格和混合網(wǎng)格在GPU上實現(xiàn)的差異；I.C.Kampolis等[10]分析了GPU架構單精度和雙精度對仿真結果的影響。馬文鵬[11]等在GPU上實現(xiàn)了基于混合網(wǎng)格的非定常流動問題的求解；劉楓[12]等建立了基于MPI+CUDA的異構并行可壓縮流求解器；李大力[13]等在天河2超級計算機上實現(xiàn)了GPU的高階格式隱式求解。國內外學者在CPU/GPU異構平臺上開展大規(guī)模CFD并行計算做出大量相關工作，但鮮有文獻對不同通量分裂方法的典型格式在CPU/GPU異構體系下的性能做出具體分析。

本文以一維激波管問題為算例，給出通量分裂方法的典型格式在GPU架構下的具體性能分析，為進一步在CPU/GPU異構平臺上開展大規(guī)模CFD工程應用提供有益的指導和參考。

1 通量分裂方法

在進行CFD數(shù)值模擬中，差分格式的構造是CFD的關鍵之一。迎風類差分格式是從20世紀80年代開始興起，在構造時就體現(xiàn)了方程在波動和流量等傳播方向上的物理特性，與中心差分格式相比具有表現(xiàn)具體流動物理特征的天然優(yōu)勢。發(fā)展至今，各種迎風類格式已經(jīng)成為工程界和學術界研究的主流。根據(jù)對方程物理特性的不同表現(xiàn)形式，迎風格式也可大致分為矢通量分裂(FVS)方法和通量差分分裂(FDS)方法。

1.1 矢通量分裂

FVS(flux vector splitting)方法，最初發(fā)展于20世紀80年代，它根據(jù)相關傳播速度的正負對通量項進行分裂。矢通量分裂方法，簡單、高效、魯棒性高、理論上不會出現(xiàn)非物理解，并且具有較強線性波(如激波等)捕捉能力。根據(jù)不同的分裂方式，構造出不同的FVS格式，最具代表性是Steger-Warming格式和Van Leer格式。其中由于Van Leer格式在駐點、聲速點等特征值變號的附近區(qū)域可以光滑過度，從而在航空航天工程領域得到廣泛應用。

一維歐拉方程進行矢通量分裂：

(1)

Van Leer格式通量分裂表達式：

Ma>1:

F+=0,F-=F

(2)

|Ma|>1:

(3)

Ma<-1:

F+=F,F-=0

(4)

1.2 通量差分分裂

通量差分分裂是迎風類格式的另一條發(fā)展路線，F(xiàn)DS以Godunov方法為基本思路，但與Godunov精確求解黎曼問題的方法不同，F(xiàn)DS近似求解黎曼問題，其對激波和接觸間斷都具有很高的分辨率。根據(jù)不同的解近似黎曼問題的方法，構造不同的FDS格式，最具代表性的是Roe格式。

將一維守恒形式歐拉方程表示為：

(5)

(6)

Roe格式通量表達式為：

(7)

其中：下標(i+1/2)為控制體表面；下標L和R為控制體表面左右兩側。

1.3 AUSM分裂方法

在20世紀90年代初，Liou和Steffen利用不同格式的優(yōu)點，提出AUSM(advection upstream splitting method) 格式。它結合了矢通量分裂在非線性波捕捉上的魯棒性和通量差分分裂在線性波上的高分辨率。從格式構造來講，AUSM是Van Leer格式的一種改進，但從其耗散項來分析，它是一種FVS與FDS的復合格式。目前，AUSM格式經(jīng)過20多年的發(fā)展，已經(jīng)衍生出一系列格式。本文以高超聲速領域使用較廣的AUSMPW+格式為基礎，研究其在GPU硬件架構下的計算性能。

AUSMPW+格式具體表達式：

(8)

(9)

(10)

(11)

(12)

(13)

(14)

式中，f為基于壓力的權函數(shù)，具體形式參考文獻[2]。

2 Riemann問題與計算模型

2.1 Riemann問題

Riemann問題是CFD的經(jīng)典算例，它包含了CFD中的各種間斷(膨脹波、激波及接觸間斷等)問題，對于這些間斷問題的求解一直是CFD發(fā)展的難點和核心問題。Riemann問題具有精確解，可以驗證數(shù)值方法的精度，由此，一般差分格式都是建立在求解Riemann問題基礎之上，之后再向多維擴展。

一維Riemann問題實質上是激波管問題。激波管是一根兩端封閉、內部充滿氣體的直管。直管中一層薄膜將管隔開，在薄膜兩側充滿均勻理想氣體，薄膜兩側氣體的壓力、密度不同。當時，氣體處理靜止狀態(tài)。當時，薄膜瞬時破裂，氣體從高壓端沖向低壓端，同時在管內形成激波、稀疏波和接觸間斷等復雜波系。激波管問題初始條件簡單，計算網(wǎng)格只需一維等距劃分，計算過程不引入復雜處理、網(wǎng)格質量等非格式因素的影響。

激波管示意圖如圖1所示。

圖1 Sod激波管問題

對于本次計算問題，計算區(qū)域為X=(0,1)，邊界條件在計算區(qū)域邊界處取內點的一階插值，初始間斷分布為：

(ul,ρl,pl)=(0,1,1) 0≤x<0.5

(ur,ρr,pr)=(0,0.125,0.1) 0.5≤x≤1

2.2 控制方程

不考慮熱源、熱傳導和體積力的影響，控制方程采用一維非定常歐拉方程的守恒形式：

式中:

(15)

補充理想氣體狀態(tài)方程使方程組封閉：

(16)

3 基于GPU的一維歐拉求解器

3.1 CPU/GPU并行模式

GPU不是一個獨立執(zhí)行的平臺，GPU的工作需要CPU協(xié)同，在異構計算體系中，CPU與GPU通過PCIE總線連接在一起協(xié)同工作。CPU稱為主機端(host)，GPU稱為設備端(device)。CPU可以實現(xiàn)復雜的邏輯運算，但其計算單元較少；與此相反，GPU包含大量的計算單元，相較于CPU，GPU線程是輕量級的，上下文切換開銷小，比較適合執(zhí)行大規(guī)模數(shù)據(jù)并行的計算密集型任務。

GPU的核心組件是流多處理器(SM,streaming multiprocessor)，SM的核心組件包括CUDA計算核心、共享內存、寄存器等，SM可以并發(fā)執(zhí)行數(shù)百上千線程，并發(fā)能力依賴SM擁有的資源數(shù)。GPU執(zhí)行模式為單指令多線程(SIMT,single instruction multiple thread)模式， 即一個warp(通常為32個線程)執(zhí)行同一條指令，遇見條件分支時，按分支順序執(zhí)行。

本文采用CUDA對GPU進行編程，CUDA 是NVIDIA公司推出的通用并行計算架構。一個典型的CUDA程序包括Host端和Device端兩部分代碼，Host代碼在CPU執(zhí)行，Device代碼調用內核在GPU上執(zhí)行。 典型CUDA程序的執(zhí)行流程如下：

1)分配host端內存，進行數(shù)據(jù)初始化；

2)分配device端內存，將host數(shù)據(jù)拷貝至device上；

3)調用CUDA核函數(shù)在device上執(zhí)行運算；

4)將device上計算結果拷貝至host；

5)釋放device端和host端內存。

3.2 基于CPU/GPU的CFD并行算法

圖2展示了一個典型的CFD算法流程，CFD程序主要包括了3個關鍵部分，時間步長計算，算法核心求解步計算和進行邊界條件處理。其中，算法的核心求解步占用了CFD計算的大量時間。在CPU/GPU異構并行計算體系如何分配這三大部分將會決定整個求解器的性能。

圖2 CFD算法流程圖

本文將時間步長計算分為兩部分，全局時間步長的計算和系統(tǒng)迭代時間步的計算。迭代時間步的計算涉及到復雜的邏輯判斷和循環(huán)過程。鑒于此，將迭代時間步的計算置于CPU上執(zhí)行。針對如何分配CFD程序三個關鍵部分的計算，本文考慮兩種并行算法。

算法一：對于算法核心步驟的計算，需要將流場變量u從CPU傳輸至GPU(主機端到設備端)；而邊界條件的計算，需要用到核心算法求解器計算過后的流場變量d_u(在GPU上分配的變量數(shù)組之前添加d，以區(qū)別CPU上具有相同功能的數(shù)組)，d_u需要從設備端傳輸至主機端。圖3給出了算法一的示意圖。在本模式中，整個計算過程只有核心求解器的計算置于GPU，其余部分全部交由CPU完成。

算法二：圖4給出了模式二的示意圖，全局時間步長的計算、核心算法求解器的求解和邊界條件的處理全部置于GPU上執(zhí)行，最小化CPU和GPU之間的數(shù)據(jù)傳輸。主機流場變量u和d_u分別僅在循環(huán)之前和之后傳輸一次，而只有聲速d_a需要從設備端單向傳輸至主機端進行迭代時間的計算。

圖3 并行算法一

圖4 并行算法二

對于GPU上的大多數(shù)應用程序，限制程序性能瓶頸的是CPU和GPU之間的帶寬[14-16]，而非GPU的浮點性能(單塊GPU RTX-3090的浮點性能已達到35.7 TFLOP/s)。此外，在許多內存綁定CFD程序中，帶寬達到極限，浮點性能遠遠小于GPU峰值的情況經(jīng)常出現(xiàn)。綜合考慮，本文采用算法二進行求解器實現(xiàn)。

4 算例分析

本文建立一維歐拉求解器，對Sod激波管問題進行計算，分析了通量分裂方法的典型差分格式在CPU/GPU異構并行體系中的計算性能。本文主要是針對格式本身性能進行研究，故Van Leer格式、Roe格式和AUSMPW+格式均采用一階計算，計算精度為雙精度。

4.1 GPU計算結果

測試采用的CPU為INTEL I7處理器，8核，主頻為3.2 GHz，主機內存32 GB，GPU為NVDIA GTX 1 660 super，擁有1 408個CUDA計算核心，單精度浮點性能約5 300 GFLOP/s，雙精度浮點性能約170 GFLOP/s。計算取時結果。

t=0.2時刻，激波管內流場結構如圖5所示。激波管內流場被分為A、B、C、D共4個區(qū)域，A區(qū)域和D區(qū)域中波未傳播，未受到擾動干擾，保持著流動初始狀態(tài)；B區(qū)域為膨脹波掃過后的區(qū)域；C區(qū)域為激波掃過區(qū)域。B區(qū)域和C區(qū)域以一個接觸間斷分開，B、C兩區(qū)域壓力相同、密度和溫度不同。

圖5 t=0.2波系結構圖

圖6 Van Leer計算結果

圖7 Roe計算結果

圖8 AUSMPW+計算結果

圖6～8分別展示了t=0.2時，3種格式在10 000個網(wǎng)格點下的計算結果。如圖所示，3種格式都能清晰地分辨出所求時刻流場的波系結構。氣體壓力差產(chǎn)生的激波從x=0.5處向右傳播，t=0.2時刻傳播至x=0.85附近，產(chǎn)生速度、密度和壓力的躍遷；x=0.25和x=0.5之間3個變量緩慢變化，形成膨脹波；x=0.68處，速度和壓力不變，密度產(chǎn)生躍遷，為接觸間斷。通過與精確值對比，3種格式在CPU/GPU異構體系下得出合理且可用的計算結果。

4.2 格式效應分析

圖9展示了t=0.2時刻，3種格式在256網(wǎng)格點下，流場的密度分布，為了便于比較格式對流場波系的分辨情況，給出t=0.2時刻的流場密度精確解進行比較。可以看出，本文所采用的3種格式在256網(wǎng)格下都能表述所求時刻間斷分布。格式本身對于間斷的捕捉有差別，圖10展示了激波附近局部放大圖， Roe 激波分辨率最高，AUSMPW+與Roe激波分辨率相當，Van Leer激波分辨率低于Roe和AUSMPW+。值得注意的是AUSMPW+格式，它是一種FVS和FDS的復合格式，兼有FVS的計算效率和FDS的間斷高分辨率。計算結果表明，它在間斷前后表現(xiàn)出與Roe格式相當?shù)男阅埽瑫r又兼有與Van Leer格式相當?shù)挠嬎阈省?/p>

圖9 3種格式計算對比

圖10 激波局部放大圖

4.3 格式加速性能分析

本文研究格式本身在CPU/GPU上性能分析，加速比定義為算法核心求解步在CPU運行時間與相應GPU時間的比值，為排除額外因素干擾，時間步長取固定值0.000 001。

表1為3種典型格式在CPU執(zhí)行時間(100、1 000、10 000共3種網(wǎng)格計算規(guī)模)，由計算結果可知，Van Leer計算效率最高，AUSMPW+計算所需時間略高于與Van Leer，Roe消耗時間最多。表2展示了3種格式在GPU執(zhí)行時間。在本文的計算的網(wǎng)格規(guī)模中，GPU硬件體系結構上，Van Leer計算效率稍高于Roe和AUSMPW+。表3體現(xiàn)了3種格式的加速比，在網(wǎng)格規(guī)模較小時，GPU計算時間大于CPU計算。主要原因是網(wǎng)格規(guī)模小，無法體現(xiàn)GPU眾核體系進行大規(guī)模計算的優(yōu)勢。加速比由高到低依次為Roe、AUSMPW+、Van Leer。其影響因素有：Roe格式構造中，需要對矩陣進行計算，其計算量遠大于Van Leer和AUSMPW+，而且Roe沒有條件語句分支，控制流指令較少，運行中warp分支少，格式適用于GPU單指令多線程運行模式，故Roe加速效果最好；AUSMPW+和Van Leer兩種格式存在大量條件分支，不利于GPU運行模式，影響了加速效果；AUSMPW+加速優(yōu)于Van Leer的原因是AUSMPW+計算量大于Van Leer(從格式構造來看，AUSMPW+是Van Leer的改進，其計算步驟更多更復雜)。從以上分析結果來看，影響格式在異構平臺的加速性能主要是格式存在的條件判斷分支以及格式的運算量。

表1 CPU不同網(wǎng)格規(guī)模格式計算時間 s

表2 GPU不同網(wǎng)格規(guī)模格式計算時間 s

表3 格式加速比 s

5 結束語

為了研究通量分裂格式在異構并行計算平臺的加速性能，本文在 GTX 1660 super上利用CUDA開發(fā)了一維歐拉求解器，并對當前航空航天領域常用的3種典型格式進行計算分析，結果表明：

1)Roe格式對流場解析效果最好，且在異構平臺加速性能最高，目前應用前景最好。

2)格式構造中存在的條件分支嚴重影響了格式的加速性能，構造數(shù)值格式時，合理地減少條件分支，將會大大提高其在異構平臺的加速效果。

本文的工作對于在CPU/GPU異構并行計算平臺開展CFD算法相關研究具有一定的指導意義。