高緩欽，陳紅全，*，賈雪松，徐圣冠，2

(1. 南京航空航天大學(xué) 航空學(xué)院，非定?？諝鈩恿W(xué)與流動控制工信部重點實驗室，南京 210016；2. 南京工業(yè)大學(xué) 機械與動力工程學(xué)院，南京 211816)

0 引言

間斷伽遼金（discontinuous Galerkin, DG）有限元方法源于求解中子輸運方程[1]，在推廣到求解非線性守恒方程（組）[2]之后才開始快速發(fā)展。前期流行的主要是顯式龍格-庫塔DG（Runge-Kutta DG, RKDG）方法[3]，后期人們追求隱式加速，著名的廣義最小殘差（generalized minimum residual, GMRES）DG[4]和下上對稱高斯-賽德爾（lower-upper symmetric Gauss-Seidel, LU-SGS）DG[5]等算法開始涌現(xiàn)。DG有限元法融合了有限元與有限體積的優(yōu)點，可在單元內(nèi)基于選定的基函數(shù)直接拓展至高階，具有所謂的任意高階可構(gòu)造性，適合利用其對工程中經(jīng)常出現(xiàn)的激波干擾波系[6]和旋渦分離演變[7]等復(fù)雜流場細(xì)節(jié)進(jìn)行高精度刻畫；但該方法因高階近似而涉及的未知量更多，同一網(wǎng)格上的控制方程求解消耗的計算量也就更大，這制約著DG方法走向?qū)嶋H工程應(yīng)用[8]。因此，很有必要開展算法的并行化研究，尤其是圍繞現(xiàn)代先進(jìn)的圖形處理器( graphics processing unit, GPU )硬件，盡可能地提高計算效率。

現(xiàn)代GPU大多擁有數(shù)以百計乃至千計的計算核，浮點運算能力強勁，適合并行處理數(shù)據(jù)密集且規(guī)模大的運算[9]?？梢宰⒁獾?，DG方法計算涉及的絕大多數(shù)操作都是在局部單元中進(jìn)行，在單元中計算密集，所需模板（Stencil）也相對較為局部緊致，與GPU硬件特性匹配度高，十分適合GPU并行化[10]。GPU編程語言已從早期的低級語言發(fā)展到具有統(tǒng)一編程框架的高級語言，目前可供編程選擇的主要有OpenCL[11]、OpenACC[12]和統(tǒng)一計算設(shè)備架構(gòu)（compute unified device architecture, CUDA）[10,13]，其中CUDA最常用。Kl?ckner等運用CUDA統(tǒng)一編程模型結(jié)合求解三維Maxwell方程，開展了Nodal DG方法GPU并行化研究[14]。之后，Siebenborn等又發(fā)展了用于機翼等三維無黏繞流的數(shù)值模擬方法[13]。與之不同的是，F(xiàn)uhry等則基于Modal DG方法實現(xiàn)了二維Euler方程的GPU并行化求解[10]，其線程結(jié)構(gòu)是依據(jù)單元或單元邊界構(gòu)建的，適合非結(jié)構(gòu)等任意網(wǎng)格，靈活性強，且這一非線性問題的并行化效率[10]與Nodal DG GPU算法求解線性問題[14]相當(dāng)。Gao等[15]也依據(jù)此類線程結(jié)構(gòu)研究了用于求解二維流動問題的RKDG GPU算法。

