陳 龍,徐添豪,田書玲

(南京航空航天大學(xué) 航空宇航學(xué)院,南京 210016)

近年來(lái),圖形處理器(GPU)在科學(xué)計(jì)算領(lǐng)域中發(fā)揮了重要的作用.由于其出色的浮點(diǎn)運(yùn)算性能以及具有優(yōu)勢(shì)的顯存帶寬,GPU往往能給科學(xué)計(jì)算帶來(lái)可觀的加速.計(jì)算流體力學(xué)(CFD)領(lǐng)域的GPU加速已得到廣泛的研究,Hagen等人[1,2]先后在GPU上實(shí)現(xiàn)了基于歐拉方程以及基于NS方程的求解器,得到了良好的加速效果.隨著GPU硬件自身的發(fā)展,更復(fù)雜的算法得到了實(shí)現(xiàn),基于RANS方程、LES方法以及DNS方法的國(guó)內(nèi)外先進(jìn)研究成果不斷被發(fā)表[3–7].為了進(jìn)一步地提高運(yùn)算能力,基于多GPU的實(shí)現(xiàn)也得到了越來(lái)越多的關(guān)注[8].劉楓等人[9]基于MPI+CUDA的異構(gòu)并行可壓縮流求解器,用于高超聲速流動(dòng)的數(shù)值模擬計(jì)算效率較CPU同構(gòu)并行計(jì)算提高10倍以上.采用GPU加速的CFD算法在航空航天工程中具有廣闊應(yīng)用前景.

CFD算法GPU加速目前的研究前沿主要集中在結(jié)構(gòu)網(wǎng)格和顯式時(shí)間推進(jìn).非結(jié)構(gòu)網(wǎng)格相比于結(jié)構(gòu)網(wǎng)格更適用于復(fù)雜的幾何外形,但由于其存儲(chǔ)結(jié)構(gòu)的特點(diǎn)難其內(nèi)存存取的連續(xù)性較差,在GPU架構(gòu)上難以獲得很高的效率.同時(shí),以LU-SGS算法為代表的隱式方法具有較顯式方法高數(shù)倍的計(jì)算效率,但其每次迭代過程中多次需要相鄰單元數(shù)據(jù),需要大量的內(nèi)存訪問,在GPU架構(gòu)上需要較好的數(shù)據(jù)局部性才能得到較高的效率,因此通常非結(jié)構(gòu)網(wǎng)格的隱式算法在GPU上的加速效果不佳.

針對(duì)這一問題,本文發(fā)展了LU-SGS算法在非結(jié)構(gòu)網(wǎng)格上的高效率的GPU實(shí)現(xiàn)方法.首先,介紹了基本控制方程和LU-SGS算法.然后,給出了算法的GPU實(shí)現(xiàn)方法.并采用RCM以及基于Metis網(wǎng)格分區(qū)的網(wǎng)格重排序方法優(yōu)化算法的加速效果,詳細(xì)討論各網(wǎng)格排序算法在GPU實(shí)現(xiàn)中的適用性.最后,通過算例驗(yàn)證了算法的正確性以及實(shí)現(xiàn)方案的高效性.

1 數(shù)值方法

慣性系下,忽略源項(xiàng)的RANS方程形式如下:

其中W,F,G分別代表守恒變量、對(duì)流通量及粘性通量.本文采用非結(jié)構(gòu)網(wǎng)格下格點(diǎn)格式離散有限體積控制體,使用Roe格式離散對(duì)流通量.通過線性重構(gòu)實(shí)現(xiàn)空間二階精度.使用收斂性良好的Venkatakrishnan限制器.采用隱式LU-SGS方法進(jìn)行時(shí)間推進(jìn),其形式如下[10]

其中,U,L和D分別為嚴(yán)格上三角矩陣,嚴(yán)格下三角矩陣和對(duì)角矩陣.

采用向前掃掠和向后掃掠求解該方程,過程可以寫成:

向前掃掠:

向后掃掠:

該方法最早由Jameson和Yoon[11]應(yīng)用于結(jié)構(gòu)網(wǎng)格,在掃掠時(shí)確定超平面,如圖1.同一超平面內(nèi)的控制體可以同步更新而不存在相互的依賴,這也為計(jì)算的數(shù)據(jù)并行提供了理論基礎(chǔ).后Nakahashi等人[12]將其應(yīng)用到非結(jié)構(gòu)網(wǎng)格上,掃掠時(shí)也需要確定類似形式的超平面,如圖2.

