一種基于塊雅可比迭代的高階FR 格式隱式方法1)

于要杰 劉 鋒 高 超 馮 毅,2)

*(西北工業(yè)大學(xué)航空學(xué)院,西安 710072)

?(加州大學(xué)爾灣分校機械與航天工程系,美國加利福尼亞州 92697-3975)

**(成都流體動力創(chuàng)新中心,成都 610071)

引言

高階精度數(shù)值格式因其高分辨率、低色散和低耗散等特點正日益受到關(guān)注[1-2].相比目前工業(yè)界最常用的二階精度格式,高階精度格式被認(rèn)為在花費相同計算代價的前提下具備獲得更高精度數(shù)值結(jié)果的巨大潛力.當(dāng)前最常見的高階精度數(shù)值格式主要包括基于結(jié)構(gòu)網(wǎng)格的高階有限差分方法[3-4]、非結(jié)構(gòu)網(wǎng)格上的有限體積法[5-6]以及基于局部單元近似的間斷有限元(discontinuous Galerkin,DG)[7-8]、譜體積(spectral volume,SV)[9]和譜差分(spectral difference,SD)[10]等方法.

最近,一種基于非結(jié)構(gòu)網(wǎng)格的高階精度通量重構(gòu)格式(flux reconstruction,FR)因其構(gòu)造簡單且更加通用而愈加流行.FR 格式最先于2007 年由Huynh[11-12]提出并證明了,至少對于一維線性方程,FR 格式通過不同修正多項式的選擇,不僅囊括了DG 和SD 這兩種格式,并且能夠產(chǎn)生新的高階格式.2009 年,Wang和Gao 等[13]提出“l(fā)ifting collocation penalty,LCP”格式并將其推廣到歐拉方程和Navier-Stokes 方程.由于FR 和LCP 這兩種格式具有高度的相似性和內(nèi)在聯(lián)系,因此后來也被Huynh 和王志堅[14]聯(lián)合命名為FR/CPR (correction procedure via reconstruction) 格式.2011 年,Vincent 等[15]提出了一種針對一維線性對流方程的能量穩(wěn)定的FR 格式,即VCJH 格式,后續(xù)又在此格式的基礎(chǔ)上開發(fā)了三維并行CFD 求解器Hi-FiLES[16]和PyFR[17].如今,FR 已不僅是一種高階格式,它還為幾種現(xiàn)有的高階格式,例如SD 和DG,提供了一個統(tǒng)一的框架.相比DG 和SD 的原有形式,FR 在概念上更加簡單且計算效率更高[18].關(guān)于FR和DG 間的聯(lián)系,Allaneau 和Jameson[19]和De Grazia等[20]已經(jīng)進行了深入研究.

盡管FR 格式已經(jīng)取得了較大的進步,但將其應(yīng)用于復(fù)雜外形流動的數(shù)值模擬時,仍然有一些難點問題需要解決,例如計算效率的問題、曲面邊界處理方法、間斷偵測方法與限制器的構(gòu)造等.其中,計算效率的問題在很大程度上制約了FR 格式在實際工程中的應(yīng)用.在當(dāng)前的FR 格式應(yīng)用中,時間離散多采用存儲需求小且易實現(xiàn)的顯式Runge-Kutta 類方法[21].但隨著格式精度的提高,FR 對應(yīng)的CFL 穩(wěn)定性條件越來越嚴(yán)格,時間推進步長會受到明顯的限制.對一些復(fù)雜流動問題(流場中存在邊界層、剪切層等小尺度流動結(jié)構(gòu))的計算更是如此.因此,發(fā)展FR 格式的加速收斂方法,具有十分重要的工程應(yīng)用價值.目前常用的加速收斂方法在算法層面主要有隱式時間推進方法、當(dāng)?shù)貢r間步長、殘差光順、多重網(wǎng)格以及物理層面的并行計算技術(shù)等.在這些方法中,隱式時間推進方法的加速效果是十分明顯的[22].相比顯式時間方法,隱式方法受CFL 穩(wěn)定性條件的約束小,可以選取更大的時間推進步長,進而能夠有效地提高CFD 求解的效率.近年來已有學(xué)者提出了多種基于FR 格式的隱式時間方法,例如LU-SGS 方法[23]、預(yù)處理GMERS 方法[24]、無矩陣(Jacobianfree) Newton Krylov 方法[25]以及多色高斯-賽德爾(multi-colored Gauss Seidel,MCGS) 方法[26]等.這些文章中以數(shù)值實驗證明,使用隱式方法的FR 格式相比顯式具有非常明顯的計算效率提升.以上這些隱式時間方法大都是在CPU 架構(gòu)上實現(xiàn),由于基于局部單元近似的FR 格式的絕大多數(shù)操作都是在局部單元中進行,計算量較為集中,因此更適合GPU 計算.當(dāng)前,基于CPU 并行策略的有限體積求解器一般只能使用并行集群性能峰值的5%～10%,而使用GPU 計算的FR 格式計算速度最高可以達到性能峰值的58%左右[27].因此,基于GPU 并行的隱式時間推進方法有望進一步提高FR 格式的計算效率.