DG隱式算法GPU并行化的最大挑戰(zhàn)是固有的數(shù)據(jù)關(guān)聯(lián)依賴性[16]。這種依賴性會不可避免地因數(shù)據(jù)訪存競爭引起線程沖突，進(jìn)而導(dǎo)致計算出錯。比較而言，為避免線程沖突，顯式格式（RKDG等）GPU并行化改造相對容易，而隱式格式實現(xiàn)起來則困難得多。本文在顯式RKDG GPU算法研究[15]的基礎(chǔ)上，致力于實現(xiàn)隱式DG算法的GPU并行化。具體圍繞求解Euler方程，用DG有限元法進(jìn)行空間離散，采用經(jīng)典的Roe格式處理數(shù)值通量，結(jié)合人工黏性方法[17]處理流場中可能存在的激波等間斷問題，并選用常用的LU-SGS隱式格式進(jìn)行時間推進(jìn)求解。為了克服傳統(tǒng)LU-SGS算法存在的數(shù)據(jù)關(guān)聯(lián)依賴問題，發(fā)展了單元著色分組技術(shù)，用于計算域網(wǎng)格單元的著色分組，確保同一顏色的組內(nèi)單元不再存在數(shù)據(jù)關(guān)聯(lián)，進(jìn)而使得LU-SGS算法時間推進(jìn)能按顏色組別依次并行，形成DG LU-SGS并行化算法。在此基礎(chǔ)上，為了把發(fā)展的算法移植到GPU上，先根據(jù)DG隱式算法的特點進(jìn)行算法任務(wù)歸類，再基于CUDA統(tǒng)一編程模型設(shè)計算法所需的核函數(shù)，并構(gòu)建對應(yīng)的線程與數(shù)據(jù)結(jié)構(gòu)，發(fā)展形成LU-SGS DG GPU并行算法。最后，結(jié)合典型流動算例對發(fā)展的算法進(jìn)行考核與性能測試，展示算法能基于隱式加速贏得進(jìn)一步的GPU加速。

1 間斷Galerkin有限元方法簡述

這里結(jié)合DG方法求解三維無黏流動問題，給出GPU并行化涉及的空間離散與LU-SGS格式等相關(guān)問題。

1.1 流動控制方程

三維無黏流動問題的控制方程為Euler方程，在笛卡爾坐標(biāo)系下可表示為：

式中：U=(ρ,ρu,ρv,ρw,e)T為守恒變量，F(xiàn)=(f,g,h)為無黏通量，

式中：u、v和w分別為x、y和z方向的速度分量；ρ和p分別為密度和壓強，與單位體積總能e滿足如下狀態(tài)方程：

式中γ為流體的比熱比，對于空氣而言，γ=1.4。

1.2 空間離散

運用DG有限元方法，將計算域Ω∈Rd劃分為Nc個互不重疊的單元Ωk的集合，定義如下有限元空間：

其中Uk是單元Ωk上的數(shù)值解。利用格林公式，對上式左端第二項進(jìn)行分部積分，得到方程求解的空間離散弱形式[5]：

其中n是 單元邊界?Ωk上的外法向，dσ代表面元?？梢钥闯觯捎谠贒G方法中單元邊界上允許存在間斷，所以邊界上無黏通量F的值是不唯一的。為此，常借鑒有限體積法中Riemann問題的處理方法，利用邊界左右單元計算值U-和U+，構(gòu)建數(shù)值唯一的通量函數(shù)H=H(U-,U+)，用以替換F[18]。H選用Roe格式[19]確定，那么式（6）可寫為：

進(jìn)一步，Uk用基函數(shù)?i(x)的線性組合[5]表示，即：

為處理流場中可能存在的激波等間斷問題，抑制激波附近的非物理振蕩，沿用文獻(xiàn)[17]附加人工耗散修正項的做法，即把方程（1）改寫為：

其中ε是人工黏性系數(shù)。結(jié)合BR2[20]格式處理梯度?U，方程（10）可離散為：

其中單元上的人工黏性系數(shù)εk通過光順指示器Sk計算確定。取守恒變量中的能量高階項反映單元光順度，定義，則

那么，εk可由下式計算確定[17]：

式中，sk=lgSk，s0=-4lgp，參數(shù)ε0=hk/p，κ為用戶自定義常數(shù)；hk為單元尺度?？梢钥闯?，εk在光滑區(qū)趨于要求的0，即無影響；而在激波等非光滑區(qū)εk則隨Sk增大而增大，起到抑制非物理振蕩的作用。

需要注意的是，式（11）中與人工黏性相關(guān)的Ha沿用BR2格式[20]的處理方式，可表示為：

最后，經(jīng)計算域單元匯總，得到：

式中：M是總體質(zhì)量矩陣，該矩陣為塊對角矩陣，每一矩陣塊對應(yīng)一個網(wǎng)格單元；為總體解系數(shù)向量；為總體殘值向量。

1.3 時間離散

