基于異構(gòu)平臺的通量分裂格式性能研究

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

(西南科技大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，四川 綿陽 621010)

0 引言

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

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

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

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

1 通量分裂方法

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

1.1 矢通量分裂

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

一維歐拉方程進(jìn)行矢通量分裂：

(1)

Van Leer格式通量分裂表達(dá)式：

Ma>1:

F+=0,F-=F

(2)

|Ma|>1:

(3)

Ma<-1:

F+=F,F-=0

(4)

1.2 通量差分分裂

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

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

(5)

(6)

Roe格式通量表達(dá)式為：

(7)

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

1.3 AUSM分裂方法

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

AUSMPW+格式具體表達(dá)式：

(8)

(9)

(10)

(11)

(12)

(13)

(14)

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

2 Riemann問題與計(jì)算模型

2.1 Riemann問題

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

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

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

圖1 Sod激波管問題

對于本次計(jì)算問題，計(jì)算區(qū)域?yàn)閄=(0,1)，邊界條件在計(jì)算區(qū)域邊界處取內(nèi)點(diǎn)的一階插值，初始間斷分布為：

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

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

2.2 控制方程

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

式中:

(15)

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

(16)

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

3.1 CPU/GPU并行模式

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

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

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

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

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

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

4)將device上計(jì)算結(jié)果拷貝至host；

5)釋放device端和host端內(nèi)存。

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

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

圖2 CFD算法流程圖

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

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

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

圖3 并行算法一

圖4 并行算法二

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

4 算例分析

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

4.1 GPU計(jì)算結(jié)果

測試采用的CPU為INTEL I7處理器，8核，主頻為3.2 GHz，主機(jī)內(nèi)存32 GB，GPU為NVDIA GTX 1 660 super，擁有1 408個(gè)CUDA計(jì)算核心，單精度浮點(diǎn)性能約5 300 GFLOP/s，雙精度浮點(diǎn)性能約170 GFLOP/s。計(jì)算取時(shí)結(jié)果。

t=0.2時(shí)刻，激波管內(nèi)流場結(jié)構(gòu)如圖5所示。激波管內(nèi)流場被分為A、B、C、D共4個(gè)區(qū)域，A區(qū)域和D區(qū)域中波未傳播，未受到擾動干擾，保持著流動初始狀態(tài)；B區(qū)域?yàn)榕蛎洸⊕哌^后的區(qū)域；C區(qū)域?yàn)榧げ⊕哌^區(qū)域。B區(qū)域和C區(qū)域以一個(gè)接觸間斷分開，B、C兩區(qū)域壓力相同、密度和溫度不同。

圖5 t=0.2波系結(jié)構(gòu)圖

圖6 Van Leer計(jì)算結(jié)果

圖7 Roe計(jì)算結(jié)果

圖8 AUSMPW+計(jì)算結(jié)果

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

4.2 格式效應(yīng)分析

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

圖9 3種格式計(jì)算對比

圖10 激波局部放大圖

4.3 格式加速性能分析

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

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

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

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

表3 格式加速比 s

5 結(jié)束語

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

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

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

本文的工作對于在CPU/GPU異構(gòu)并行計(jì)算平臺開展CFD算法相關(guān)研究具有一定的指導(dǎo)意義。