近年來,GPU 計算在CFD 領(lǐng)域中的應(yīng)用已經(jīng)越來越普遍.已有的研究結(jié)果表明[27-29],對于相同規(guī)模的問題,GPU 并行能夠較大程度地提升高階CFD 程序的計算效率.例如,Romero 等[30]利用顯式時間方法的VCJH 格式計算非定常流動問題,發(fā)現(xiàn)GPU 并行相比基于MPI 并行的CPU 節(jié)點具有更高的計算效率.Vermeire 等[31]使用PyFR 求解器在多個GPU 卡上以三維Taylor-Green 渦和SD7003 翼型為例進行隱式大渦模擬計算,發(fā)現(xiàn)與使用CPU 節(jié)點并行的商業(yè)CFD 軟件Star-CCM+相比,在計算成本相當(dāng)?shù)那疤嵯?5 階FR 計算得到的模擬結(jié)果更精細(xì).最近,Jourdan 和Wang 等[32]在美國“Summit”超級計算機上使用FR/CPR 高階求解器,分別利用了800 個節(jié)點的CPU 和4800 個GPU 卡進行了大渦模擬計算,發(fā)現(xiàn)使用顯式時間推進方法時,在六面體和四面體這兩種網(wǎng)格單元類型上,GPU 相比CPU 并行分別有一個和兩個量級左右的加速效果.

在上述的大規(guī)模CFD 數(shù)值模擬中,基于FR 格式的GPU 計算使用的都為顯式時間方法.為了探索如何進一步提高FR 格式的求解效率,本文以二維歐拉方程為例,提出一種基于高階FR 格式的單GPU并行隱式時間求解方法.通過塊Jacobi 迭代的方式,改變了經(jīng)過隱式離散后全局線性方程組左端矩陣的特征,克服了影響求解并行性的相鄰單元依賴.最終,將解全局線性方程組轉(zhuǎn)化解一系列局部單元線性方程組,進而又通過LU 分解法并行求解這些局部線性方程組.數(shù)值實驗表明,相比顯式方法而言,該隱式方法在計算效率方面至少有一個量級的提升.并且,該隱式方法不局限于FR 格式,對DG 和SD 等基于局部單元近似的高階格式也具有參考意義.

1 控制方程

無黏流動的控制方程為歐拉方程,其二維形式如下

其中,U=[ρ,ρu,ρv,ρE]為守恒變量,F和G分別為x和y方向的無黏通量

式中,ρ 是密度,u和v分別表示x和y方向的速度,p為壓強,E代表單位質(zhì)量流體的總能量。本文中假設(shè)流體為理想氣體,則有

其中,γ=Cp/Cv為比熱容比,本文中其值為1.4.

2 高階通量重構(gòu)格式