可選用顯式或隱式對式（15）進(jìn)行時間離散求解。相較于顯式格式，隱式格式中推進(jìn)時間步穩(wěn)定性條件限制不嚴(yán)苛，步長大、收斂快，可贏得所謂的隱式加速。因此，人們追求基于隱式格式發(fā)展并行化算法，以求進(jìn)一步加速[16,21]。

采用隱式時間離散，式（15）可改寫為：

其中Λ為譜半徑。為穩(wěn)定計算，在初始ns迭代步采用指數(shù)增長形式的變CFL數(shù)。令CFLmin和CFLmax分別對應(yīng)最小和最大的CFL數(shù)，并引入?yún)?shù)：

那么，第n時間步的CFL數(shù)為：

本文后續(xù)計算中，CFLmin、CFLmax和ns統(tǒng)一設(shè)置為1、100和2000。

其中Ci為單元i的關(guān)聯(lián)單元（二維共邊，三維共面）編號集。假設(shè)Ci中的元素可根據(jù)大于與小于單元本身編號i，進(jìn)行LU自然分割，那么式（20）用LU-SGS算法[24]求解過程可分為前掃和后掃兩步：

前掃步：

后掃步：

本文將基于上式構(gòu)建基于GPU并行加速的隱式DG算法，為此必須先對算法進(jìn)行并行化改造。

2 LU-SGS算法并行化改造

2.1 數(shù)據(jù)關(guān)聯(lián)依賴問題

因此隱式算法的并行化改造一直是學(xué)者們所關(guān)注的。在結(jié)構(gòu)網(wǎng)格上，人們利用網(wǎng)格結(jié)構(gòu)化的特征進(jìn)行單元分組，實現(xiàn)算法的并行化改造。如SOR算法，通過結(jié)構(gòu)網(wǎng)格單元紅（Red）黑（Black）自然分組，組內(nèi)單元互不關(guān)聯(lián)，形成所謂的Red-Black SOR并行算法[25]。但對于無結(jié)構(gòu)化特征的非結(jié)構(gòu)網(wǎng)格，此類做法不能直接運用。Malte等利用四色定理[26]，將二維非結(jié)構(gòu)三角形網(wǎng)格分為組內(nèi)單元互不關(guān)聯(lián)的至多4個顏色組，實現(xiàn)了隱式LU-SGS算法的GPU并行。本文則依據(jù)LU-SGS DG算法關(guān)聯(lián)單元數(shù)據(jù)依賴的特征，發(fā)展一種面向任意網(wǎng)格的單元著色分組技術(shù)，用于式（22）和式（23）的并行化改造。

2.2 單元著色分組與LU-SGS算法并行化改造

注意計算域Ω∈Rd已離散為互不重疊的網(wǎng)格單元：。為了克服LU-SGS算法的數(shù)據(jù)關(guān)聯(lián)依賴問題，根據(jù)式（21）中的關(guān)聯(lián)單元編號集Ci，構(gòu)建辨識區(qū)分關(guān)聯(lián)特征的著色數(shù)組color(:)，使得對于?j∈Ci,color(i)≠color(j)，其中color(i)∈N+代表單元Ωi的顏色。這就是所謂的單元著色過程，圖1給出了具體的實施方法，記為算法1。

圖1 算法1：網(wǎng)格單元著色Fig. 1 Steps for grid cell coloning (Algorithm 1)

算法1中第1行為數(shù)組color(:)的初始化；第2行選定起始單元c0并著色。一旦c0選定，最終著色結(jié)果也唯一確定。選取物面邊界相連的單元作為起始單元，符合物面擾動向流場傳播的特征。測試表明，不同起始單元對計算的影響可忽略不計。為方便，本文算例起始單元均選為編號最小的物面單元。第3～6行是著色循環(huán)體，遍歷著色所有Ci中的單元。在最后第8行，提取所用顏色數(shù)Ncolor供后續(xù)計算使用。

