應(yīng)用多GPU的可壓縮湍流并行計(jì)算*

曹文斌，李 樺，謝文佳，張 冉

(國(guó)防科技大學(xué)航天科學(xué)與工程學(xué)院，湖南長(zhǎng)沙410073)

隨著硬件性能的提高及編程技術(shù)的改進(jìn)，圖形處理器(Graphical Processing Unit，GPU)加速器在高性能計(jì)算領(lǐng)域逐漸得到廣泛的應(yīng)用。在最新公布的超級(jí)計(jì)算機(jī)Top500名單中共有62套系統(tǒng)采用了加速器/協(xié)處理器，其中采用GPU加速器有46套，而在最新的Green500名單中前10位的超級(jí)計(jì)算機(jī)均采用了GPU加速器，由此可見(jiàn)GPU加速器在高性能計(jì)算中的優(yōu)勢(shì)。在通用計(jì)算中NVIDA GPU具有低成本、低功耗、高性能及易編程等優(yōu)點(diǎn)，非常適合計(jì)算密集型應(yīng)用。發(fā)布的CUDA并行編程架構(gòu)［1］使開(kāi)發(fā)人員能夠在不了解圖形學(xué)知識(shí)的條件下進(jìn)行高級(jí)語(yǔ)言環(huán)境下的GPU編程以高效快速地解決許多復(fù)雜計(jì)算問(wèn)題。目前在流體力學(xué)領(lǐng)域，已有大量的數(shù)值模擬方法針對(duì)GPU進(jìn)行了改造，這其中有計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)、格子玻爾茲曼、直接模擬蒙特卡洛、分子模擬等方法。

湍流一直是CFD領(lǐng)域的難點(diǎn)與熱點(diǎn)，近年來(lái)，國(guó)內(nèi)外已有許多學(xué)者將GPU計(jì)算應(yīng)用到湍流的數(shù)值模擬當(dāng)中，取得了良好的加速效果。Aaron［2］針對(duì)不可壓縮湍流問(wèn)題，基于Tesla C1060 GPU發(fā)展了大渦模擬求解器并給出了實(shí)現(xiàn)多GPU并行計(jì)算的方法。Cao與Xu等［3］應(yīng)用高精度求解器在天河－1A系統(tǒng)上實(shí)現(xiàn)了多GPU與多CPU的協(xié)同并行計(jì)算，討論了提高異構(gòu)系統(tǒng)并行效率的方法。Lin等［4］基于GeForce 560Ti GPU發(fā)展了多塊結(jié)構(gòu)網(wǎng)格湍流求解器，利用MPI加OpenMP實(shí)現(xiàn)了多GPU的并行計(jì)算。Ali等［5］在Tesla S1070多GPU集群上實(shí)現(xiàn)了湍流的直接數(shù)值模擬，分析了利用零拷貝技術(shù)與固定內(nèi)存技術(shù)對(duì)提高系統(tǒng)效率的影響?？梢园l(fā)現(xiàn)，在現(xiàn)有的計(jì)算流體力學(xué)應(yīng)用中，大多數(shù)研究針對(duì)的是計(jì)算能力較低的設(shè)備。而最新的設(shè)備(計(jì)算能力3.5)更新了硬件架構(gòu)，例如引入了Hyper－Q，提升了顯存帶寬，增加了可用的寄存器個(gè)數(shù)等，這就需要研究人員根據(jù)最新架構(gòu)的功能與特點(diǎn)對(duì)算法做出相應(yīng)的調(diào)整與改進(jìn)，以發(fā)揮設(shè)備的最佳性能，實(shí)現(xiàn)更加復(fù)雜更大計(jì)算規(guī)模任務(wù)的模擬。研究發(fā)展了基于CUDA Fortran的三維定?？蓧嚎s湍流求解器。

1 控制方程

當(dāng)不考慮外部加熱和徹體力的影響時(shí)，三維曲線(xiàn)坐標(biāo)系下守恒形式的Navier-Stokes方程如下:

此處Q為守恒變量;J代表曲線(xiàn)坐標(biāo)系(ξ，η，ζ)與笛卡爾坐標(biāo)系(x，y，z)之間的雅克比轉(zhuǎn)換行列式;F，G，H為曲線(xiàn)坐標(biāo)系三個(gè)方向的無(wú)黏通量;Fv，Gv，Hv代表黏性通量?？紤]k－ωSST湍流模型時(shí)，除了添加兩個(gè)包含源項(xiàng)的湍流方程之外，方程(1)的形式不變，此時(shí)黏性系數(shù)μ由層流黏性系數(shù)μL及渦黏性系數(shù)μT兩部分組成，相應(yīng)的μ/pr關(guān)系式如下:

式中，μL通過(guò)Sutherland公式得到，μT由湍流模型給出，prL為層流普朗特?cái)?shù)，對(duì)于空氣可取0.72，prT為湍流普朗特?cái)?shù)，一般取0.9。當(dāng)不考慮湍流模型時(shí)，渦黏性系數(shù)為0。湍流模型的具體公式及參數(shù)設(shè)置見(jiàn)文獻(xiàn)［6］。

對(duì)控制方程(1)運(yùn)用基于單元中心型的有限體積法求解，無(wú)黏通量采用AUSMPW+格式［7］，通過(guò)MUSCL方法選用Vanleer平均限制器對(duì)原始變量進(jìn)行空間重構(gòu)達(dá)到二階精度，黏性通量采用高斯積分公式進(jìn)行求解，時(shí)間推進(jìn)為單步顯式方法。

2 GPU算法優(yōu)化方法

由于GPU程序優(yōu)化與硬件的體系結(jié)構(gòu)及軟件環(huán)境密切相關(guān)，因此本節(jié)首先給出計(jì)算平臺(tái)配置情況。計(jì)算軟件環(huán)境如下:CPU代碼用Fortran語(yǔ)言編寫(xiě)，GPU代碼為CUDA Fortran語(yǔ)言;兩個(gè)版本的代碼均由PGI Fortran 13.10編譯器［8］加OpenMPI庫(kù)編譯，優(yōu)化參數(shù)設(shè)置使兩種設(shè)備的性能達(dá)到最佳;GPU計(jì)算使用單精度。硬件配置如下:每個(gè)節(jié)點(diǎn)配置一塊CPU與一塊GPU及64G內(nèi)存，節(jié)點(diǎn)間通信采用鏈路匯聚的千兆以太網(wǎng)絡(luò);運(yùn)行的CPU為Intel Xeon E5－2670@2.6 GHz，超線(xiàn)程關(guān)閉;測(cè)試使用的GPU為6G顯存的NVIDA GTX TITAN Black。

2.1 單GPU算法優(yōu)化

CFD迭代求解過(guò)程主要包括邊界條件更新、通量計(jì)算、時(shí)間步長(zhǎng)與源項(xiàng)計(jì)算及流場(chǎng)更新四個(gè)步驟。其中通量計(jì)算是決定求解器整體運(yùn)行速度的主要因素，通常是優(yōu)化的重點(diǎn)。不管是無(wú)黏通量還是黏性通量，在GPU上它們的計(jì)算都是首先從全局內(nèi)存中輸入網(wǎng)格單元的幾何量和原始變量;其次通過(guò)空間格式得到單元界面通量;最后將通量或通量差結(jié)果輸出至全局內(nèi)存中。因此從優(yōu)化方面考慮，可以將通量計(jì)算方法分為兩類(lèi)算法:輸入優(yōu)化(算法1)與輸出優(yōu)化(算法2)。為了分析對(duì)比，將不進(jìn)行優(yōu)化的算法定義為全局內(nèi)存方法。

算法1的具體流程如下所示。為了節(jié)省全局內(nèi)存開(kāi)銷(xiāo)，界面通量數(shù)組f被多次重復(fù)使用。

算法1輸入優(yōu)化Alg.1 Input optimization