本節(jié)將以二維標(biāo)量守恒律為例來簡要地介紹FR格式的基本思想.在區(qū)域Ω ∈R2上的一般形式的二維標(biāo)量守恒律形式為

其中,u=u(x,y,t)為標(biāo)量守恒量,f=(f,g)為通量項,f=f(u)和g=g(u)分別為f在x和y方向的分量.

下面將區(qū)域Ω 分割為N個互不重疊的四邊形子區(qū)域且子區(qū)域滿足如下條件

為了數(shù)學(xué)書寫的方便以及計算效率方面考慮,首先通過如圖1 所示的變換將物理域(x,y)上任意Ωn轉(zhuǎn)換到計算域(ξ,η)上的一個標(biāo)準(zhǔn)單元ΩS(為了敘述簡潔,下面將省略下標(biāo)n).

圖1 將任意四邊形單元變換到一個標(biāo)準(zhǔn)單元Fig.1 Mapping of an arbitrary physical space quadrilateral element to the standard reference quadrilateral element

經(jīng)過以上變換后,子區(qū)域Ωn上的方程(4) 變換為

其中,頂標(biāo)“?”代表其為計算域(ξ,η)上的變量,計算域(ξ,η)和物理域(x,y)上的變量存在如下關(guān)系

其中,右上標(biāo)“δ”代表近似解,“D”表示是單元上的局部解,是求解點(ξi,ηj)處的近似解,li和lj是相對應(yīng)的拉格朗日插值多項式.在單元內(nèi)部是連續(xù)的,但在相鄰單元之間是間斷的.圖2 中給出了p=2 時在四邊形單元上求解點和后續(xù)用于通量修正的邊界面通量點的分布.

圖2 兩個四邊形單元上求解點(圓圈)和通量點(方塊)的分布(p=2)Fig.2 Arrangement of solution points(circles)and flux points(squares)on a pair of quadrilateral elements(p=2)

在計算散度項之前,首先進行通量項的構(gòu)造和計算.參照方程(9) 中近似解的構(gòu)造,同樣可以獲取局部通量.但由于僅在單元內(nèi)部連續(xù),而在單元之間并不連續(xù).為了滿足格式守恒性的要求,FR 通過增加修正項的方式構(gòu)造在單元界面處連續(xù)的通量.因此,連續(xù)的全局通量由單元內(nèi)的局部通量和修正通量構(gòu)成上式中,局部通量通過p階的插值多項式進行近似

綜上,經(jīng)過FR 格式空間離散之后的半離散方程為

上述常微分方程組可進一步通過顯式或隱式的時間推進方法進行求解.

3 隱式時間推進方法

二維歐拉方程經(jīng)過空間離散之后得到的半離散方程可進一步寫作如下形式

在上式中為M×N×4 維度的全局解向量,其中,M和N分別代表網(wǎng)格單元總數(shù)以及每個單元內(nèi)求解點的個數(shù),UM(t)為第M個單元上的局部解,是第M個單元第N個求解點上包含守恒變量的解向量.

代表右端殘差項,RM(U)是第M個單元上的局部殘差,為第M個單元第N點上殘差.

由于本文計算的是定常流動,無需考慮時間精度,因此時間項的離散使用一階向后隱式差分格式.方程(15)經(jīng)過該方法離散后可得

式中,Un=U(tn)和Un+1=U(tn+1)分別為tn和tn+1時刻的解.ΔTn為一個包括時間步長的對角矩陣,其與Un具有相同的維數(shù),其具體形式為

由于方程(16)中的R(Un+1)是非線性的,通過泰勒展開對其線性化,則有

上式中,?R/?U|Un為全局雅可比矩陣.將方程(18)代入方程(16)并整理得

上式中,tn+1時的解Un+1需要先求解線性方程組才能獲取.

3.1 塊雅可比迭代

線性方程組的求解可分為直接法和迭代法.當(dāng)網(wǎng)格規(guī)模較大時,使用直接法求解效率低下且內(nèi)存占用較大.因此本文通過迭代法求解方程(19),所用的迭代方法為雅可比迭代.

