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

(1.南京航空航天大學 航空學院,南京 210016;2.非定?？諝鈩恿W與流動控制工信部重點實驗室(南京航空航天大學),南京 210016)

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

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

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

1 RKDG有限元并行算法

1.1 流動控制方程

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

(1)

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

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

(2)

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

(3)

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

1.2 間斷Galerkin有限元離散

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

(4)

(5)

(6)

(7)

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

(8)

(9)

(10)

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

1.3 基于GPU的RKDG算法

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

圖1 算法的GPU并行化架構(gòu)

1.3.1 核函數(shù)的構(gòu)建說明

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

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

圖2 k_line核函數(shù)偽代碼

對于算法涉及的諸如龍格-庫塔時間推進等其他計算子任務,創(chuàng)建對應的核函數(shù)相對簡單,但為了節(jié)省GPU內(nèi)存,可從算法層面壓縮物理量的存儲,個別物理量可考慮用到時直接計算不再存儲,本文不再一一贅述?？梢钥闯?與線程結(jié)構(gòu)相關(guān)的核函數(shù),其創(chuàng)建還必須注意數(shù)據(jù)結(jié)構(gòu)與存儲問題,這些因素與GPU訪存效率密切相關(guān)。

1.3.2 數(shù)據(jù)結(jié)構(gòu)與存儲問題

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

(11)

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

(12)

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

2 重排加速技術(shù)

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

本文面向?qū)嵱眯詮姷姆墙Y(jié)構(gòu)網(wǎng)格,借鑒無網(wǎng)格點云分層排序思想[17],針對性地提出一種單元分層排序方法,力求獲得層化結(jié)構(gòu),提升計算效率。具體基于任意初始網(wǎng)格單元拓撲,創(chuàng)建單元或單元邊界對應的分層結(jié)構(gòu),逐層重排,匯總形成適合GPU對齊合并訪問的單元排序。方法的具體實施過程可從圖3、4排序?qū)嵗闯觥?/p>

圖3 排序方法實施實例

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

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

3 算例與分析

表1 i5-3450 CPU和GTX TITAN GPU的參數(shù)

3.1 庫埃特流動

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

(13)

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

表2 RKDG GPU算法不同階精度測試結(jié)果

圖5 漸次加密的結(jié)構(gòu)((a)～(c))與非結(jié)構(gòu)((d)～(f))網(wǎng)格

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

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

圖6 網(wǎng)格構(gòu)建示意

3.2 NACA0012翼型亞聲速黏性繞流

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

表4 阻力系數(shù)計算結(jié)果

圖7 NACA0012翼型計算網(wǎng)格

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

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

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

圖9 初始結(jié)構(gòu)網(wǎng)格(170×30)

3.3 圓柱繞流

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

表6 斯特哈爾數(shù)

圖10 圓柱繞流計算網(wǎng)格

圖11 圓柱繞流渦量等值線

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

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

表7 圓柱繞流不同階數(shù)算法GPU加速特性測試結(jié)果

圖13 結(jié)構(gòu)與非結(jié)構(gòu)混合網(wǎng)格單元關(guān)聯(lián)點陣圖

3.4 NLR7301多段翼型亞聲速繞流

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

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

圖14 NLR7301翼型繞流的計算網(wǎng)格

圖15 NLR7301翼型繞流的計算結(jié)果

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

4 結(jié) 論

1)算法涉及的線程結(jié)構(gòu)是依據(jù)單元和單元邊界構(gòu)建的,使得階數(shù)不同的算法可采用統(tǒng)一的程序結(jié)構(gòu),線程資源能得以有效利用。

2)算法結(jié)構(gòu)網(wǎng)格上結(jié)構(gòu)有序的數(shù)據(jù)依賴區(qū),已能較好滿足GPU對齊合并訪問的要求,無需重排加速。

3)算法在非結(jié)構(gòu)等混合網(wǎng)格上,數(shù)據(jù)依賴區(qū)非結(jié)構(gòu),重排加速是很有必要的。重排后的網(wǎng)格,層化結(jié)構(gòu)明顯,有助于GPU對齊合并訪問,使得算法贏得GPU加速基礎(chǔ)上的進一步重排加速。

4)網(wǎng)格規(guī)模大的數(shù)值模擬問題更能贏得大的GPU加速。

5)算法是依據(jù)任意類型有限元單元設(shè)計的,適合涉及非結(jié)構(gòu)等各類混合網(wǎng)格的復雜氣動外形數(shù)值模擬問題。