算法1中每個(gè)內(nèi)核在讀取數(shù)據(jù)階段，均使用了下面介紹的優(yōu)化方法。由于每個(gè)線(xiàn)程都需要用到相鄰線(xiàn)程的信息，例如在計(jì)算ξ方向無(wú)黏通量Fi－1/2時(shí)涉及i－2，i－1，i，i+1共4個(gè)網(wǎng)格單元的原始變量值，若應(yīng)用全局內(nèi)存方法(將數(shù)據(jù)從全局內(nèi)存直接讀入寄存器)讀取這些信息，將會(huì)產(chǎn)生大量的重復(fù)訪(fǎng)問(wèn)。為了提高帶寬利用率降低訪(fǎng)問(wèn)延遲，常用的辦法就是使用共享內(nèi)存來(lái)緩解對(duì)全局內(nèi)存的重復(fù)訪(fǎng)問(wèn)，具體使用方法見(jiàn)文獻(xiàn)［9］。文獻(xiàn)［10］給出了無(wú)黏通量計(jì)算的一種更高效的方法，即使用寄存器移位技術(shù)達(dá)到減少重復(fù)訪(fǎng)問(wèn)的目的，同時(shí)還給出了黏性通量的優(yōu)化方法。本文求解器在算法1中采用文獻(xiàn)［10］的方法進(jìn)行優(yōu)化。圖1給了大小為1923的網(wǎng)格在不同GPU架構(gòu)下的通量計(jì)算時(shí)間對(duì)比，與全局內(nèi)存方法相比，算法1在計(jì)算能力2.0的Tesla C2050上獲得了25.1%的性能提升，而在計(jì)算能力3.5的GTX TITAN Black上僅有8.3%的提升?？梢钥闯觯瑢?duì)于重復(fù)數(shù)據(jù)的讀取，在顯存帶寬有了顯著提升并引入了L2緩存的計(jì)算能力3.5的設(shè)備上，使用優(yōu)化加速技術(shù)并不能取得明顯的加速效果。

算法2對(duì)結(jié)果輸出做優(yōu)化，具體的流程如下

算法2輸出優(yōu)化Alg.2 Output optimization

可以看出，算法2在輸出階段利用共享內(nèi)存消除了中間變量f的使用，這樣既節(jié)省了全局內(nèi)存消耗又避免了對(duì)全局變量的額外讀寫(xiě)操作，同時(shí)所有內(nèi)核不存在數(shù)據(jù)依賴(lài)性可并發(fā)執(zhí)行。因?yàn)榇鎯?chǔ)于網(wǎng)格中心的右端項(xiàng)RHS是通過(guò)兩個(gè)界面通量做差值得到，所以每個(gè)線(xiàn)程塊中計(jì)算得到的RHS個(gè)數(shù)少于線(xiàn)程塊中的線(xiàn)程數(shù)，這樣在每個(gè)線(xiàn)程塊中除了需要引入新的分支語(yǔ)句，還需要在相鄰線(xiàn)程塊的邊界處進(jìn)行冗余計(jì)算。分支語(yǔ)句的增多與冗余計(jì)算的存在使得算法1中的輸入優(yōu)化技術(shù)在算法2中不但不能發(fā)揮應(yīng)有的加速效果而且會(huì)引起性能的下降。因此算法2在讀入數(shù)據(jù)階段不進(jìn)行優(yōu)化，即應(yīng)用全局內(nèi)存方法讀取數(shù)據(jù)。由圖1可以看出，算法2非常適合計(jì)算能力3.5的GTX TITAN Black，與全局內(nèi)存方法相比，獲得了18%的速度提升。

除了上述優(yōu)化方法外，一些常規(guī)的優(yōu)化方法在本文中均有應(yīng)用，如全局變量按數(shù)組結(jié)構(gòu)體方式儲(chǔ)存，盡量避免束內(nèi)分支，盡可能減少中間變量以降低寄存器使用量，減少PCIe數(shù)據(jù)傳輸，采用樹(shù)形規(guī)約算法執(zhí)行殘值加和操作。所有計(jì)算配置的線(xiàn)程塊大小為32×8×1。

圖1 不同算法計(jì)算時(shí)間對(duì)比Fig.1 Computing time of different algorithms

2.2 多GPU并行算法

