基于非結構網(wǎng)格隱式算法的GPU加速研究①

陳 龍,徐添豪,田書玲

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

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

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

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

1 數(shù)值方法

慣性系下,忽略源項的RANS方程形式如下:

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

其中,U,L和D分別為嚴格上三角矩陣,嚴格下三角矩陣和對角矩陣.

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

向前掃掠:

向后掃掠:

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

圖1 結構網(wǎng)格上的超平面形式

圖2 非結構網(wǎng)格上的超平面形式

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

2 算法 GPU 并行加速

2.1 GPU并行方法

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

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

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

可以寫成對應的由CPU調用、GPU執(zhí)行的函數(shù)(Kernel函數(shù))的形式

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

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

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

圖4 CPU 及 GPU 算法流程圖

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

2.3 性能優(yōu)化

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

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

圖5 SoA 和 AoS 內存組織方式對比

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

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

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

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

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

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

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

3 算例分析

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

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

3.1 數(shù)值結果分析

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

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

3.2 加速性能分析

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

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

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

表1 計算時間及相應加速比

4 結論

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

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

參考文獻

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 上計算流體力學的加速.計算機系統(tǒ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 的異構并行可壓縮流求解器.國防科技大學學報,2014,36(1): 6–10.

10陳龍.基于CFD/CSD耦合的旋翼氣動彈性數(shù)值模擬[博士學位論文].南京: 南京航空航天大學,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.