首先對方程(19) 左端的矩陣進行分解.令A(yù)=[(ΔTn)-1-?R/?U|Un].則A=D-N可看作是由分塊對角矩陣D和剩余非對角部分N組成.其中,分塊對角矩陣D為

與此同時,N的矩陣形式是

式中的N為一稀疏矩陣.經(jīng)過上述分解之后,對全局線性方程組(19)進行迭代求解,可得

在上式中,k為雅可比迭代指標(biāo).當(dāng)k→∞時,雅可比迭代解Uk+1趨近于下一個時刻的解Un+1.對于定常問題的計算,k只需要設(shè)置為有限的值,當(dāng)n→∞,Uk+1和Un最終都將趨向于所求的定常解.

為了避免非對角部分N的計算.對方程(23)兩側(cè)同時減去D(Uk-Un),則有

上式經(jīng)進一步整理后可得

和方程(18)類似,將R(Uk)在R(Un)處進行泰勒展開

上式可進一步寫成

將方程(27)代入到方程(25)中消去全局矩陣相關(guān)的一項?R/?U|Un(Uk-Un),并最終整理可得

式中,D=diag{D1,D2,···,DM}.

方程(28)和方程(19)相比,它的左端矩陣D是一個分塊對角矩陣,這使得求解全局線性方程組(19)的問題解耦為求解一系列局部單元上的小型線性方程組.

值得注意的是,全局矩陣A的分解方式并不唯一,如將其分解為分塊上三角矩陣U和分塊下三角矩陣L,即為Gauss-Seidel 迭代法.此時需要求解的方程的左端矩陣為分塊下三角矩陣L.由于此時Uk+1存在數(shù)據(jù)依賴問題,無法同時更新,為了并行地求解該方程,需要使用圖著色技術(shù),所得方法即為更復(fù)雜的MCGS 法[26].

3.2 局部線性方程組的求解

方程(28)可進一步寫成一系列局部單元線性方程組.對于第i(1 ≤i≤M)個單元,則有上式中,i為單元指標(biāo),n是時間迭代指標(biāo),k為雅可比迭代指標(biāo).為第i個單元對應(yīng)的單元雅可比矩陣,在下一小節(jié)中將討論其具體計算方法.以上獲取的M個局部線性方程組(29),其左端矩陣規(guī)模較小,在程序中通過調(diào)用cuBLAS 函數(shù)庫使用LU 分解法實現(xiàn)同步求解.

3.3 單元雅可比矩陣的計算

本文中所用是解析法與自動微分相結(jié)合的方法.由式(13)可知,第i個單元上,FR 格式離散后的殘差Ri是不連續(xù)通量的散度和修正通量的散度之和.則單元雅可比矩陣為

4 數(shù)值實驗

本節(jié)中以二維Bump 流動和NACA0012 翼型無黏繞流兩個算例來展示隱式方法的精度和效率.為了進行計算效率的對比研究,所有算例均同時使用顯式和隱式兩類方法進行.顯式方法為具有TVD 性質(zhì)的三級Runge-Kutta(TVD RK3)方法[21],同時使用局部時間步長和多重網(wǎng)格(pmulti-grid)方法加速其收斂.隱式方法為本文中所提出的塊雅可比迭代的方法.隱式方法的時間步長ΔTn由CFL 數(shù)決定,CFL數(shù)初始值為2,并隨迭代步數(shù)指數(shù)增長,最大值為104.

所有算例均使用高階并行通量重構(gòu)軟件HAFR(heterogeneous architecture flux reconstruction) 在型號為英偉達GeForce RTX 2080ti 的單個GPU 卡上計算.為了保證FR 方法在邊界的精度不至于下降,算例中用到的網(wǎng)格是由開源網(wǎng)格生成軟件Gmsh 生成的二階曲邊網(wǎng)格.計算的收斂性通過密度殘差的L2范數(shù)判斷,當(dāng)其下降到10-14以下或平穩(wěn)而不再下降時便認(rèn)為計算已穩(wěn)態(tài)收斂并終止.

