RKDG有限元GPU算法及其重排加速技術

高緩欽,陳紅全,張加樂,賈雪松

(1.南京航空航天大學 航空學院,南京 210016;2.非定常空氣動力學與流動控制工信部重點實驗室(南京航空航天大學),南京 210016)

間斷Galerkin(Discontinuous Galerkin method,DG)有限元方法源于20世紀70年代求解中子輸運方程[1]。后期流行的Runge-Kutta DG(RKDG)有限元方法[2-3],是結合求解非線性守恒律方程(組)提出的。該方法的主要特點是能夠構造任意階精度的格式(基于非結構網格等),但相較于傳統的低階精度格式,同一網格上控制方程的求解,因高階近似而涉及的未知量更多,消耗的計算量也就更多,影響到實際復雜工程問題的推廣應用[4-5]。因此,很有必要提高算法的計算效率。實現算法的并行化是可考慮的可行途徑。

現代圖形處理器(GPU)通常擁有數以百計或千計的計算核,浮點運算能力強勁,適合并行地處理數據密集且規模大的運算[6-8]。DG有限元方法具有很好的局部緊致性,構造的高階格式不需要非常寬的模板(stencil),因此,十分適合算法的GPU并行化[9]。早在2009年,Kl?ckner等[9]就把Nodal DG有限元方法GPU化,求解了三維Maxwell方程。之后,Siebenborn等[10]又發展用于求解三維Euler方程。Karakus等[11]則拓展用于求解二維不可壓Navier Stokes方程。與上述基于Nodal DG有限元方法不同,Xia等[12]基于Modal DG有限元方法,實現了三維Euler方程GPU并行加速。Fuhry等[13]則結合求解二維Euler方程,依據單元或單元邊界構建所需的線程結構,開展了算法GPU并行化的深入研究,其并行化效率能與Nodal DG GPU算法求解線性問題[9]相當。Gao等[14]進一步結合Navier Stokes方程,發展用于求解二維層流流動。不難看出,上述工作大多關注的是數據存儲、線程結構和核函數構建等具體的GPU化問題,尚未涉及網格拓撲存在的單元不規則分布(非結構等)對并行效率的影響問題。就低階精度格式的GPU并行化而言,這種影響不僅存在且是負面的,可采用單元[15]或點重排[16-17]等技術加以改善。因此有理由結合高階精度DG格式,對此類問題加以深入研究。

本文結合求解Navier Stokes方程,選用Modal DG有限元方法,先構建出對應的RKDG算法,然后移植到GPU架構,發展出GPU加速的RKDG并行算法。該算法在結構網格上涉及的數據依賴區結構有序,但在非結構網格上,數據依賴區雜亂無序,影響到數據訪存效率。為此,本文針對性地提出一種適合高階DG有限元算法架構(任意階)的單元分層排序方法,致力于改善并行效率。具體基于初始網格拓撲,創建單元或單元邊界對應的分層結構,逐層重排,匯總形成適合GPU對齊合并訪問(coalesced access)的數據存儲結構。結合RKDG GPU算法的算例驗證,展示出算法本身的GPU加速效果。同時,分析了網格拓撲對并行效率的影響,算例證實所提排序方法有助于進一步提升GPU并行化效率。

1 RKDG有限元并行算法

1.1 流動控制方程

二維守恒形式可壓縮Navier-Stokes方程為

(1)

其中守恒變量U、無黏通量Fi=(fi,gi)和黏性通量Fv=(fv,gv),可分別寫為:

式中:密度ρ、單位質量總能E、壓強p與沿坐標x和y方向的速度分量u和v,滿足如下狀態方程:

(2)

式中γ為流體的比熱比,對于空氣而言,γ=1.4。黏性通量項中,即:

(3)

式中:i、j分別為x或y方向,δij為Kronecker德爾塔函數,T為溫度,cP為質量定壓熱容,Pr為普朗特數,μ為黏性系數,可由Sutherland公式計算確定。

1.2 間斷Galerkin有限元離散

采用間斷Galerkin有限元方法,將計算域Ω劃分為Nc個互不重疊單元Ωi的集合。定義有限元空間:

(4)

(5)

(6)

(7)

有關re和η的具體求解過程詳見文獻[18]。單元上的近似解Uh可表示為

(8)

(9)

(10)

可以看出,上述構建的RKDG算法,相較于傳統的低階方法,相同網格上待求的未知量與Np成正比,耗費的計算量往往難以承受[21]。因此有必要對算法進行并行化加速。