圖1 結(jié)構(gòu)網(wǎng)格上的超平面形式

圖2 非結(jié)構(gòu)網(wǎng)格上的超平面形式

湍流模型選用Spalart-Allmaras一方程模型[13].該模型廣泛應(yīng)用于飛機(jī)機(jī)翼、翼身組合體等繞流的數(shù)值模擬中.

2 算法 GPU 并行加速

2.1 GPU并行方法

單個(gè)GPU,從硬件上來(lái)說(shuō),由多個(gè)流多處理器SM組成.每個(gè)SM可以同時(shí)處理上千個(gè)線程.在英偉達(dá)公司提出的CUDA框架下,線程以32為一束進(jìn)行調(diào)度,即同一線程束(Warp)內(nèi)的線程在同一時(shí)間執(zhí)行相同的指令,這也對(duì)應(yīng)著并行的最細(xì)粒度.而若干個(gè)Warp又組成線程塊(Block)被分發(fā)到不同SM上執(zhí)行.每個(gè)計(jì)算任務(wù)會(huì)對(duì)應(yīng)一個(gè)線程網(wǎng)格(Grid)進(jìn)行線程組織,而一個(gè)Grid是由多個(gè)Block組成的.從Grid到Block到Thread都可以在編程時(shí)進(jìn)行一維、二維乃至三維的索引,如圖3所示.

運(yùn)行于GPU上的并行程序是基于數(shù)據(jù)并行的,每個(gè)線程執(zhí)行相同的任務(wù),程序員也無(wú)需過多地關(guān)心線程的調(diào)度方式,它由 GPU 自行處理.因此,在 CFD 領(lǐng)域內(nèi),涉及大量數(shù)據(jù)并行的CPU程序很容易移植到GPU上.在傳統(tǒng)的CPU求解器中,非結(jié)構(gòu)網(wǎng)格結(jié)合格點(diǎn)格式會(huì)涉及兩種形式的循環(huán),即基于點(diǎn)的和基于邊的.對(duì)于這樣兩種循環(huán),我們?cè)贕PU上最為直接地將單個(gè)點(diǎn)或者單條邊映射到單個(gè)線程上,而每個(gè)線程則執(zhí)行原來(lái)循環(huán)體內(nèi)的指令.

圖3 線程、線程塊、線程塊網(wǎng)絡(luò)的組織形式

可以寫成對(duì)應(yīng)的由CPU調(diào)用、GPU執(zhí)行的函數(shù)(Kernel函數(shù))的形式

代碼區(qū)域則幾乎相同,這為快速的代碼移植提供了相當(dāng)?shù)谋憬菪?可見GPU的細(xì)粒度并行方式通常與CPU上的多線程、多進(jìn)程實(shí)現(xiàn)的任務(wù)并行不同(將單個(gè)任務(wù)拆分成多個(gè)子任務(wù)).當(dāng)然在單個(gè)或者多個(gè)GPU上進(jìn)行類似于CPU的任務(wù)并行也是實(shí)現(xiàn)大規(guī)模GPU運(yùn)算的必要手段.

2.2 算法的GPU實(shí)現(xiàn)

涉及算法在GPU上的流程圖如圖4.

圖4 CPU 及 GPU 算法流程圖

GPU流程與CPU上的實(shí)現(xiàn)幾乎完全相同.在GPU上進(jìn)行計(jì)算的前后,需要將數(shù)據(jù)拷貝到GPU上或者從GPU拷貝回CPU,但是所有的求解過程是在GPU上進(jìn)行的,如圖中紅色虛線框.可以認(rèn)為GPU上實(shí)現(xiàn)的每一個(gè)步驟對(duì)應(yīng)一個(gè)Kernel函數(shù).不同的是,LU-SGS向前、向后掃掠步驟中涉及了數(shù)百甚至上千個(gè)Kernel函數(shù)(與超平面?zhèn)€數(shù)有關(guān)).LU-SGS算法本身就適合于進(jìn)行超平面內(nèi)的數(shù)據(jù)并行.而進(jìn)行掃掠的時(shí)候,下一個(gè)超平面內(nèi)的變量依賴于更新過后的當(dāng)前面內(nèi)的變量.因此,可以采用每一層對(duì)應(yīng)一個(gè)Kernel函數(shù),而這些Kernel函數(shù)串行執(zhí)行來(lái)實(shí)現(xiàn)對(duì)應(yīng)語(yǔ)義.