在求解器進(jìn)行基于區(qū)域分解的并行計(jì)算過(guò)程中，相鄰網(wǎng)格區(qū)域在邊界處需要對(duì)原始變量進(jìn)行數(shù)據(jù)交換。當(dāng)多區(qū)網(wǎng)格位于同一塊GPU中，此時(shí)為內(nèi)部邊界，可以在GPU中實(shí)現(xiàn)直接數(shù)據(jù)交換(求解器支持一個(gè)MPI進(jìn)程處理多區(qū)結(jié)構(gòu)網(wǎng)格);當(dāng)相鄰網(wǎng)格區(qū)域位于不同GPU時(shí)，需要先將一塊GPU的邊界數(shù)據(jù)通過(guò)PCIe總線(xiàn)傳輸至主機(jī)端，然后調(diào)用MPI進(jìn)行交換，最后再通過(guò)PCIe傳輸至另一塊GPU。雖然已有GPUDirect技術(shù)實(shí)現(xiàn)GPU間的直接通信，但是現(xiàn)有系統(tǒng)的通信設(shè)備并不具有該功能，而升級(jí)設(shè)備成本又很高，因此提高現(xiàn)有系統(tǒng)的并行計(jì)算效率具有很大的實(shí)用價(jià)值。計(jì)算能力3.5的GPU引入了一個(gè)名為Hyper－Q的重要功能，它允許32個(gè)獨(dú)立的CUDA流隊(duì)列同時(shí)工作。利用該功能，可以實(shí)現(xiàn)GPU計(jì)算與邊界數(shù)據(jù)交換的高度重疊。具體實(shí)現(xiàn)過(guò)程如圖2所示。

首先，創(chuàng)建N個(gè)CUDA流(N≤32)。流1與流2分別用于所有區(qū)域的ζ方向內(nèi)點(diǎn)無(wú)黏通量計(jì)算與內(nèi)部邊界更新等內(nèi)核的啟動(dòng)。流3至流N負(fù)責(zé)所有需要進(jìn)行數(shù)據(jù)交換的邊界處理，這其中包括GPU上數(shù)據(jù)整理(將邊界數(shù)據(jù)整理至預(yù)先分配好的連續(xù)儲(chǔ)存空間中)，設(shè)備端至主機(jī)端的異步傳輸(Device to Host，D2H)，主機(jī)端至設(shè)備端的異步傳輸(Host to Device，H2D)及GPU邊界數(shù)據(jù)更新。其次，所有流啟動(dòng)相應(yīng)的內(nèi)核執(zhí)行，流3至流N的每個(gè)流在啟動(dòng)MPI非阻塞發(fā)送之前應(yīng)用CUDA流查詢(xún)函數(shù)確保D2H傳輸已完成，在啟動(dòng)H2D傳輸之前需等待MPI非阻塞接受完成。由于Hyper-Q特性允許流與流之間的操作互不阻塞，所以流1與流2對(duì)應(yīng)的多個(gè)內(nèi)核從一開(kāi)始就與數(shù)據(jù)整理內(nèi)核同時(shí)啟動(dòng)，內(nèi)核的執(zhí)行與PCIe傳輸及MPI數(shù)據(jù)交換同步進(jìn)行，多個(gè)D2H或H2D傳輸可以同時(shí)存在，MPI非阻塞通信與D2H傳輸在不同流之間重疊。最后調(diào)用流同步完成所有流的內(nèi)核執(zhí)行。圖中流1雖然只包含了ζ方向的內(nèi)點(diǎn)無(wú)黏通量計(jì)算，若它的計(jì)算時(shí)間小于數(shù)據(jù)交換所需要的時(shí)間，可以將η方向的相應(yīng)計(jì)算加入流隊(duì)列中，以實(shí)現(xiàn)GPU計(jì)算與邊界數(shù)據(jù)交換的高度重疊。

圖2 多GPU并行算法Fig.2 Multi-GPU parallel algorithm

3 算例分析

3.1 超聲速進(jìn)氣道

Reinartz等［11］對(duì)混合壓縮進(jìn)氣道進(jìn)行了大量的數(shù)值與實(shí)驗(yàn)研究，獲得了清晰的流場(chǎng)結(jié)構(gòu)與進(jìn)氣道壁面壓力分布。圖3所示為擴(kuò)張段以前的進(jìn)氣道構(gòu)型，上壓縮面的壓縮角δ2為21.5°，下楔面的傾斜角δ3為9.5°，擴(kuò)張段的擴(kuò)張角δ4為5°，隔離段的長(zhǎng)度l為79.3mm，高度h為15mm，進(jìn)氣道總長(zhǎng)為400mm，寬為52mm。

圖3 進(jìn)氣道模型示意圖Fig.3 Schematic of the inlet model