1.3 基于GPU的RKDG算法

采用NVIDIA公司發布的CUDA C編程語言,對RKDG算法進行GPU并行化。這一過程也可理解為算法的CPU程序移植到GPU上。GPU內嵌平臺大多擁有數以百計乃至千計的CUDA核(core),適合并行地處理數據密集且規模大的運算,但在涉及邏輯判斷和分支結構等的計算時表現不佳[6]。為此,可結合算法任務歸類,與CPU協同運行。如圖1所示,本文把前、后處理等保留在CPU,而把最耗時的部分(迭代運算等)移植到GPU。一旦GPU上指派的任務計算完成,就再把結果傳回CPU。可以看出,每一步迭代運算由計算時間步、積分點插值、解系數更新等子任務構成。每一個子任務對應一個由用戶開發的核函數(kernel)。要實現算法的GPU化,就必須對各類子任務創建核函數。可以看到,就RKDG方法而言,待創建的核函數主要與高斯積分相關。

圖1 算法的GPU并行化架構

1.3.1 核函數的構建說明

在GPU上運行一個并行任務(調用一個核函數),是基于Grid-Block-Thread線程結構進行的。所有的線程(thread)被組織成若干個線程塊(block),再由Block匯聚成對應的線程網格(grid)。在線程結構中,塊中線程號(threadIdx.x)、線程塊號(blockIdx.x)是內置給定的,而線程塊維度(blockDim.x)可由用戶調節確定[6]。此外,內嵌的各類存儲器,訪存速度不一,其合理利用問題也是核函數創建時所必須考慮的。按訪存速度的快慢排列,GPU上的存儲器主要有寄存器、常數存儲器、共享存儲器和全局存儲器。全局存儲器容量最大,適合存放算法產生的主要數據。共享存儲器適合存儲線程塊中線程之間共享的數據,但容量相對不大。具有緩存加速等特點的常數存儲器容量十分有限,適合存放規模不大且使用頻繁的常量。寄存器在線程塊中的數目有限,但訪存速度最快,是需要考慮最大限度利用的對象。

如前述,DG方法主要涉及單元和邊界等積分計算。本文根據DG方法局部緊致的特點,沿用Fuhry等[13]的做法,對每一計算子任務,依據單元或單元邊界構建所需的線程結構,使得任一單元積分(或邊界積分)都有對應的一個線程負責其積分運算。依據這一原則,本文已對算法涉及的相關積分運算構建了對應的核函數。以殘值計算為例(見式(6)右端項),本文構建了單元積分核函數k_surf和單元邊界積分核函數k_line。k_surf僅依賴于單元內部信息,積分運算可簡單直接。比較而言,k_line涉及關聯單元,相對復雜。這里結合k_line的偽代碼,就核函數構建給出實例說明。

圖2 k_line核函數偽代碼

對于算法涉及的諸如龍格-庫塔時間推進等其他計算子任務,創建對應的核函數相對簡單,但為了節省GPU內存,可從算法層面壓縮物理量的存儲,個別物理量可考慮用到時直接計算不再存儲,本文不再一一贅述。可以看出,與線程結構相關的核函數,其創建還必須注意數據結構與存儲問題,這些因素與GPU訪存效率密切相關。

1.3.2 數據結構與存儲問題

如前述,本文算法涉及的線程結構是依據單元或單元邊界構建的,線程運行相對獨立,因此上述開發的核函數僅涉及全局存儲器、常數存儲器和寄存器。具體來說,常數存儲器用于存放參考單元基函數相關的物理量,寄存器用于處理循環重復利用(訪存頻率高)的數據(如k_line中的殘值項處理),而全局存儲器則用于存放計算產生的大容量數據。必須指出,減少算法對最慢的全局存儲器的訪存次數,有助于提升整體計算效率。通常可通過合并訪問優化大容量數據的訪存效率[22]。要做到這一點,就要結合算法構建合適的數據結構。為盡可能地滿足合并訪問的優化要求,本文把通常的多維數組映射到連續的一維數組,要求相同物理特征的變量依單元或單元邊界編號連續存放。不失一般性,如果算法涉及的物理量為如下二維數組:

(11)

式中:m為變量維度(如守恒變量等),N為單元數或單元邊界數。把二維數組a映射到一維數組A,可直接表示為

(12)

對于高階DG算法特有的總體對角塊質量矩陣(見式(9)),每一對角塊相當于一個單元上的二維數組(Np×Np),可按上述原則,先把所有對角塊映射成一維數組,再按單元順序排列形成最終的整體一維數組。可以看出,這樣的數據結構是依賴于網格拓撲中單元或邊排列順序的,還不能確保算法依賴區信息的連續提取,進而影響到所謂的GPU合并訪問。