2.3 性能優(yōu)化

對(duì)應(yīng)于上述流程圖的實(shí)現(xiàn)并不一定立即獲得顯著的加速比.GPU上程序的性能優(yōu)化需要更多地關(guān)注硬件底層.通常認(rèn)為會(huì)有基于硬件利用率、內(nèi)存以及指令三個(gè)層面的優(yōu)化方法.在進(jìn)行了基本的包括上述三個(gè)層面的常見優(yōu)化[14]后,本文著重關(guān)注非結(jié)構(gòu)網(wǎng)格下全局內(nèi)存訪問效率的提升,屬于內(nèi)存層面的深度優(yōu)化.

全局內(nèi)存訪問的優(yōu)化關(guān)鍵在于合并訪問.合并訪問指的是同一個(gè)Warp內(nèi)的Thread訪問全局內(nèi)存時(shí),內(nèi)存訪問指令可以合并成若干個(gè)32、64或者128字節(jié)的內(nèi)存讀取指令.當(dāng)相鄰Thread內(nèi)存訪問可以合并時(shí),說(shuō)明其訪問的內(nèi)存在地址上是連續(xù)的.一般在結(jié)構(gòu)網(wǎng)格下,合理安排數(shù)據(jù)形式就能實(shí)現(xiàn)合并的內(nèi)存訪問,采用 SoA(Structure of Array)而不是CPU實(shí)現(xiàn)常用的AoS(Array of Structure)形式,如守恒變量安排如圖5所示.

圖5 SoA 和 AoS 內(nèi)存組織方式對(duì)比

在非結(jié)構(gòu)網(wǎng)格中,根據(jù)本文使用的Thread映射方式,連續(xù)的Thread對(duì)應(yīng)連續(xù)的邊,但是每條邊需要訪問它相連的左右節(jié)點(diǎn),而左右節(jié)點(diǎn)由于非結(jié)構(gòu)網(wǎng)格的不規(guī)則性,在內(nèi)存上卻是不連續(xù)的.得益于GK110核心對(duì)應(yīng)的 L2 Cache 以及 Read-only Data Cache,對(duì)于內(nèi)存連續(xù)性的要求并沒有之前的硬件那么高,不再需要保證完全按照順序在地址上連續(xù),僅需保證有良好的數(shù)據(jù)局部性 (Data Locality).因此,本文采用對(duì)網(wǎng)格重排序的方法來(lái)改善基于非結(jié)構(gòu)網(wǎng)格的內(nèi)存合并訪問性能.

考慮由邊描述的節(jié)點(diǎn)連接關(guān)系的鄰接關(guān)系矩陣(nNodes×nNodes),可以通過 RCM(Reverse Cuthill-McKee algorithm)排序[14]來(lái)減小該稀疏矩陣帶寬,也意味著空間中鄰近的點(diǎn)在內(nèi)存上也是鄰近的,效果如圖6所示.

圖6 RCM排序用于減小稀疏矩陣帶寬

對(duì)網(wǎng)格點(diǎn)進(jìn)行RCM排序,然后根據(jù)排序后的網(wǎng)格點(diǎn)順序,按照每個(gè)網(wǎng)格點(diǎn)遍歷連接邊的貪婪算法進(jìn)行邊的重排序,完成排序后加速效果顯著.

RCM排序可以在線性的時(shí)間復(fù)雜度下完成,這是該方法的優(yōu)勢(shì).本文再提出另一種基于Metis網(wǎng)格分區(qū)方法[15]的排序方法用于改善數(shù)據(jù)局部性,旨在通過時(shí)間代價(jià)換取更高的加速比.Metis方法常用于并行計(jì)算過程的中的網(wǎng)格分區(qū)及負(fù)載均衡.針對(duì)GPU結(jié)構(gòu)優(yōu)化時(shí),基于Metis的排序方法較RCM排序更為直觀,語(yǔ)義上直接對(duì)應(yīng): 靠的近的點(diǎn)分進(jìn)同一組.本文通過將網(wǎng)格使用Metis分成N個(gè)小區(qū)域,每個(gè)區(qū)域包含100到1000個(gè)點(diǎn)(不同點(diǎn)的個(gè)數(shù)在這個(gè)區(qū)域內(nèi)對(duì)效果影響不明顯).如簡(jiǎn)化了的只包含四種顏色的示意圖7.