進(jìn)氣道流動(dòng)計(jì)算條件為:攻角α為－10°，馬赫數(shù)M∞為2.41，總壓540kPa，總溫305K。由于實(shí)驗(yàn)時(shí)間較長(zhǎng)，可認(rèn)為壁面已達(dá)到熱傳導(dǎo)平衡狀態(tài)，因此計(jì)算時(shí)壁面采用絕熱條件假設(shè)。采用單區(qū)結(jié)構(gòu)網(wǎng)格對(duì)進(jìn)氣道流場(chǎng)區(qū)域進(jìn)行離散。為了測(cè)定GPU在不同網(wǎng)格量下的性能表現(xiàn)，設(shè)置了多套網(wǎng)格，其中最粗的網(wǎng)格由256×48×64(流向×法向×展向)個(gè)網(wǎng)格點(diǎn)組成，在此基礎(chǔ)上將前兩個(gè)方向的網(wǎng)格點(diǎn)數(shù)量分別加密得到其余的網(wǎng)格，最密的網(wǎng)格其網(wǎng)格點(diǎn)個(gè)數(shù)達(dá)5033萬(wàn)(2048×384×64)。各套網(wǎng)格均在壁面處進(jìn)行了加密。

圖4 計(jì)算得到的密度梯度圖Fig.4 Density gradient magnitude of the computation

圖5 下表面壓力分布Fig.5 Lower surface pressure distribution

圖4給出了最密網(wǎng)格下湍流模型計(jì)算所得到的流場(chǎng)對(duì)稱(chēng)面密度梯度云圖，可以發(fā)現(xiàn)數(shù)值計(jì)算結(jié)果清晰地捕獲了進(jìn)氣道內(nèi)復(fù)雜的波系結(jié)構(gòu)。圖5給出了進(jìn)氣道下壁面無(wú)量綱壓力分布計(jì)算與實(shí)驗(yàn)結(jié)果的對(duì)比。從圖中結(jié)果可以看出湍流模型的計(jì)算值與實(shí)驗(yàn)值在整體分布規(guī)律上都比較符合，這表明本文所發(fā)展的GPU求解器對(duì)可壓縮湍流定常問(wèn)題的模擬是可行的。

圖6給出了不同網(wǎng)格下單塊GPU與單個(gè)CPU核心的計(jì)算時(shí)間對(duì)比及相應(yīng)的加速比。圖中所示的時(shí)間為每迭代步所消耗的平均時(shí)間，可以看出GPU計(jì)算具有明顯的速度優(yōu)勢(shì)，相對(duì)于單個(gè)CPU核心而言，隨著網(wǎng)格量的增加，GPU的加速比由最初的107倍逐漸增大至125倍。由加速比曲線(xiàn)可發(fā)現(xiàn)，小網(wǎng)格量下的加速比明顯低于大網(wǎng)格量的情況，即GPU非常適合處理大網(wǎng)格量的計(jì)算任務(wù)。同時(shí)，隨著網(wǎng)格量的增加，加速比趨于平緩，這是由于GPU中用于并行計(jì)算的CUDA核心數(shù)及可并發(fā)執(zhí)行的線(xiàn)程塊數(shù)是有限的，當(dāng)CUDA核心的利用率達(dá)到飽和時(shí)一些線(xiàn)程將會(huì)采用串行方式執(zhí)行，此時(shí)加速比將不再增加而趨于穩(wěn)定。對(duì)于不同架構(gòu)的設(shè)備，這個(gè)飽和值是不同的，可以看出，當(dāng)前使用的GPU負(fù)載飽和值為500萬(wàn)以上的網(wǎng)格量。

圖6 不同網(wǎng)格的計(jì)算時(shí)間與加速比Fig.6 Computing time and speedup for different grid sizes