這里結(jié)合著色實例（見圖2），直觀展示所提單元著色分組技術(shù)。圖2（a）為結(jié)構(gòu)網(wǎng)格的著色結(jié)果，可以看出，對于二維結(jié)構(gòu)網(wǎng)格，只需兩種顏色即可，其效果相當(dāng)于Red-Black分組技術(shù)[25]。圖2（b）則給出了非結(jié)構(gòu)網(wǎng)格的著色結(jié)果，其網(wǎng)格的非結(jié)構(gòu)性要求使用更多的顏色。本例圖中為4種顏色，分別以1～4標(biāo)記。

圖2 二維網(wǎng)格著色示例Fig. 2 Examples of two-dimensional mesh coloring

在網(wǎng)格單元著色后，再按著色總數(shù)Ncolor構(gòu)建如下顏色組：

由于組內(nèi)單元已互不關(guān)聯(lián)，可以按組別依次并行。圖3列出了算法2基于單元著色分組的LUSGS并行過程?？梢钥闯?，并行化改造后的LUSGS算法關(guān)聯(lián)單元編號集Ci是以顏色標(biāo)記的大小進(jìn)行LU分割，從第一個顏色到最后一個顏色完成算法的前掃步，后掃步則反過來，即從最后一個顏色循環(huán)到第一個顏色。需要注意的是，傳統(tǒng)的LU-SGS算法（見式（22）和式（23））是依單元串行計算，改造后的算法是以顏色組群依次并行。下一節(jié)將進(jìn)一步討論該算法的GPU化，即最大程度地把算法移植到高性能的GPU架構(gòu)上，贏得GPU并行化加速。

圖3 算法2：LU-SGS算法并行過程Fig. 3 LU-SGS algorithm parallel preocess (Algorithm 2)

3 DG有限元GPU隱式并行算法

3.1 算法的整體框架

算法的GPU并行化就是算法從CPU到GPU的移植過程。為實現(xiàn)高效移植、贏得GPU加速，有必要兼顧算法特點與GPU硬件特征。目前，GPU硬件大多內(nèi)嵌幾百乃至數(shù)千的計算核，適合并行處理規(guī)模大且計算密集的運算，但在涉及邏輯判斷和分支結(jié)構(gòu)等計算時表現(xiàn)不佳[15]。為此，對算法2進(jìn)行任務(wù)歸類，并與CPU協(xié)同運行。如圖4所示，前、后處理等相關(guān)部分保留在CPU上執(zhí)行，時間推進(jìn)迭代部分最耗時，因此移植到GPU上計算。一旦指派到GPU上的各項任務(wù)計算完成，再把計算結(jié)果傳回CPU?？梢钥闯?，每一時間推進(jìn)步由殘值計算、全局系統(tǒng)矩陣計算、前掃和后掃等子任務(wù)構(gòu)成，每一個子任務(wù)則對應(yīng)一個GPU核函數(shù)（Kernel）。因此，從某種意義上，算法的GPU化就是對各類子任務(wù)創(chuàng)建對應(yīng)的核函數(shù)。本文選用廣泛使用的CUDA C編程語言構(gòu)建核函數(shù)。必須指出的是，就LU-SGS DG隱式算法而言，全局系統(tǒng)矩陣更新、算法中的前掃步和后掃步運算等相關(guān)核函數(shù)是最為至關(guān)重要的，必須結(jié)合算法的特點和GPU訪存等硬件特征來構(gòu)建。

圖4 LU-SGS DG算法GPU并行架構(gòu)Fig. 4 General procedure of GPU-based LU-SGS DG algorithm

3.2 核函數(shù)構(gòu)建問題