4.1 無黏亞聲速Bump 算例

本算例是 international workshop on high-order CFD methods[39]給定的一個二維無黏驗證算例.該問題的物理域為[-1.5,1.5]×[0,0.8],底邊的幾何描述為y=0.062 5e-25x2.左側(cè)為亞聲速入口邊界條件(給定入口總壓和總溫),右側(cè)出口為出口邊界條件(給定靜壓),上下邊界均為滑移壁面邊界條件.來流的馬赫數(shù)Ma=0.5 且攻角為零度.盡管該流動的精確解未知,但由于流動為亞聲速無黏流動且流場光滑,因此可通過熵來衡量解的精度,流場中熵的計算公式如下

式中的積分在HAFR 程序中通過10 階高斯積分法進行計算.

計算中所使用的網(wǎng)格由粗到細(xì)(即h-refinement)為：R1(12×4),R2(24×8),R3(48×16)和R4(96×32),圖3 是R2 網(wǎng)格的示例.對于同一網(wǎng)格,同時還將使用p=2,3,4 階的插值多項式(即p-refinement)進行計算.

圖3 本算例中所用的24×8 的四邊形網(wǎng)格Fig.3 Quadrangular mesh(24×8)used in the case

首先從數(shù)值精度的角度分析計算結(jié)果.表1 為HAFR 程序使用塊雅可比迭代隱式方法在四組網(wǎng)格上,插值多項式為p=2,3,4 階時的熵以及數(shù)值精度.從表中可知,對任意一階插值多項式階數(shù)而言,網(wǎng)格加密的過程中熵也隨之降低.同時,給定某組網(wǎng)格,提高插值多項式的階數(shù)也同樣會導(dǎo)致熵的降低.正如圖4 所示,在R2 網(wǎng)格上使用p=3 和p=4 階插值多項式得到的下壁面馬赫數(shù)等值線相比p=2 的更為光滑.此外還觀察到,對于p=2,3,4 階插值多項式,在4 組網(wǎng)格上獲取的數(shù)值精度,都近似達到p+1階這一理論值,這也驗證了HAFR 的計算精度.

圖4 在24×8 的網(wǎng)格上使用p=2,3,4 階插值多項式計算所得的馬赫數(shù)云圖Fig.4 Mach number contours calculated on the 24×8 mesh with p=2,3,4 polynomial orders

表1 在四組網(wǎng)格上,使用p=2,3,4 階多項式時的熵和數(shù)值精度Table 1 Entropies and numerical orders on four different grids with p=2,3,4 polynomial orders

下面從計算效率的角度進行分析.圖5 是HAFR程序使用顯式和隱式兩種方法在4 組網(wǎng)格上、插值多項式為p=2,3,4 階時密度殘差隨迭代次數(shù)的收斂情況.

在圖5 中,左右兩列子圖分別為使用塊雅可比迭代隱式和TVD RK3 顯式時間方法的計算結(jié)果.從迭代次數(shù)來看,對同一階插值多項式,顯式方法和隱式方法收斂所需的迭代次數(shù)都隨著網(wǎng)格的加密而增大.對于同一網(wǎng)格,使用更高階的插值多項式都會導(dǎo)致迭代次數(shù)的增多.但同時發(fā)現(xiàn),隱式方法的迭代次數(shù)隨自由度(degree of freedom,DoF)增幅較小,并且對于同一網(wǎng)格或同一階插值多項式,隱式方法所需的迭代次數(shù)要遠(yuǎn)遠(yuǎn)小于顯式方法.例如,在R1 網(wǎng)格上,隱式方法收斂的的迭代次數(shù)均在80 步以內(nèi),而顯式方法最快也需要1000 步左右才能收斂.