圖7 基于 Metis的網(wǎng)格排序方法示意圖

基于Metis的排序方法同樣也需要根據(jù)新的點(diǎn)順序使用貪婪算法對(duì)邊進(jìn)行重排序.這種方法在絕大多數(shù)情況下性能優(yōu)于RCM方法,但是問題在于需要占用大量的內(nèi)存并且需要較長(zhǎng)的時(shí)間.對(duì)于CPU內(nèi)存滿足條件的設(shè)備,使用這種排序算法能夠?qū)崿F(xiàn)更高的加速比.

3 算例分析

本節(jié)通過ONERA M6機(jī)翼算例驗(yàn)證求解器求解三維繞流問題的正確性.該算例是CFD中的經(jīng)典算例,幾何模型簡(jiǎn)單卻又包含了流動(dòng)的復(fù)雜性,涉及跨音速流動(dòng)、局部超音速流動(dòng)、激波以及湍流邊界層分離等現(xiàn)象.并于實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比.計(jì)算條件設(shè)置為: 來(lái)流馬赫數(shù) 0.8395,攻角 3.06°,無(wú)側(cè)滑角,雷諾數(shù)11.72×106.網(wǎng)格由四面體網(wǎng)格和附面層區(qū)域的棱柱網(wǎng)格組成,包含 210 萬(wàn)格點(diǎn),467 萬(wàn)單元,877 萬(wàn)條邊,241萬(wàn)棱柱單元,第一層網(wǎng)格高度計(jì)算得到的采用音速無(wú)量綱化后的壓力云圖與物面網(wǎng)格如圖8所示.

圖8 ONERA M6 機(jī)翼表面網(wǎng)格及壓力分布

3.1 數(shù)值結(jié)果分析

實(shí)驗(yàn)結(jié)果給出的數(shù)據(jù)為展向七個(gè)截面處壓力系數(shù)分布,以翼展b為基準(zhǔn),本節(jié)選出20%和90%翼展處截面的壓力系數(shù)分布作為對(duì)比,結(jié)果與實(shí)驗(yàn)數(shù)據(jù)基本吻合,如圖9所示.

圖9 機(jī)翼展向不同截面處壓力系數(shù)分布

3.2 加速性能分析

本節(jié)所用 CPU 設(shè)備為 Intel(R) Xeon(R) CPU E3-1230 v3,GPU 設(shè)備為 NVIDIA Tesla K40c.操作系統(tǒng)環(huán)境為 Ubuntu14.04.C++編譯器選用 gcc-4.8,CUDA 版本為7.5.GPU端以及CPU端代碼均采用-O2等級(jí)優(yōu)化.加速比通過執(zhí)行并對(duì)比算法相同的CPU串行程序和GPU并行程序得到.計(jì)算采用雙精度變量.

采用本章算例網(wǎng)格,分別在CPU和GPU上進(jìn)行1000步迭代步計(jì)算,并在GPU實(shí)現(xiàn)中結(jié)合本文提出的基于RCM及基于Metis的網(wǎng)格排序方法同樣進(jìn)行1000步迭代,所需計(jì)算時(shí)間及相應(yīng)加速比如表1所示.

結(jié)合基于Metis的網(wǎng)格排序的隱式算法求解過程在GPU上相較于CPU得到了27倍的加速比.表1也反應(yīng)出了本文提出的網(wǎng)格排序算法的高效性,分別得到了63%和69%的性能提升.如果原始網(wǎng)格過于規(guī)整,該提升也將會(huì)因?yàn)閷?duì)合并訪問改善幅度的下降而有所下降.此外,通過不斷改善網(wǎng)格質(zhì)量,或者使用結(jié)構(gòu)網(wǎng)格,本文實(shí)現(xiàn)的GPU并行程序可以實(shí)現(xiàn)50倍的加速比.可見,對(duì)內(nèi)存合并訪問的改善仍然有很大空間.另外,測(cè)試所用網(wǎng)格應(yīng)避免過小以防GPU未被充分利用,這一點(diǎn)不再展開.

表1 計(jì)算時(shí)間及相應(yīng)加速比