核函數(shù)的構(gòu)建與GPU上各類存儲器的利用和線程結(jié)構(gòu)的創(chuàng)建密切相關(guān)。計算所需的線程（Thread）被集合為若干個線程塊（Block），再由Block匯集成線程網(wǎng)格（Grid），形成所謂的Grid-Block-Thread線程結(jié)構(gòu)。Block中的線程號（threadIdx.x）、Block號（blockIdx.x）以及自定義的Block維度（blockDim.x）等內(nèi)置變量，用來計算確定Grid中的每個具體的線程編號[27]。依據(jù)上述創(chuàng)建的線程結(jié)構(gòu)實現(xiàn)核函數(shù)的調(diào)用，其調(diào)用效率還與GPU上各類存儲器（訪存速度與存儲空間大小不一）的合理利用相關(guān)。這些存儲器按存儲空間從大到小排列，主要有全局存儲器、寄存器、共享存儲器和常數(shù)存儲器。容量最大的全局存儲器適合存放計算產(chǎn)生的主要數(shù)據(jù)，但訪存速度是所列存儲器中最慢的，訪存數(shù)據(jù)時要盡可能地利用對齊合并訪問[27]，減少訪存次數(shù)。訪存速度最快的寄存器，數(shù)目相對有限，需要最大限度使用。共享存儲器可用于線程塊中線程間數(shù)據(jù)共享，但容量相對不大。具有緩存加速等特點的常數(shù)存儲器容量十分有限，往往存放頻繁使用且規(guī)模較小的常量。

如前所述，系統(tǒng)矩陣的迭代更新、前掃和后掃涉及的相關(guān)運算都是隱式DG算法所特有的，主要涉及單元和單元邊界等積分計算（見式（11））。為了盡可能地滿足對齊合并訪問的要求，本文依單元設(shè)計創(chuàng)建積分運算核函數(shù)，并借鑒Fuhry等[10]的做法，依據(jù)單元構(gòu)建相關(guān)核函數(shù)所需的線程結(jié)構(gòu)，使得線程與單元相關(guān)的運算一一對應(yīng)。依據(jù)這一原則，本文已對算法涉及的各子任務(wù)構(gòu)建了對應(yīng)的核函數(shù)。其中，最耗時的時間推進(jìn)迭代相關(guān)的核函數(shù)已在圖5中列出，并對代表性的核函數(shù)k_Diagvol和k_Fsweep給出具體說明。

圖5 LU-SGS DG GPU算法時間推進(jìn)調(diào)用過程Fig. 5 Time marching of LU-SGS DG GPU algorithm

核函數(shù)k_Diagvol用于計算單元積分對于全局系統(tǒng)矩陣中對角塊Aii（見式（21））的貢獻(xiàn)，其偽代碼和必要的說明已在圖6中列出。圖中前兩行計算確定線程號id并界定其范圍；接著兩行進(jìn)行初始化；從第5行到第16行為核函數(shù)循環(huán)運算的主體，循環(huán)遍歷所有單元積分點，并計算累加各積分點的貢獻(xiàn)來更新對角塊。

圖6 核函數(shù)k_Diagvol的偽代碼Fig. 6 Pseudocode of the kernel k_Diagvol

圖7則給出了前掃核函數(shù)k_Fsweep 的偽代碼和必要說明。該核函數(shù)用于完成算法前掃步運算，以顏色組群依次并行。圖中前兩行完成線程號id的計算與范圍界定，涉及的Ns為顏色組Bs中的單元總數(shù)。接著為遍歷單元所有邊界的外循環(huán)（第4～17行），先根據(jù)關(guān)聯(lián)單元顏色標(biāo)記LU分割，確定需要計算的單元邊界（第5行），再進(jìn)入邊界積分點內(nèi)循環(huán)，遍歷積分點計算累加更新殘值；跳出外循環(huán)體后，計算更新解系數(shù)未知量（第20行）。

圖7 核函數(shù)k_Fsweep的偽代碼Fig. 7 Pseudocode of the kernel k_Fsweep

算法涉及的后掃核函數(shù)k_Bsweep也可以類似構(gòu)建。至于殘值、時間步計算等相關(guān)核函數(shù)，顯式隱式類似，可參照作者先前有關(guān)顯式GPU算法的工作[15]，這里不再一一贅述。不難發(fā)現(xiàn)，上述核函數(shù)的創(chuàng)建兼顧了數(shù)據(jù)結(jié)構(gòu)構(gòu)建與各類存儲器利用，這些因素與GPU訪存效率密切相關(guān)。

3.3 存儲器利用與數(shù)據(jù)結(jié)構(gòu)問題