圖5 在四組網(wǎng)格上,插值多項式為p=2,3,4 階時,使用隱式和顯式方法所得到的密度殘差||Res(ρ)||2 隨迭代次數(shù)的變化Fig.5 Density residual||Res(ρ)||2 versus the iteration numbers for the explicit and implicit time-marching schemes on four different grids with p=2,3,4 polynomial orders

圖5 在四組網(wǎng)格上,插值多項式為p=2,3,4 階時,使用隱式和顯式方法所得到的密度殘差||Res(ρ)||2 隨迭代次數(shù)的變化(續(xù))Fig.5 Density residual||Res(ρ)||2 versus the iteration numbers for the explicit and implicit time-marching schemes on four different grids with p=2,3,4 polynomial orders(continued)

在密度殘差收斂方面,使用隱式方法的計算的密度殘差均能收斂到10-14這一量級.但顯式方法的收斂性隨著自由度的增大變差,其值下降到某個量級后將無法下降,例如使用p=2 階插值多項式時,在R1 網(wǎng)格上的密度殘差能夠收斂10-14量級,但在R4 網(wǎng)格上僅能下降到10-13這一量級.

在GPU 計算時間方面,表2 為HAFR 程序采用塊雅可比迭代隱式方法在4 組網(wǎng)格上,使用p=2,3,4 階插值多項式計算收斂所用的迭代次數(shù)和GPU 墻鐘時間.從表2 可知,隱式方法計算收斂所用時間隨著自由度的增大而增大.

表2 在四組網(wǎng)格上,插值多項式為p=2,3,4 階時,使用隱式時間方法時的迭代次數(shù)和墻鐘時間Table 2 Iteration numbers and wall-clock times for the implicit time-marching scheme on four different grids with p=2,3,4 polynomial orders

圖6 是HAFR 程序分別使用塊雅可比隱式和TVD RK3 顯式時間方法,在不同網(wǎng)格上,使用p=2,3,4 階插值多項式時所用的墻鐘時間對比.從圖6中可看到,對每一組計算,隱式方法收斂所用的時間明顯小于顯式方法.在所有算例中,隱式方法相比顯式至少有一個量級的加速效果.

圖6 在四組網(wǎng)格上,插值多項式為p=2,3,4 階時,使用隱式和顯式方法所用的計算時間Fig.6 Wall-clock times for the explicit and implicit time-marching schemes on four different grids with p=2,3,4 polynomial orders

4.2 NACA0012 翼型無黏繞流

在上一個Bump 算例中所使用的網(wǎng)格較為均勻,為了進一步考察在非均勻網(wǎng)格條件下塊雅可比迭代隱式方法的計算效率,下面對NACA0012 翼型無黏繞流進行數(shù)值模擬.在此算例中,自由來流的馬赫數(shù)為0.5 且攻角為零度.NACA0012 翼型表面為滑移壁面條件,遠(yuǎn)場為特征邊界條件,遠(yuǎn)場與翼型的距離大約為100 個弦長.

計算中用到的網(wǎng)格為Gmsh 生成的O 型網(wǎng)格,網(wǎng)格分為3 組：R1(32×32),R2(44×64)和R3(128×128).圖7 為R1(32×32) 網(wǎng)格的示意圖.對于同一網(wǎng)格,還將分別使用p=2,3,4 階的插值多項式(即prefinement)進行計算.

圖7 本算例中所用的32×32 的O 型四邊形網(wǎng)格Fig.7 O-type quadrangular mesh(32×32)used in this case

下面對 NACA0012 數(shù)值模擬結(jié)果進行分析.表3 為HAFR 程序采用塊雅可比迭代隱式方法在三組網(wǎng)格上,使用p=2,3,4 階插值多項式計算收斂所用的迭代次數(shù)和GPU 墻鐘時間.和上個算例中類似,可觀察到無論是網(wǎng)格加密還是提高插值多項式階數(shù),這兩種增大自由度的方式都會導(dǎo)致迭代次數(shù)和GPU計算時間的增多.

表3 在三組網(wǎng)格上,插值多項式為p=2,3,4 階時,使用隱式時間方法時的迭代次數(shù)和墻鐘時間Table 3 Iteration numbers and wall-clock times for the implicit time-marching schemes on three different grids with p=2,3,4 polynomial orders