4 結(jié)論

針對(duì)非結(jié)構(gòu)網(wǎng)格格點(diǎn)格式在GPU上計(jì)算加速效果不佳的問題,提出了采用RCM以及基于Metis網(wǎng)格分區(qū)的網(wǎng)格重排序方法,大大地改善了基于非結(jié)構(gòu)網(wǎng)格的隱式算法在GPU上的并行加速效果.針對(duì)文中算例,兩種方法分別得到了63%和69%的加速效果.基于Metis的排序方法在效果上往往優(yōu)于RCM方法,但對(duì)于規(guī)模較大的網(wǎng)格,其耗時(shí)較長(zhǎng).采用了基于Metis的網(wǎng)格重排序方法后,文中算例結(jié)果得到了相較于CPU串行程序27倍的加速比,充分證明了本文提出的方法的可行性及優(yōu)越性.此外,本文的結(jié)果也意味著基于非結(jié)構(gòu)網(wǎng)格格點(diǎn)格式的隱式算法在GPU上也能得到良好的加速性能,相關(guān)應(yīng)用擁有著廣闊的前景.

致謝: 英偉達(dá)公司為本研究工作的展開贊助了Tesla K40計(jì)算卡,我們?yōu)榇吮硎居芍缘母兄x.

參考文獻(xiàn)

1Hagen TR,Lie KA,Natvig JR.Solving the Euler equations on graphics processing units.Proceedings of the 6th International Conference Computational Science-ICCS 2006.Berlin,Heidelberg,Germany.2006.220–227.

2Cohen J,Molemaker MJ.A fast double precision CFD code using CUDA.Proceedings of Parallel Computational Fluid Dynamics: Recent Advances and Future Directions.Lancaster,UK.2009.414–429.

3Zimmerman BJ, Wie B. Graphics-processing-unitaccelerated multiphase computational tool for asteroid fragmentation/pulverization simulation.AIAA Journal,2017,55(2): 599–609.

4Croaker P,Kessissoglou N.Aeroacoustic scattering using a particle accelerated computational fluid dynamics/boundary element technique.AIAA Journal,2016,54(7): 2116–2133.

5Mostafazadeh B,Marti F,Pourghassemi B,et al.Unsteady navier-stokes computations on GPU architectures.Proceedings of the 23rd AIAA Computational Fluid Dynamics Conference.Denver,CD,USA.2017.4508.

6Jude D,Baeder J.Extending a three-dimensional GPU RANS solver for unsteady grid motion and free-wake coupling.Proceedings of the 54th AIAA Aerospace Sciences Meeting.San Diego,CA,USA.2016.1811.

7董廷星,李新亮,李森,等.GPU 上計(jì)算流體力學(xué)的加速.計(jì)算機(jī)系統(tǒng)應(yīng)用,2011,20(1): 104–109.

8Watkins J,Romero J,Jameson A.Multi-GPU,implicit time stepping for high-order methods on unstructured grids.Proceedings of the 46th AIAA Fluid Dynamics Conference.Washington,DC,USA.2016.3956.

9劉楓,李樺,田正雨,等.基于 MPI + CUDA 的異構(gòu)并行可壓縮流求解器.國(guó)防科技大學(xué)學(xué)報(bào),2014,36(1): 6–10.

10陳龍.基于CFD/CSD耦合的旋翼氣動(dòng)彈性數(shù)值模擬[博士學(xué)位論文].南京: 南京航空航天大學(xué),2011.

11Jameson A,Yoon S.Lower-upper implicit schemes with multiple grids for the Euler equations.AIAA Journal,1987,25(7): 929–935.

12Sharov D,Nakahashi K.Reordering of 3-D hybrid unstructured grids for vectorized LU-SGS Navier-Stokes computations.Proceedings of the 13th Computational Fluid Dynamics Conference.Snowmass Village,CO,USA.1997.

13Spalart PR,Allmaras SR.A one-equation turbulence model for aerodynamic flows.Proceedings of the 30th Aerospace Sciences Meeting and Exhibit.Reno,NV,USA.1992.

14George A,Liu JW.Computer solution of large sparse positive definite.Biometrical Journal,1981,26(2): 218–225.

15Karypis G,Kumar V.Metis-unstructured graph partitioning and sparse matrix ordering system,version 5.1.0.University of Minnesota,Department of Computer Science,2013.