并行效率是并行計(jì)算必須關(guān)注的性能指標(biāo)。受顯存容量的限制，單塊GPU所能處理的最大網(wǎng)格量是有限的(對(duì)于定常湍流問(wèn)題，6GB顯存至多可處理5000萬(wàn)網(wǎng)格)，考慮到網(wǎng)格規(guī)模大小對(duì)GPU的性能發(fā)揮影響較大，為了真實(shí)反映同步與通信開(kāi)銷(xiāo)對(duì)并行效率的影響，本文定義并行效率如下:對(duì)于某個(gè)計(jì)算任務(wù)，假設(shè)它在多塊GPU上并行處理的時(shí)間為tp，而對(duì)每個(gè)GPU來(lái)說(shuō)都負(fù)載一定數(shù)量的網(wǎng)格(單個(gè)或多個(gè)分區(qū))，設(shè)所有GPU單獨(dú)處理相應(yīng)的負(fù)載所需時(shí)間的平均值為ts，則并行效率為ts與tp之比。上述的GPU單獨(dú)處理既不考慮節(jié)點(diǎn)間通信又不考慮節(jié)點(diǎn)內(nèi)通信。對(duì)于本算例，每塊GPU的負(fù)載網(wǎng)格量都相等，且網(wǎng)格拓?fù)浣Y(jié)構(gòu)簡(jiǎn)單無(wú)節(jié)點(diǎn)內(nèi)通信。表1給出了4塊GPU在不同網(wǎng)格量下的計(jì)算時(shí)間與并行效率，可以看出，較小網(wǎng)格量下GPU的并行效率不高，隨著網(wǎng)格量的增加GPU的并行效率逐漸提高。這是因?yàn)殡S著網(wǎng)格量的增加，節(jié)點(diǎn)間的通信時(shí)間在總的時(shí)間中所占比重逐漸減小所以并行效率得到提高。由表中數(shù)據(jù)可知，應(yīng)用本文所發(fā)展的計(jì)算與通信的重疊算法，可使并行效率得到較大提高，相對(duì)于無(wú)重疊的并行計(jì)算，四套網(wǎng)格下內(nèi)點(diǎn)通量計(jì)算與數(shù)據(jù)傳輸及通信的重疊可使并行效率提高10%以上。

表1 計(jì)算時(shí)間與并行效率Tab.1 Computing time and parallel efficiency

3.2 超聲速進(jìn)氣道

針對(duì)文獻(xiàn)［12］所述的空天飛機(jī)模型，本文計(jì)算了馬赫數(shù)為4.94，攻角為－5°的實(shí)驗(yàn)狀態(tài)，對(duì)應(yīng)的雷諾數(shù)為5.26×107m?？仗祜w機(jī)的長(zhǎng)度為0.29m，文獻(xiàn)提供了高精度的實(shí)驗(yàn)結(jié)果。計(jì)算方法如前所述，按定常湍流狀態(tài)考慮，壁面采用絕熱條件假設(shè)。

圖7所示為空天飛機(jī)的表面及對(duì)稱(chēng)面網(wǎng)格分布，網(wǎng)格總量為1.34億，共有26個(gè)分區(qū)，壁面第一層網(wǎng)格高度為1×10－6m。采用4塊GPU并行計(jì)算，每塊GPU負(fù)載網(wǎng)格量為3350萬(wàn)左右。圖8給出空天飛機(jī)的壓力及流線(xiàn)分布，因外形復(fù)雜，流場(chǎng)中不可避免地會(huì)出現(xiàn)激波與激波干擾及流動(dòng)分離等現(xiàn)象。

圖7 空天飛機(jī)計(jì)算網(wǎng)格Fig.7 Computational grid for the space shuttle

圖9給出了機(jī)身對(duì)稱(chēng)面中心線(xiàn)上表面的壓力分布，可以看出，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好，能夠反映壓力的變化趨勢(shì)，這說(shuō)明了求解器對(duì)復(fù)雜外形模擬具有很好的適應(yīng)性。

4塊GPU并行計(jì)算時(shí)，每迭代步所消耗的時(shí)間為0.667秒，每塊GPU單獨(dú)處理的平均時(shí)間為0.611秒，計(jì)算得到的并行效率為91.6%，與超聲速進(jìn)氣道算例的最高并行效率相比有所降低。這主要是因?yàn)閱螇KGPU負(fù)載的網(wǎng)格分區(qū)過(guò)多(尾流場(chǎng)包含18個(gè)分區(qū)，分配在兩個(gè)GPU上)引起同步時(shí)間消耗增多。此外，負(fù)載不平衡也造成了一定的影響，雖然4個(gè)GPU所負(fù)載的網(wǎng)格量基本相等，但是內(nèi)部邊界及物理邊界條件(如壁面條件、遠(yuǎn)場(chǎng)條件等)更新個(gè)數(shù)并不相等，分區(qū)過(guò)多導(dǎo)致了節(jié)點(diǎn)內(nèi)的邊界更新操作過(guò)多從而影響了整體并行效率。

圖8 空天飛機(jī)壓力云圖及空間流線(xiàn)分布Fig.8 Pressure contours and stream lines distribution of the space shuttle