圖8 是HAFR 程序分別使用隱式和顯式時間方法,在不同網(wǎng)格上,使用插值多項式為p=2,3,4 時所用的墻鐘時間對比.從圖8 中可知,對于每一組計算,使用塊雅可比迭代隱式方法收斂所用的時間相比TVD RK3 顯式方法至少有一個量級的加速效果.

圖8 在三組網(wǎng)格上,插值多項式為p=2,3,4 階時,使用隱式和顯式方法所用的墻鐘時間Fig.8 Wall-clock times for the explicit and implicit time-marching schemes on three different grids with p=2,3,4 polynomial orders

圖9 為HAFR 程序使用顯式和隱式兩種方法在三組網(wǎng)格上、插值多項式為p=2,3,4 階時密度殘差隨迭代次數(shù)的收斂情況.在圖9 中,左右兩列子圖分別為使用塊雅可比迭代隱式和TVD RK3 顯式時間方法的結(jié)果.從迭代次數(shù)來看,對同一階插值多項式,顯式方法和隱式方法收斂所需的迭代次數(shù)都隨著網(wǎng)格的加密而增大.對于同一網(wǎng)格,使用更高階的插值多項式會導(dǎo)致迭代次數(shù)的增多.但同時,相比顯式方法,隱式方法的迭代次數(shù)隨自由度增幅較小,并且對于同一網(wǎng)格或同一階插值多項式,隱式方法所需的迭代次數(shù)要遠(yuǎn)遠(yuǎn)小于顯式方法.

圖9 在三組網(wǎng)格上,插值多項式為p=2,3,4 階時,使用隱式和顯式方法所得到的密度殘差||Res(ρ)||2 隨迭代次數(shù)的變化Fig.9 Density residual||Res(ρ)||2 versus the iteration numbers for the explicit and implicit time-marching schemes on three different grids with p=2,3,4 polynomial orders

在密度殘差收斂方面,在3 組網(wǎng)格計算中,使用塊雅可比迭代隱式方法計算的密度殘差均能收斂到10-14這一量級.但TVD RK3 顯式時間方法的收斂性隨著自由度的增大變差,其值下降到某個量級后將無法下降.因此從殘差收斂角度,基于塊雅可比迭代的隱式方法要優(yōu)于TVD RK3 顯式方法.

綜上,隨著自由度的增大,顯式和隱式方法收斂所需的迭代次數(shù)都會增多.但對于同一網(wǎng)格或同一階插值多項式,隱式方法所需的迭代次數(shù)和計算時間要明顯小于顯式方法,且隱式方法的殘差下降幅度要比顯式方法更大.

5 結(jié)論

為了提高FR 方法的求解效率,本文提出一種基于雅可比塊迭代的高階FR 格式求解定常二維歐拉方程的單GPU 隱式時間推進方法.通過雅可比塊迭代的方式,克服了影響求解并行性的相鄰單元依賴,使得只需存儲和計算對角塊矩陣.最終將求解全局線性方程組轉(zhuǎn)化為求解一系列局部單元線性方程組.這些局部的線性方程組在GPU 上又可利用LU 分解法并行求解.

文中通過二維無黏Bump 流動和NACA0012 無黏繞流兩個算例表明:隨著自由度的增大,塊雅可比迭代隱式方法收斂所需的迭代次數(shù)和計算時間都會增多.但對于同一網(wǎng)格或同一階插值多項式,該隱式方法所需的迭代次數(shù)要明顯小于TVD RK3 顯式方法,且此隱式方法的殘差下降幅度要比TVD RK3 顯式方法更大.在計算時間方面,對于文中的每一組計算,雅可比塊迭代隱式方法收斂所用的GPU 計算時間明顯小于TVD RK3 顯式方法,且計算效率至少有一個量級的提高.

附錄

二維歐拉方程的雅克比矩陣為