2 重排加速技術

如果并發的32個線程(線程束[22])訪存的數據處在同一數據段,則只需要一次數據傳輸即可完成訪存,否則需要多次傳輸。研究表明,訪存效率一般與傳輸次數成反比[17]。將線程束需訪問的數據集中置于盡可能少的數據段,可有效減少數據傳輸次數。可以看出,上述算法GPU化構建的核函數,單元相關的積分運算,是依賴于自身及其關聯單元(相鄰共邊單元)的。對于依據單元編號構建的數據結構,提高訪存效率的關鍵在于同一線程束訪存單元的編號應盡可能接近或集中。這一點在復雜計算域非結構網格生成時,通常是難以顧及的。可以預見,在結構網格上,因單元排序的結構化,數據依賴單元相對集中,有助于線程束所需數據合并訪問,訪存效率會相對較高。

本文面向實用性強的非結構網格,借鑒無網格點云分層排序思想[17],針對性地提出一種單元分層排序方法,力求獲得層化結構,提升計算效率。具體基于任意初始網格單元拓撲,創建單元或單元邊界對應的分層結構,逐層重排,匯總形成適合GPU對齊合并訪問的單元排序。方法的具體實施過程可從圖3、4排序實例看出。

圖3 排序方法實施實例

圖4 排序前、后單元編號順序

上述重排加速技術是基于任意非結構網格提出的,可作為一般網格生成后處理應用。本文將結合Delaunay網格生成,給出排序算法加速效果算例分析。

3 算例與分析

表1 i5-3450 CPU和GTX TITAN GPU的參數

3.1 庫埃特流動

選用經典庫埃特模型化流動[23]對發展的算法進行考核計算。該流動是兩塊無限長平行平板間的層流流動,下面的平板水平放置且固定不動,上面的平板則以速度U水平向右運動。假定黏性系數μ為常數,板間距為2,則該流動有如下形式的解析解:

(13)

沿用文獻[23]的參數設置,U=1,下壁面溫度T0=0.80,上壁面溫度T1=0.85和馬赫數Ma1=0.1,黏性系數μ=0.001時的雷諾數Re=100。用結構或非結構網格離散計算域(0≤x≤4,0≤y≤2),并用精確的解析解施加所需的邊界條件。為了驗證算法計算精度所能達到的階數,兩種計算網格均按比例進行了剖分加密(如圖5所示),并用密度的絕對誤差(數值解與解析解的差值)的L2范數指示計算誤差,進行了不同階數算法的測試。測試結果已在表2中列出。可從表2中看出,數值解逼近的階數在密網格上更符合預期。比較而言,結構網格優于非結構網格,高階格式如預期比低階格式更精確(表現為L2誤差更小)。

表2 RKDG GPU算法不同階精度測試結果

圖5 漸次加密的結構((a)～(c))與非結構((d)～(f))網格

接著,基于結構與非結構網格,構建了非結構區域占比約為0%、33%、67%和100%的一組網格(對應圖6中網格G到J),用于測試網格不規則程度對GPU加速的影響。注意,為排除網格單元類型的影響,結構和非結構區域均統一用三角形單元離散。表3列出了這一測試結果。從表3中可看出,發展的不同階精度算法都能實現期望的GPU加速,加速比最高可達58.46;加速比大小是與非結構區域占比相關的,占比越高則加速比越小,體現出網格不規則程度對GPU加速的不利影響。

表3 算法不同階精度加速特性測試結果

圖6 網格構建示意

3.2 NACA0012翼型亞聲速黏性繞流

NACA0012翼型亞聲速黏性繞流[24]經常被用來考核發展的高階算法。本文先用2階、3階和4階算法對該繞流進行了數值模擬。計算網格由2 240個四邊形單元構成[24],翼型表面僅有41個點(如圖7所示)。計算來流馬赫數為0.5,攻角為1°,雷諾數為5 000。隨同文獻計算結果[25],表4給出了計算得到的阻力系數。可以注意到,隨著算法階數的提高,得到的阻力系數能與文獻5階結果接近,對應的馬赫數等值線也變得清晰光滑(如圖8所示)。

表4 阻力系數計算結果

圖7 NACA0012翼型計算網格

圖8 馬赫數云圖(2階與4階)