如前所述，本文算法涉及的線程結(jié)構(gòu)是利用了DG算法依單元計算密集的特征，依據(jù)單元或單元邊界構(gòu)建，線程運行相對獨立，因此上述開發(fā)的核函數(shù)僅涉及全局存儲器、常數(shù)存儲器和寄存器。必須明確的是，為了提高核函數(shù)的并行效率，要盡可能地減少對慢速全局存儲器的訪存次數(shù)，最大限度地利用資源有限的快速寄存器。本文依據(jù)這一原則，常數(shù)存儲器用于存放參考單元基函數(shù)相關(guān)的物理量，寄存器用于處理多次重復(fù)利用的數(shù)據(jù)（如k_Diagvol中的基函數(shù)梯度處理），而全局存儲器則用于存放計算產(chǎn)生的大容量數(shù)據(jù)。但相比于CPU內(nèi)存，全局存儲器容量仍相對有限。因此，在設(shè)計核函數(shù)前，對運算所需數(shù)據(jù)進(jìn)行了分析篩選，尤其是系統(tǒng)矩陣（見式（21））相關(guān)的大規(guī)模數(shù)據(jù)。本文只計算并存儲了前掃和后掃中都需使用的對角塊逆矩陣；將非對角塊Aij相關(guān)計算嵌入前掃或后掃核函數(shù)（見圖7），不再進(jìn)行存儲，以減少對慢速全局存儲器的訪存。

數(shù)據(jù)存放與調(diào)用涉及的數(shù)據(jù)結(jié)構(gòu)，要盡可能地滿足合并訪問的優(yōu)化要求[16]。這一點作者已結(jié)合RKDG GPU算法的工作進(jìn)行了研究，發(fā)現(xiàn)將多維數(shù)組映射為連續(xù)的一維數(shù)組，使得具有相同物理特征的變量依單元或單元邊界編號連續(xù)存放，有助于滿足數(shù)據(jù)的合并訪問，進(jìn)而提高訪存效率。本文發(fā)展的隱式算法沿用了這一做法，具體可參見文獻(xiàn)[15]。

4 算例考核與性能測試

先選用二維單段/多段翼型和三維機翼等典型繞流，考核發(fā)展的LU-SGS DG GPU算法；再結(jié)合考核算例，進(jìn)行GPU加速性能的測試與分析。

4.1 NACA 0012翼型跨聲速繞流

NACA 0012翼型跨聲速繞流算例常用來考核算法捕捉激波的能力[28]。本文對馬赫數(shù)Ma=0.8和攻角α=1.25°的來流條件進(jìn)行了數(shù)值模擬，計算得到的馬赫數(shù)云圖和物面壓強系數(shù)分布見圖8。從圖8（a）的馬赫數(shù)云圖中看出，計算采用了較粗的混合網(wǎng)格，僅包含185個三角形單元和1706個四邊形單元。圖8（b）同時畫出了二階（p=1）到四階（p=3）的計算結(jié)果，可從激波強度與位置注意到，由于網(wǎng)格較粗，二階結(jié)果捕捉的激波橫跨多個數(shù)據(jù)點，隨著階數(shù)的提高改善明顯，表現(xiàn)為本文四階結(jié)果與文獻(xiàn)[29]更密網(wǎng)格（256×256）的四階DG結(jié)果接近。這一點也可從表1的升阻力系數(shù)變化看出。

表1 NACA 0012翼型繞流的升阻力系數(shù)Table 1 Lift and drag coefficients of NACA 0012 airfoil

圖8 繞NACA 0012翼型跨聲速流動計算結(jié)果Fig. 8 Computational results of the flow over NACA 0012 airfoil

為了直觀地展示加速效果，以四階計算為例，圖9同時給出了顯式RKDG GPU算法和隱式LU-SGS DG GPU算法的收斂歷程。如預(yù)期，發(fā)展的隱式算法比顯式收斂速度更快（圖9（a）），計算更省時（圖9（b））。

圖9 NACA 0012翼型繞流的收斂歷程Fig. 9 Convergence histories of flow over NACA 0012 airfoil

4.2 NLR 7301多段翼型繞流