圖9 空天飛機(jī)對(duì)稱(chēng)面上表面的壁面壓力分布Fig.9 Pressure coefficient distribution along the upper surface in the symmetry plane

4 結(jié)論

本文建立的基于GPU的可壓縮湍流求解器，可以快速準(zhǔn)確地對(duì)超聲速流動(dòng)中的湍流問(wèn)題進(jìn)行求解。為了達(dá)到良好的加速效果，根據(jù)最新的GPU架構(gòu)特點(diǎn)對(duì)GPU程序的算法結(jié)構(gòu)進(jìn)行了調(diào)整與優(yōu)化，采用兩個(gè)超聲速流動(dòng)算例驗(yàn)證了GPU求解器的準(zhǔn)確性與適應(yīng)性，實(shí)現(xiàn)了多GPU上復(fù)雜外形上億量級(jí)網(wǎng)格的快速計(jì)算，在最新的GTX TITAN Black GPU上取得了可觀的加速比與并行效率。根據(jù)研究結(jié)果，得出以下結(jié)論:

1)不同架構(gòu)的GPU對(duì)應(yīng)不同的優(yōu)化算法，對(duì)于重復(fù)數(shù)據(jù)的讀取，在計(jì)算能力3.5的設(shè)備上，使用共享內(nèi)存或寄存器移位等優(yōu)化加速技術(shù)并不能取得明顯的加速效果。

2)利用Hyper－Q特性，可以實(shí)現(xiàn)GPU計(jì)算與PCIe數(shù)據(jù)傳輸、MPI通信的高度重疊，能有效掩蓋邊界交換所帶來(lái)的時(shí)間延遲。

3)從GPU的性能發(fā)揮及并行效率方面考慮，GPU負(fù)載的網(wǎng)格量應(yīng)越大越好。對(duì)于本文使用的GTX TITAN Black GPU，最佳的負(fù)載網(wǎng)格量應(yīng)在500萬(wàn)以上。

References)

［1］NVIDIA Corporation.NVIDIA CUDA C Programming Guide Version6.5［EB/OL］.［2014－08－01］.http://docs.nvidia.com/cuda/pdf/CUDA-C-Programming-Guide.pdf.

［2］Shinn A F.Large eddy simulations of turbulent flows on graphics processing units:application to film cooling flows［D］.USA:University of Illinois at Urbana Champaign，2011.

［3］Cao W，Xu C F，Wang Z H，et al.CPU/GPU computing for a multi-block structured grid based high-order flow solver on a large heterogeneous system［J］.CLUSTER COMPUT，2014，17(2):255－270.

［4］Fu L，Gao Z H，Xu K，et al.A multi-block viscous flow solver based on GPU parallel methodology［J］.Computers＆Fluids，2014，95:19－39.

［5］Khajeh-Saeed A，Perot J B.Direct numerical simulation of turbulence using GPU accelerated supercomputers［J］.Journal of Computational Physics，2013，235:241－257.

［6］Menter F R.Influence of freestream values on k－ω turbulence model prediction［J］.AIAA Journal，1993，30(6):1657－1659.

［7］Kim K H，Kim C，Rho O H.Methods for the accurate computations of hypersonic flows I.AUSMPW+scheme［J］.Journal of Computational Physics，2001，174(1):38－80.

［8］The Portland Group.CUDA Fortran Programming Guide and Reference Release 2013，version 13.10［EB/OL］.［2013］.http://www.pgroup.com/resources/docs.htm.

［9］Ruetsch G，F(xiàn)atica M.CUDA fortran for scientists and engineers:best practices for efficient CUDA fortran programming［M］.Holland:Elsevier，2013.

［10］Salvadore F，Bernardini M，Botti M.GPU accelerated flow solver for direct numerical simulation of turbulent flows［J］.Journal of Computational Physics，2013，235:129－142.

［11］Reinartz B U，Herrmann C D，Ballmann J，et al.Aerodynamic performance analysis of a hypersonic inlet isolator using computation and experiment［J］.Journal of Propulsion and Power，2003，19(5):868－875.

［12］李素循.典型外形高超聲速流動(dòng)特性［M］.北京:國(guó)防工業(yè)出版社，2007.LI Suxun.Hypersonic flow characteristics around typical con Figuration［M］.Beijing:National Defense Industry Press，2007.(in Chinese)