接著,本文依據如圖9所示的初始網格,構建了一套漸次加密的結構網格網格規模依次為170×30、340×60、680×120和1 360×240,進行了算法加速性能的測試。大體上,在同一階數下,加速比隨著網格規模增大而提高(表5第6列),最大加速比可達67.47。表5同時列出了重排加速比,如預期,接近于1且幾乎不隨網格變化,表明網格的結構化已很適合GPU對齊合并訪問,無需進一步重排加速。

表5 不同階算法NACA0012翼型繞流加速特性測試結果

圖9 初始結構網格(170×30)

3.3 圓柱繞流

圓柱繞流卡門渦街[26]常用來考核高階算法的低耗散性。計算來流馬赫數為0.2,雷諾數為100。計算網格由整體非結構網格(20 376個三角形單元)和局部結構網格(圓柱附近,1 560個四邊形單元)混合構成(如圖10所示)。本文先用不同階數的算法對該繞流進行了數值模擬。總體上,階數不同的算法都能捕捉到卡門渦街,階數越高則渦結構越清晰,表現為渦量等值線越光滑(如圖11所示)。圖12給出了對應的升、阻力系數周期解,表6則列出了對應的斯特哈爾數(Strouhal)。可從中看出,對應3階和4階的計算斯特哈爾數能與文獻[26]計算和實驗[27]結果接近。

表6 斯特哈爾數

圖10 圓柱繞流計算網格

圖11 圓柱繞流渦量等值線

圖12 計算得到的升、阻力系數周期解

接著,對圖10所示的網格進行粗化和細化操作,由此形成一套漸次加密的結構與非結構混合網格(單元數分別為5 300、21 936和79 428),進行算法加速特性測試。圖13給出了本算例單元關聯點陣圖,圖13中縱、橫向坐標都為單元編號,點代表兩單元關聯。從圖13中可以看出,本例初始網格對應點陣無序散布,重排后,關聯單元編號已十分接近,表現為點陣向主對角集中,呈清晰的條帶狀,展示出所提重排算法對網格排序的影響。大體上,在同一階數下,所有加速比隨著網格規模增大而提高,GPU加速比最大可達45.68(表7第6列),最大重排加速比達1.28(表7第7列),實現了近30%的進一步GPU加速。表明整體非結構化的混合網格,是與上述結構網格算例不同的,有必要進一步結合重排加速。

表7 圓柱繞流不同階數算法GPU加速特性測試結果

圖13 結構與非結構混合網格單元關聯點陣圖

3.4 NLR7301多段翼型亞聲速繞流

為了展示發展的算法處理復雜氣動外形繞流的能力,這里給出了NLR7301多段翼型[28]亞聲速無黏繞流算例,來流條件為馬赫數0.187和攻角6°。計算網格為一套漸次加密的非結構網格,單元數依次分別為5 030、19 850、78 370和314 704。兩段翼型和遠場邊界在初始網格(如圖14所示)中,表面網格點數相同,均取為65個點。本文先基于初始網格,用發展的3階RKDG GPU算法對該繞流問題進行了數值模擬。計算得到的壓力云圖和表面壓強系數分布已在圖15中給出,圖中同時給出了文獻[28]計算結果和實驗結果[29]供比較。接著,以主翼、襟翼和遠場邊界作為重排算法的起始邊界,進行了不同分層方向對重排加速比的影響測試,大體上影響不大(見表8)。

表8 不同起始邊界對重排加速比的影響

圖14 NLR7301翼型繞流的計算網格

圖15 NLR7301翼型繞流的計算結果

不失一般性,表9列出了以主翼作為起始邊界的測試結果,表中結果再次證實,對于非結構網格,GPU加速比和重排加速比都隨著網格規模增大而提高,這在一定程度上表明,對于涉及大規模非結構網格的數值模擬問題,算法的GPU化和重排加速是很有必要的。

4 結 論

1)算法涉及的線程結構是依據單元和單元邊界構建的,使得階數不同的算法可采用統一的程序結構,線程資源能得以有效利用。

2)算法結構網格上結構有序的數據依賴區,已能較好滿足GPU對齊合并訪問的要求,無需重排加速。

3)算法在非結構等混合網格上,數據依賴區非結構,重排加速是很有必要的。重排后的網格,層化結構明顯,有助于GPU對齊合并訪問,使得算法贏得GPU加速基礎上的進一步重排加速。

4)網格規模大的數值模擬問題更能贏得大的GPU加速。

5)算法是依據任意類型有限元單元設計的,適合涉及非結構等各類混合網格的復雜氣動外形數值模擬問題。