為了展示發(fā)展的算法處理復(fù)雜氣動外形的能力，本算例選用NLR 7301多段翼型亞聲速無黏繞流，來流條件為馬赫數(shù)0.185和攻角6°。計算網(wǎng)格見圖10（a），圓形遠(yuǎn)場邊界半徑為30倍弦長，共包含5274個三角形單元，兩段翼型表面網(wǎng)格點數(shù)相同，均取為65個點。在該網(wǎng)格上，采用發(fā)展的二階（p= 1）到四階（p=3）DG方法計算該繞流問題，計算得到的壓力云圖和表面壓強系數(shù)分布已在圖10中給出；圖中同時給出了實驗結(jié)果[30]和參考的SU2開源軟件[31]密網(wǎng)格（包含51932個三角形單元）的無黏結(jié)果供比較。需要注意的是，SU2采用的密網(wǎng)格，其單元數(shù)已相當(dāng)于DG所用網(wǎng)格的10倍。表2中則列出了對應(yīng)的升阻力系數(shù)CL和CD。算例表明，相較于低階SU2密網(wǎng)格結(jié)果，用四階DG計算可獲得精度相當(dāng)?shù)慕Y(jié)果，且網(wǎng)格量可粗化到十分之一左右。

表2 NLR 7301翼型繞流的升阻力系數(shù)Table 2 Lift and drag coefficients of NLR 7301 airfoil

圖10 NLR 7301翼型繞流的計算結(jié)果Fig. 10 Computational results of flows over NLR 7301 airfoil

4.3 M6機翼跨聲速繞流

選用三維M6機翼實際外形繞流問題[32]考核算法的實用性。算例來流馬赫數(shù)為0.84，攻角為3.06°。采用二階（p= 1）和三階（p= 2）DG算法，計算在僅有62897個四面體網(wǎng)格單元（相對稀疏）的網(wǎng)格上進(jìn)行（見圖11）。圖12給出了三階算法計算得到的物面壓強系數(shù)云圖和對稱面壓強系數(shù)等值線。不同算法計算的沿機翼展向非對稱截面壓強系數(shù)分布在圖13中給出。從3個典型截面捕捉的激波位置和強度來看，隨階數(shù)提升，本文三階結(jié)果總體上與實驗結(jié)果[33]和文獻(xiàn)[32]的無黏計算結(jié)果（四階DG）較為吻合，展示出發(fā)展的算法模擬實際氣動外形的能力。以三階計算為例，三維算例再次證實，相較于顯式GPU算法，發(fā)展的隱式GPU算法收斂速度的確更快（圖14（a）），計算更省時（圖14（b））。因三維算例的規(guī)模相對較大，這種特性更為顯著，本算例消耗的時間僅為對應(yīng)顯式算法的十三分之一。

圖11 M6機翼計算網(wǎng)格Fig. 11 Computational mesh of M6 wing

圖13 M6機翼典型展向截面壓強系數(shù)分布Fig. 13 Surface pressure coefficients at typical spanwise sections of M6 wing

圖14 M6機翼繞流的收斂歷程Fig. 14 Convergence histories of flow over M6 wing

4.4 GPU加速特性測試

如上所述，發(fā)展的隱式LU-SGS DG GPU算法是由CPU上的隱式LU-SGS DG算法移植到GPU上所得。為定量分析算法的GPU加速性能，設(shè)T為算法單一步時間步的計算耗時，并用下標(biāo)區(qū)分算法在CPU和GPU上的耗時，那么隱式GPU加速比可定義為TCPU/TGPU。依照傳統(tǒng)做法，T統(tǒng)一取為1000個時間步迭代耗時的平均值[9]。所有測試算例均以雙精度運行在Windows工作平臺上。該平臺配備有英特爾i5-3450 CPU和英偉達(dá)GTX TITAN GPU等硬件，硬件性能參數(shù)可參見文獻(xiàn)[15]。

4.4.1 規(guī)模效應(yīng)

