GPU加速的高精度可壓縮流動(dòng)計(jì)算

1）國(guó)家自然科學(xué)基金資助項(xiàng)目（項(xiàng)目編號(hào)：11272325）2）聯(lián)系作者Email：pengjun@第八屆全國(guó)流體力學(xué)學(xué)術(shù)會(huì)議2014年9月18～21日甘肅蘭州文章編號(hào)：CSTAM2014-B01-0275標(biāo)題：GPU加速的高精度可壓縮流動(dòng)計(jì)算作者：彭峻，申義慶單位：中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室Copyright?2014版權(quán)所有中國(guó)力學(xué)學(xué)會(huì)地址:北京市北四環(huán)西路15號(hào)郵政編碼:100190Address:No.15BeisihuanxiRoad,Beijing1001901）國(guó)家自然科學(xué)基金資助項(xiàng)目（項(xiàng)目編號(hào)：11272325）2）聯(lián)系作者Email：pengjun@第八屆全國(guó)流體力學(xué)學(xué)術(shù)會(huì)議2014年9月18-21日甘肅蘭州 CSTAM2014-A26-BS09057GPU加速的高精度可壓縮流動(dòng)計(jì)算1）彭峻*，2），申義慶**（中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京海淀區(qū)100190）摘要本文在GPU上實(shí)現(xiàn)了高精度可壓縮流動(dòng)計(jì)算的加速，獲得了可觀的加速效果。在可壓縮流動(dòng)計(jì)算中，激波的出現(xiàn)要求所使用的格式具有捕捉激波的能力。加權(quán)緊致格式是一種很好的激波捕捉格式。相比于其他激波捕捉格式，加權(quán)緊致格式具有高分辨率和低耗散的特點(diǎn)。然而，由于需要求解多對(duì)角方程組，加權(quán)緊致格式難于實(shí)現(xiàn)并行計(jì)算。本文針對(duì)三對(duì)角型加權(quán)緊致格式，采用Tomas方法，實(shí)現(xiàn)了在GPU上的并行求解，并將所實(shí)現(xiàn)的GPU程序應(yīng)用于可壓縮流動(dòng)計(jì)算中。GPU程序相比于相同網(wǎng)格數(shù)下的CPU程序獲得了18倍的加速效果。同GPU版本的WENO算法相比，本文所實(shí)現(xiàn)CRWENO格式的方法增加了一倍的計(jì)算時(shí)間，這與文獻(xiàn)中CPU程序的規(guī)律接近。關(guān)鍵詞GPU計(jì)算加權(quán)緊致格式可壓縮流動(dòng)引言GPU（GraphicsProcessingUnits）俗稱(chēng)顯卡,是計(jì)算機(jī)中用于渲染顯示畫(huà)面的設(shè)備。由于圖像的高度動(dòng)態(tài)性，顯卡需要以極快的速度完成圖像計(jì)算并進(jìn)行顯示。為此，顯卡芯片結(jié)構(gòu)被設(shè)計(jì)成具有高度并行計(jì)算能力，能夠承擔(dān)高強(qiáng)度計(jì)算。由于顯卡的這種特點(diǎn)，人們嘗試用顯卡進(jìn)行科學(xué)計(jì)算，并獲得了出其意料的加速效果18。然而由于顯卡編程難度大，早期GPU計(jì)算并未得到廣泛關(guān)注。隨著NVIDIA公司推出CUDA編程語(yǔ)言，大大降低了GPU編程的難度，同時(shí)還推出了具有更強(qiáng)浮點(diǎn)運(yùn)算能力的芯片，使得GPU在科學(xué)計(jì)算領(lǐng)域有了更大的應(yīng)用空間。圖1GPU與CPU的浮點(diǎn)運(yùn)算速度計(jì)算流體力學(xué)（CFD）屬于計(jì)算密集型領(lǐng)域，Navier-Stokes方程的求解需要進(jìn)行大量的浮點(diǎn)運(yùn)算。同時(shí)由于流動(dòng)結(jié)構(gòu)的復(fù)雜性，要捕捉足夠細(xì)致的流場(chǎng)結(jié)構(gòu)需要足夠密集的網(wǎng)格。具有實(shí)際應(yīng)用價(jià)值的計(jì)算甚至需要108數(shù)量級(jí)的網(wǎng)格，其計(jì)算量是驚人的，普通的計(jì)算機(jī)難以勝任，必須進(jìn)行并行計(jì)算。考慮到GPU具有極強(qiáng)的并行計(jì)算能力，很多學(xué)者嘗試著將GPU引進(jìn)到計(jì)算流體力學(xué)中。Cohen3等實(shí)現(xiàn)了雙精度的CFD計(jì)算。Tolke20實(shí)現(xiàn)了格子波爾茲曼方法在GPU上的加速。AliKhajeh-Saeed15使用GPU進(jìn)行了各項(xiàng)同性湍流直接數(shù)值模擬。針對(duì)存在激波間斷的可壓縮流動(dòng)，Antoniou1在GPU上實(shí)現(xiàn)了WENO格式的加速，獲得了較大的加速比例。VahidEsfahanian9實(shí)現(xiàn)了一維和二維WENO格式的GPU程序，其一維的加速效果達(dá)到驚人的1000倍。為實(shí)現(xiàn)大規(guī)模的直接數(shù)值模擬，F(xiàn)rancescoSalvadore19用多塊GPU進(jìn)行了混合層流動(dòng)的直接數(shù)值模擬。Karantasis14則比較了不同形式并行對(duì)GPU計(jì)算效率的影響。對(duì)于更復(fù)雜的結(jié)構(gòu)，Elsen7首次在在GPU上實(shí)現(xiàn)了有限體積方法對(duì)復(fù)雜外形超音速流動(dòng)的計(jì)算。Fu10在GPU上實(shí)現(xiàn)了多種格式的有限體積方法計(jì)算，并進(jìn)行了大量驗(yàn)算。就數(shù)值方法而言，以上提到的文獻(xiàn)中多采用顯示迎風(fēng)格式進(jìn)行空間離散。相比于迎風(fēng)格式，緊致格式具有更高的譜分辨率和低耗散的特點(diǎn)，因此更適合于湍流計(jì)算。然而，使用緊致格式時(shí)，需要求解代數(shù)方程組。在GPU上實(shí)現(xiàn)代數(shù)方程組的求解并不簡(jiǎn)單。BulentTutkun21在GPU上實(shí)現(xiàn)了緊致格式并將其用于流動(dòng)計(jì)算中，其計(jì)算時(shí)間的40%花費(fèi)在代數(shù)方程求解上。為提高緊致格式的計(jì)算效率，VahidEsfahanian8使用了一種改進(jìn)的三對(duì)角方程求解器，并獲得了一定的加速效果，然而其方法受到網(wǎng)格點(diǎn)數(shù)的限制。從數(shù)學(xué)角度來(lái)說(shuō)，三對(duì)角方程并不難以求解，經(jīng)典的Tomas方法是現(xiàn)行最有效的求解方法。然而由于GPU的高度并行結(jié)構(gòu)及有限的內(nèi)存帶寬等，人們希望并行的求解方程組從而進(jìn)一步提高計(jì)算效率。從這一點(diǎn)出發(fā)，CyclicReduction被提出來(lái)，由該算法并行的特點(diǎn)，被認(rèn)為非常適合使用在GPU上。Zhang23采用全局共享內(nèi)存的模擬在GPU上實(shí)現(xiàn)了并行CR方法。隨后一種混合PCR-Tomas方法被提出來(lái)并被證明性能優(yōu)于PCR方法。Davidson6提出了一種基于寄存器而非共享內(nèi)存的CR方法，其性能也有一定的提升。Kim等16提出了一種新的PCR-Tomas方法，其方法在方程數(shù)較大時(shí)，完全采用Tomas方法求解。然而，從Kim等的結(jié)果可以看出來(lái)，在方程組數(shù)量較大，方程中未知數(shù)較多時(shí)，PCR-Tomas方法的計(jì)算效率并不如單純的Tomas方法?？紤]到計(jì)算流體力學(xué)中，計(jì)算網(wǎng)格通常較大，因此，直接采用Tomas方法而非CR方法更為有效。本文在GPU上實(shí)現(xiàn)了二維可壓縮Navier-Stokes方程的數(shù)值計(jì)算。將Steger-Warming分裂、五階WENO格式13和加權(quán)緊致格式（CRWENO）格式19以及三階TVD龍格-庫(kù)塔格式移植到GPU上。緊致格式所產(chǎn)生的三對(duì)角方程組采用直接Tomas方法求解。第二節(jié)中簡(jiǎn)要介紹了GPU的構(gòu)架及其對(duì)計(jì)算效率的影響；第三節(jié)中描述了如何在GPU上實(shí)現(xiàn)WENO格式及CRWENO；第四節(jié)中展示了采用GPU進(jìn)行的幾個(gè)數(shù)值算例并與采用相同方法的CPU串行/并行程序進(jìn)行了比較。最后一節(jié)對(duì)全文進(jìn)行了總結(jié)并對(duì)未來(lái)的研究方向進(jìn)行了簡(jiǎn)述。此外，本文中所實(shí)現(xiàn)的GPU程序全部采用雙精度數(shù)據(jù)。1控制方程本文中控制方程為無(wú)量綱化的Navier-Stokes方程：其中：，為求解控制方程，對(duì)流項(xiàng)采用Steger-Warming方法進(jìn)行通量分裂，而后采用5階WENO/CRWENO格式離散，粘性項(xiàng)采用四階中心差分格式。時(shí)間項(xiàng)推進(jìn)采用3階TVD龍格庫(kù)塔方法。2數(shù)值方法對(duì)于對(duì)流問(wèn)題，數(shù)值格式的迎風(fēng)性既能保證解的物理性又能保證數(shù)值計(jì)算的穩(wěn)定性，本文的計(jì)算中采用Seger-Warming格式對(duì)對(duì)流項(xiàng)進(jìn)行分裂。對(duì)于通量分裂后的對(duì)流項(xiàng)，為求其一階倒數(shù)，本文采用Jiang等12提出的WENO-JS格式。以正通量為例，考慮網(wǎng)格單元交界面處的流動(dòng)通量：其中，與三階迎風(fēng)格式比較可以發(fā)現(xiàn)，WENO格式可以看做9次三階格式計(jì)算，因此其計(jì)算量至少是三階格式的九倍。如果考慮二維N-S方程的正負(fù)通量以及兩個(gè)通量矢量和四個(gè)方程，其增加的計(jì)算量是可觀的。相比于WENO格式，CRWENO格式具有高的普分辨率和低耗散的特點(diǎn)。其流動(dòng)通量的線(xiàn)性表達(dá)式為：考慮到五階線(xiàn)性格式可以表達(dá)為一下三個(gè)三階格式的加權(quán)疊加：其系數(shù)分別為。將這些加權(quán)系數(shù)替換成WENO格式中的加權(quán)函數(shù)便可獲得激波捕捉的能力。為保證數(shù)值解的無(wú)振蕩性，時(shí)間推進(jìn)采用三階TVD龍格庫(kù)塔格式：3GPU結(jié)構(gòu)相比于CPU，GPU有著完全不同的硬件、執(zhí)行和內(nèi)存構(gòu)架，為了設(shè)計(jì)出高效率的GPU程序，必須對(duì)此有比較深入的理解。以下簡(jiǎn)單介紹一下NVIDIAGPU的一些基本概念和結(jié)構(gòu)。在硬件層面，每個(gè)GPU由若干個(gè)SM（多線(xiàn)程處理器）組成，這些SM相當(dāng)于CPU中核（Core）的概念。每個(gè)SM還擁有若干個(gè)SP（流處理器或CUDA核），這些處理器是GPU運(yùn)算的基本單元。以本文中所使用的NVIDIAGTX750Ti為例，總共有五個(gè)多線(xiàn)程處理器，每個(gè)處理器包含了128個(gè)CUDA核。圖2CPU的結(jié)構(gòu)與GPU的結(jié)構(gòu)從運(yùn)行角度來(lái)說(shuō)，GPU運(yùn)算的最基本結(jié)構(gòu)成為線(xiàn)程（Thread），一個(gè)或多個(gè)線(xiàn)程組成一個(gè)線(xiàn)程塊（ThreadBlock），GPU中的線(xiàn)程以塊為單元，被分配到各個(gè)處理器上。在每個(gè)線(xiàn)程塊中，所有的線(xiàn)程同時(shí)執(zhí)行相同的運(yùn)算，而不同的線(xiàn)程塊雖然也是執(zhí)行同樣的運(yùn)算但是實(shí)際上并不是同時(shí)執(zhí)行的。在分配Block的過(guò)程中，GPU并不按照邏輯順序執(zhí)行分配，而是將隨機(jī)的Block分配到各個(gè)處理器上。也就是說(shuō)GPU實(shí)際的并行是在Block內(nèi)部，認(rèn)識(shí)到這一點(diǎn)才能設(shè)計(jì)正確的程序。對(duì)于NVIDIAGTX750Ti來(lái)說(shuō)，每個(gè)線(xiàn)程塊最多包含2048個(gè)線(xiàn)程，這也是每個(gè)SM能同時(shí)處理的最大線(xiàn)程數(shù)。在執(zhí)行過(guò)程中，每個(gè)Block中的線(xiàn)程又以32個(gè)為一組分成不同的warp。在處理器中，這些warp被同時(shí)執(zhí)行，與線(xiàn)程塊中線(xiàn)程的同時(shí)執(zhí)行不同的是，線(xiàn)程塊中的線(xiàn)程同時(shí)開(kāi)始但可能并不是同時(shí)完成某項(xiàng)任務(wù)，而warp中的所有線(xiàn)程同時(shí)開(kāi)始并同時(shí)結(jié)束。如果線(xiàn)程執(zhí)行過(guò)程中存在邏輯分支，如if-else結(jié)構(gòu)。那么執(zhí)行if中運(yùn)算的線(xiàn)程將繼續(xù)執(zhí)行而不滿(mǎn)足條件的進(jìn)程將被擱置直到前面這些線(xiàn)程完成運(yùn)算。隨后等待的這些線(xiàn)程執(zhí)行else結(jié)構(gòu)中的運(yùn)算而其他的線(xiàn)程繼續(xù)等待。從以上可以看出，由于存在等待機(jī)制，GPU不擅長(zhǎng)處理存在大量邏輯分支的運(yùn)算。之所以這樣設(shè)計(jì)是出于節(jié)能的考慮，具體原理可參考NVIDIA的相關(guān)說(shuō)明書(shū)18。圖3GPU的線(xiàn)程結(jié)構(gòu)GPU的內(nèi)存有三個(gè)層次，包括全局內(nèi)存（Globalmemory）、共享內(nèi)存(Sharedmemory)和寄存器(Register)。全局內(nèi)存是整個(gè)GPU共享的內(nèi)存空間，既我們平時(shí)稱(chēng)的顯存。相比于其他兩種內(nèi)存其特點(diǎn)是存儲(chǔ)量大但是讀取速度慢。由于全局內(nèi)存讀取速度較慢，為了提高效率，GPU中有一種獨(dú)特的內(nèi)存讀取方式成為一致讀取，既當(dāng)一個(gè)Warp中的所有線(xiàn)程讀取內(nèi)存中連續(xù)的數(shù)據(jù)時(shí)，GPU將通過(guò)某種機(jī)制快速的將這些數(shù)據(jù)發(fā)送到各個(gè)線(xiàn)程。考慮到GPU內(nèi)存讀取的這一特點(diǎn)，在設(shè)計(jì)數(shù)據(jù)結(jié)構(gòu)時(shí)應(yīng)該到使用數(shù)組結(jié)構(gòu)（structureofarrays）而不是通常所使用的結(jié)構(gòu)數(shù)組(arrayofstructures)。共享內(nèi)存屬于片上內(nèi)存,是每個(gè)SM所獨(dú)有的，每個(gè)SM擁有48k的共享內(nèi)存。之所以稱(chēng)為共享內(nèi)存是因?yàn)閷?duì)于同一個(gè)Block中的線(xiàn)程，這個(gè)存儲(chǔ)空間是被共享的。一個(gè)Block中的線(xiàn)程不能訪(fǎng)問(wèn)另一個(gè)Block中的共享內(nèi)存。由于這一個(gè)特點(diǎn)，如果程序中大量使用到了共享內(nèi)存，同時(shí)還要權(quán)衡一個(gè)線(xiàn)程塊的大小和每個(gè)SM上最多分配的線(xiàn)程數(shù)量。如果一個(gè)Block使用了SM中的全部共享內(nèi)存，則即使這個(gè)block中的線(xiàn)程數(shù)量小于上線(xiàn)，其它block也不能被執(zhí)行直到該Block執(zhí)行完畢。寄存器是在運(yùn)算過(guò)程中用于存儲(chǔ)變量的結(jié)構(gòu)，是幾種內(nèi)存中讀取速度最快的。當(dāng)然每個(gè)SM上的寄存器數(shù)量也是有限的，NVIDIAGTX750Ti的每個(gè)SM總共有65536個(gè)寄存器可以使用。雖然整體的數(shù)量較大，但是考慮到每個(gè)SM最多處理2048個(gè)線(xiàn)程，在此情況下分配到每個(gè)線(xiàn)程上的寄存器數(shù)量?jī)H為32個(gè)?，F(xiàn)有的GPU構(gòu)架中，一個(gè)寄存器的大小為32-bit，也就是一個(gè)單精度浮點(diǎn)數(shù)據(jù)的大小，如果使用雙精度數(shù)據(jù)，所使用的寄存器將翻倍。寄存器的分配不由程序員控制，而是由編譯器分配和優(yōu)化。因此，要限制寄存器的使用量需要一定的經(jīng)驗(yàn)。圖4GPU的內(nèi)存結(jié)構(gòu)考慮到寄存器的讀取速度要快于共享內(nèi)存和全局內(nèi)存，因此除了需要反復(fù)更新數(shù)據(jù)的地方用到共享內(nèi)存外，所有數(shù)據(jù)從全局內(nèi)存中獲取后將一致被存儲(chǔ)在寄存器中。由于使用到了大量的寄存器，而GPU每一個(gè)處理器的寄存器數(shù)量是有限的，因此會(huì)降低SM的占有率。然而通過(guò)比較使用不同方式的數(shù)據(jù)存儲(chǔ)方式發(fā)現(xiàn)，即使在使用共享內(nèi)存達(dá)到100%占用率時(shí)，其速度仍然要慢于大量使用寄存器時(shí)的速度。另一個(gè)角度來(lái)說(shuō)，由于雙精度數(shù)據(jù)需要使用到兩個(gè)32-bit寄存器，在進(jìn)行大量浮點(diǎn)運(yùn)算時(shí)，難以將寄存器數(shù)量降低到實(shí)現(xiàn)100%占用率的水平。因此與其花費(fèi)大量時(shí)間優(yōu)化使用共享內(nèi)存的程序，不如直接使用寄存器。4程序構(gòu)架基于上一節(jié)所介紹的GPU結(jié)構(gòu)，作者設(shè)計(jì)了具有較高效率的CUDA程序。本節(jié)將介紹本文中所使用程序的結(jié)構(gòu)。由于存在一致讀取的機(jī)制，本文程序中將不同網(wǎng)格點(diǎn)上同一變量的數(shù)據(jù)存儲(chǔ)在一個(gè)數(shù)組中，并按照網(wǎng)格的下標(biāo)順序排列。程序整體上分為通量分裂、通量重構(gòu)和時(shí)間推進(jìn)三個(gè)模塊，每一個(gè)模塊對(duì)應(yīng)一組kernel函數(shù)。對(duì)于通量分裂模塊，由于大量重復(fù)的使用一些變量，例如密度、聲速等，為避免重復(fù)讀取內(nèi)存到來(lái)的延遲，程序中只使用一個(gè)kernel函數(shù)進(jìn)行全部通量的分裂。在重構(gòu)部分，由于正負(fù)通量存儲(chǔ)在不同的數(shù)組中，并且重構(gòu)所使用到的數(shù)據(jù)范圍不同——對(duì)于正通量WENO格式使用i-2到i+2五個(gè)點(diǎn)，而負(fù)通量使用i-1到i+3個(gè)點(diǎn)，每個(gè)方向（x和y方向）分辨有兩個(gè)kernel。為了充分利用GPU，在kernel中，每個(gè)線(xiàn)程需要處理流動(dòng)通量矢量中的四個(gè)分量，這相當(dāng)于每個(gè)線(xiàn)程進(jìn)行四次重構(gòu)。對(duì)于CRWENO格式，為節(jié)約顯存，計(jì)算所得到的代數(shù)方程系數(shù)被存儲(chǔ)在之后才會(huì)被其它模塊使用的內(nèi)存中去（如粘性通量數(shù)組）。三對(duì)角方程的Tomas求解分為兩個(gè)步驟，向前消元和回帶過(guò)程，其中消元過(guò)程為：回帶過(guò)程為：可以看到，在整個(gè)計(jì)算過(guò)程中需要反復(fù)的讀取和存放數(shù)據(jù)。在計(jì)算中，每個(gè)方程組被映射給一個(gè)線(xiàn)程。由于擁有大量的方程組，為了減小數(shù)據(jù)存取帶來(lái)的延遲，需要盡量使用一致讀取的方式來(lái)獲得數(shù)據(jù)。然而，對(duì)于x方向來(lái)說(shuō)，方程組的各個(gè)系數(shù)是存儲(chǔ)在連續(xù)的內(nèi)存地址中的（本文中所有數(shù)據(jù)按照x變化方向存儲(chǔ)，既所有y相同的數(shù)據(jù)存儲(chǔ)在連續(xù)分部的內(nèi)存中），如圖a所示。為了實(shí)現(xiàn)一致內(nèi)存訪(fǎng)問(wèn)，在存儲(chǔ)x方向方程組各個(gè)系數(shù)時(shí)，將系數(shù)的下標(biāo)（i，j）對(duì)換，這樣可以得到如圖b所示的存儲(chǔ)結(jié)構(gòu)。在回帶過(guò)程中由于讀取的數(shù)據(jù)量小，可以將計(jì)算所得的數(shù)據(jù)按照原來(lái)的方式存儲(chǔ)到流通量的數(shù)組中去。（a）非一致內(nèi)存訪(fǎng)問(wèn)（b）一致內(nèi)存訪(fǎng)問(wèn)圖5Tomas方法求解示意圖相比于文獻(xiàn)中的各種CR方法，本文直接采用Tomas方法具有簡(jiǎn)單有效、可擴(kuò)展性強(qiáng)不受共享內(nèi)存大小限制的特點(diǎn)。對(duì)于大網(wǎng)格的計(jì)算，直接Tomas方法具有更高的計(jì)算效率。時(shí)間推進(jìn)由于沒(méi)有復(fù)雜的算法，因此其程序的結(jié)構(gòu)與CPU版本基本一致。4數(shù)值結(jié)果與比較4.1歐拉方程為驗(yàn)證程序的正確性，本文首先測(cè)試了程序?qū)τ诠饣鲌?chǎng)的計(jì)算。所用測(cè)試初始條件為：既渦對(duì)流問(wèn)題，其中來(lái)流速度為u=0.5，v=0,b=0.5。在本問(wèn)題中一個(gè)漩渦周期性的通過(guò)計(jì)算區(qū)域。在時(shí)間t=20時(shí)考察流場(chǎng)的密度分布。由于流動(dòng)是無(wú)粘的，因此理論上渦及其周?chē)鲌?chǎng)結(jié)構(gòu)不應(yīng)發(fā)生變化，但由于格式存在耗散性，渦的強(qiáng)度會(huì)衰減。如圖所示，在一個(gè)周期后渦的結(jié)構(gòu)和密度分部保持良好，可以認(rèn)為程序是正確的。圖6t=20時(shí)的密度分部（CRWENO）在驗(yàn)證了程序正確的基礎(chǔ)上，比較了采用相同算法的GPU程序和CPU程序。測(cè)試中所使用的是IntelXeonE5-26092.4G處理器，并行程序采用openmp方式。圖7所示為不同版本程序在不同網(wǎng)格點(diǎn)數(shù)下時(shí)間推進(jìn)一次所需的計(jì)算時(shí)間?？梢钥闯?，單核CPU程序計(jì)算時(shí)間增加呈線(xiàn)性增長(zhǎng)，而多核時(shí)由于在小網(wǎng)格時(shí)未充分利用CPU的大cache的特點(diǎn)計(jì)算時(shí)間增長(zhǎng)不太線(xiàn)性，當(dāng)網(wǎng)格點(diǎn)大于256時(shí)，表現(xiàn)為線(xiàn)性的時(shí)間增長(zhǎng)。對(duì)于GPU程序，在小網(wǎng)格數(shù)量時(shí)，GPU性能不能被充分釋放，隨著網(wǎng)格數(shù)的增加，計(jì)算時(shí)間逐漸趨于線(xiàn)性，表明GPU被充分利用。可以看出，WENO算法和CRWENO算法隨消耗的時(shí)間隨著網(wǎng)格數(shù)量的增加逐漸接近。在N=1024時(shí)，CRWENO算法所需要的時(shí)間大概是WENO算法的1.7倍。圖7每個(gè)時(shí)間步長(zhǎng)的計(jì)算時(shí)間圖8比較了在不同網(wǎng)格點(diǎn)下GPU相對(duì)于CPU的加速效果。從圖中可以看出，相比于串行的CPU程序，對(duì)于WENO格式，本文中的GPU程序達(dá)到了23倍的加速效果，而相比于8核CPU并行程序，達(dá)到了4倍的加速效果。對(duì)于CRWENO格式，比較了在單核下的加速比，可以看到加速比對(duì)于網(wǎng)格數(shù)量呈現(xiàn)出較強(qiáng)的線(xiàn)性標(biāo)度，當(dāng)網(wǎng)格點(diǎn)數(shù)足夠大的時(shí)候，GPU的WENO算法與CRWENO算法加速效果接近，大約20倍。圖8不同CPU核數(shù)和網(wǎng)格點(diǎn)下的加速比(WENO)對(duì)于流場(chǎng)中存在激波的情況，本文計(jì)算了雙馬赫反射問(wèn)題。由于在上卷區(qū)域和接觸線(xiàn)上存在小尺度結(jié)構(gòu)，計(jì)算中采用較大的網(wǎng)格數(shù)量（NxxNy=2048*512）。圖給出了t=0.2時(shí)上卷區(qū)域的密度分布。可以看出，CRWENO格式給出了分辨率更高的結(jié)構(gòu)。在接觸線(xiàn)上，CRWENO格式的結(jié)果有更多的小尺度結(jié)構(gòu)，說(shuō)明其耗散小。a．WENOb.Crweno圖9雙馬赫反射整體密度分部WENOCRWENO圖10雙馬赫反射上卷區(qū)域密度分部4.2Navier-Stokes方程當(dāng)存在粘性時(shí)，流動(dòng)容易失穩(wěn)成為湍流。二維湍流與三維湍流在物理上有著很大的區(qū)別，但他們都存在小尺度結(jié)構(gòu)，因此網(wǎng)格尺度必須足夠小才能準(zhǔn)確的將這些結(jié)構(gòu)捕捉。本文測(cè)試了激波剪切層干擾問(wèn)題。一個(gè)馬赫數(shù)為0.6，雷諾數(shù)為500的剪切層與一道斜激波發(fā)生相互作用，剪切穿過(guò)斜激波及其在壁面上的反射波。由于初始條件中包含了擾動(dòng)，剪切層在與激波發(fā)生作用前已失穩(wěn)產(chǎn)生渦，這些渦隨后與激波發(fā)生碰撞。WENOCRWENO圖11激波剪切層干擾問(wèn)題密度云圖計(jì)算中網(wǎng)格點(diǎn)數(shù)為NxxNy=1024*256。圖11給出了t=120時(shí)的流暢結(jié)構(gòu)，從圖中可以看到精細(xì)的渦結(jié)構(gòu)和明顯的激波間斷現(xiàn)象。相比于CPU程序，GPU計(jì)算此算例總用時(shí)約為1個(gè)小時(shí)，如果采用單核CPU計(jì)算將需要20個(gè)小時(shí)左右，而如果采用8核并行則仍然需要4小時(shí)左右的時(shí)間。由此可以看出采用GPU進(jìn)行流動(dòng)計(jì)算能夠大大算短計(jì)算周期，甚至使很多之前由于計(jì)算時(shí)間過(guò)長(zhǎng)而無(wú)法實(shí)現(xiàn)的算例變的可行。5結(jié)論本文在單塊GPU上實(shí)現(xiàn)了可壓縮流動(dòng)計(jì)算，獲得了可觀的加速效果。本文中所使用到的GPU價(jià)格低廉，4塊NVIDIAGTX750TI的價(jià)格才與8核CPU相當(dāng)，而其計(jì)算效率是后者的20倍。現(xiàn)有的技術(shù)已經(jīng)能實(shí)現(xiàn)多塊GPU的并行計(jì)算，其加速效果更為明顯。從本文的比較結(jié)果可以看出，在未來(lái)CFD計(jì)算中GPU有著很大的潛力和市場(chǎng)，是流體計(jì)算的一種強(qiáng)有力的工具。參考文獻(xiàn)Antoniou,A.S.,Karantasis,K.I.,Polychronopoulos,E.D.,&Ekaterinaris,J.A.(2010).AccelerationofaFinite-DifferenceWENOSchemeforLarge-ScaleSimulationsonMany-CoreArchitectures.Chang,L.-W.,Stratton,J.A.,Kim,H.-S.,&Hwu,W.W.(2012).Ascalable,numericallystable,high-performancetridiagonalsolverusingGPUs.Cohen,J.,&Molemake,J.(2009).AFastDoublePrecisionCFDCodeUsingCUDA.Paperpresentedatthe21stInternationalConferenceonParallelComputationalFluidDynamics(ParCFD2009).Corrigan,A.,Camelli,F.,L?hner,R.,&Mut,F.(2012).Semi-automaticportingofalarge-scaleFortranCFDcodetoGPUs.InternationalJournalforNumericalMethodsinFluids,69(2),314-331.doi:10.1002/fld.2560Corrigan,A.,&L?hner,R.(2012).Semi-automaticportingofalarge-scaleCFDcodetomulti-graphicsprocessingunitclusters.InternationalJournalforNumericalMethodsinFluids,69(11),1786-1796.doi:10.1002/fld.2664Davidson,A.,&Owens,J.D.(2011).Registerpackingforcyclicreduction:acasestudy.Elsen,E.,LeGresley,P.,&Darve,E.(2008).LargecalculationoftheflowoverahypersonicvehicleusingaGPU.JournalofComputationalPhysics,227(24),10148-10161.doi:10.1016/j.jcp.2008.08.023Esfahanian,V.,Baghapour,B.,Torabzadeh,M.,&Chizari,H.(2014).AnefficientGPUimplementationofcyclicreductionsolverforhigh-ordercompressibleviscousflowsimulations.Computers&Fluids,92,160-171.doi:10.1016/pfluid.2013.12.011Esfahanian,V.,Darian,H.M.,&Gohari,S.M.I.(2013).AssessmentofWENOschemesfornumericalsimulationofsomehyperbolicequationsusingGPU.Computers&Fluids,80,260-268.doi:10.1016/pfluid.2012.02.031Fu,L.,Gao,Z.,Xu,K.,&Xu,F.(2014).Amulti-blockviscousflowsolverbasedonGPUparallelmethodology.Computers&Fluids,95,19-39.doi:10.1016/pfluid.2014.02.005Ghosh,D.,&Baeder,J.D.(2012).CompactReconstructionSchemeswithWeightedENOLimitingforHyperbolicConservationLaws.SIAMJournalonScientificComputing,34(3),A1678-A1706.doi:10.1137/110857659ImanGohari,S.M.,Esfahanian,V.,&Moqtaderi,H.(2013).CoalescedcomputationsoftheincompressibleNavier–Stokesequationsoveranairfoilusinggraphicsprocessingunits.Computers&Fluids,80,102-115.doi:10.1016/pfluid.2012.04.022Jiang,G.-S.,&Shu,C.-W.(1996).EfficientImplementationofWeightedENOSchemes.JournalofComputationalPhysics,126(1),202-228.doi:/10.1006/jcph.1996.0130Karantasis,K.I.,Polychronopoulos,E.D.,&Ekaterinaris,J.A.(2014).HighorderaccuratesimulationofcompressibleflowsonGPUclustersoverSoftwareDistributedSharedMemory.Computers&Fluids,93,18-29.doi:10.1016/pfluid.2014.01.005Khajeh-Saeed,A.,&BlairPerot,J.(2013).DirectnumericalsimulationofturbulenceusingGPUacceleratedsupercomputers.JournalofComputationalPhysics,235,241-257.doi:10.1016/j.jcp.2012.10.050Kim,H.-S.,Wu,S.,Chang,L.-w.,&Hwu,W.-m.W.(2011).AScalableTridiagonalSolverforGPUs.444-453.doi:10.1109/icpp.2011.41Kuo,F.-A.,Smith,M.R.,Hsieh,C.-W.,Chou,C.-Y.,&Wu,J.-S.(2011).GPUaccelerationforgeneralconservationequationsanditsapplicationtoseveralengineeringproblems.Computers&Fluids,45(1),147-154.doi:10.1016/pfluid.2010.10.007Owens,J.D.,Luebke,D.,Govindaraju,N.,Harris,M.,Krüger,J.,Lefohn,A.E.,&Purcell,T.J.(2007).ASurveyofGeneral-PurposeComputationonGraphicsHardware.ComputerGraphicsForum(1).Salvadore,F.,Bernardini,M.,&Botti,M.(2013).GPUacceleratedflowsolverfordirectnumericalsimulationofturbulentflows.JournalofComputationalPhysics,235,129-142.doi:10.1016/j.jcp.2012.10.012T?lke,J.,&Krafczyk,M.(2008).TeraFLOPcomputingonadesktopPCwithGPUsfor3DCFD.InternationalJournalofComputationalFluidDynamics,22(7),443-456.doi:10.1080/10618560802238275Tutkun,B.,&Edis,F.O.(2012).AGPUapplicationforhigh-ordercompactfinitedifferencescheme.Computers&Fluids,55,29-35.doi:10.1016/pfluid.2011.10.016Waltz,J.(2013).Performanceofathree-dimensionalunstructuredmeshcompressibleflowsolveronNVIDIAFermi-classgraphicsprocessingunithardware.InternationalJournalforNumericalMethodsinFluids,72(2),