研究表明，GPU加速比與構(gòu)建的線程塊Block大小密切相關(guān)[34]。作者對顯式RKDG GPU算法的測試也表明，Block中線程數(shù)取為64或128時，加速效果最好[15]。本文以64或128作為基準(zhǔn)，考察了Block中更小或更大的線程數(shù)，按從小到大排列，線程數(shù)為8、16、32、64、128、256、512或1024。具體測試工作結(jié)合二維NACA 0012翼型和三維M6機翼繞流算例進(jìn)行。同時必須注意，算例網(wǎng)格單元數(shù)應(yīng)遠(yuǎn)大于GPU計算核數(shù)，即排除所謂的規(guī)模影響[34]，確保測試可信。為此，本文將網(wǎng)格相對較粗的二維算例加密至包含48639個三角形單元，三維算例網(wǎng)格規(guī)模已足夠，保持不變。圖15給出了對應(yīng)的測試結(jié)果，圖中橫坐標(biāo)代表線程塊中的線程數(shù)，縱坐標(biāo)為計算耗時。二維和三維不同階數(shù)的測試結(jié)果都表明，當(dāng)每個Block包含64或128個線程時，計算耗時最小。這一結(jié)果與文獻(xiàn)[15]顯式算法得到的結(jié)論一致。

圖15 線程塊大小對計算耗時的影響Fig. 15 Effects of the thread block size on computational time TGPU

為了不失一般性，設(shè)置Block中的線程數(shù)為64，測試網(wǎng)格規(guī)模對加速比的影響。圖16給出了加速比隨網(wǎng)格規(guī)模變化的曲線?？梢钥闯觯篌w上，在達(dá)到GPU性能極限前，加速比均隨網(wǎng)格規(guī)模增大而增大，最大可達(dá)18.64。

圖16 網(wǎng)格規(guī)模對加速比的影響Fig. 16 Speedups with respect to different mesh sizes

4.4.2 隱式GPU算法總體加速性能

匯總4.1～4.3節(jié)中的算例，可以得到發(fā)展算法的總體加速性能。相關(guān)數(shù)據(jù)已在表3中列出?？梢园l(fā)現(xiàn)，不同階數(shù)的算例測試都能贏得預(yù)期的GPU加速，二維算例加速比最高可達(dá)11.68，三維算例加速比則高達(dá)14.15。需要說明的是，此加速比是在隱式加速的基礎(chǔ)上得到的，此外，規(guī)模較大的三維機翼算例能贏得更大的加速比，展示出算法處理此類大規(guī)模問題的潛力。

表3 隱式GPU算法的總體加速性能Table 3 GPU speedups of the developed algorithm

5 結(jié)論

本文遵循結(jié)構(gòu)網(wǎng)格紅黑著色分組思想，依據(jù)LUSGS DG算法關(guān)聯(lián)單元數(shù)據(jù)依賴的特征，發(fā)展形成了一種面向任意網(wǎng)格的單元著色分組技術(shù)，并成功地用于隱式算法的并行化改造；在此基礎(chǔ)上，基于CUDA統(tǒng)一編程模型，依據(jù)單元或單元邊界創(chuàng)建了核函數(shù)與對應(yīng)的線程結(jié)構(gòu)，成功地把并行化改造后的算法移植到了GPU，形成了隱式LU-SGS DG GPU算法。發(fā)展的算法通過算例驗證與性能測試，得到的結(jié)論主要有：

1） 所提單元著色分組方法可用于結(jié)構(gòu)與非結(jié)構(gòu)混合等任意網(wǎng)格，當(dāng)用于結(jié)構(gòu)網(wǎng)格時，該方法與傳統(tǒng)的紅黑著色法等價，因此可為其他隱式算法并行化改造提供借鑒。

2） 并行化LU-SGS DG算法是以顏色組群依次并行的，所需顏色數(shù)越少，并行化程度越高，因此，可用的任意網(wǎng)格應(yīng)盡可能地選擇結(jié)構(gòu)化程度高的，以減少顏色數(shù)。

3） 發(fā)展的算法都能在隱式加速的基礎(chǔ)上贏得進(jìn)一步的GPU加速，較顯式GPU算法更省時，且對三維等規(guī)模大的數(shù)值模擬加速效果會更好。

4） 算法是依單元面向任意網(wǎng)格構(gòu)建的，可以靈活地處理涉及復(fù)雜計算域的離散問題。

當(dāng)前算法中線程結(jié)構(gòu)主要是依據(jù)單元構(gòu)建的，下一步可結(jié)合DG算法的特征，考慮不同的線程結(jié)構(gòu)，進(jìn)一步完善本文發(fā)展的算法。