適用于非定常流模擬的分布式并行GMRES方法

陳 龍,夏 健,田書玲

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

0 引 言

隨著計算流體力學(CFD)的高速發展,其在航空航天工程中的應用越來越廣泛。為了適應航空航天工程對于精細模擬的更高要求,目前CFD已經發展到千萬量級的計算網格階段。非常密的網格給數值計算方法的收斂性帶來巨大的挑戰,尤其是高雷諾數湍流問題和非定常問題。同時,高性能計算機硬件水平也不斷提高,使上萬CPU核心的并行計算成為可能。這就需要CFD高效且魯棒地處理航空航天工程中的復雜外形、適用于最新的高性能計算機及未來的發展[1-2]。綜上所述,提高CFD收斂性和并行計算效率的研究有著重要理論研究意義和工程應用價值。

全隱式方法是CFD提高收斂性和計算效率的重要途徑,全隱式方法的并行化是CFD研究難點之一。Luo等[3]最早將無矩陣LU-SGS方法作為GMRES的預處理算子,在共享內存系統中采用Open MP實現了算法的并行化。國內閻超等[4]在混合網格上采用Open MP實現了GMRES＋LU-SGS的并行化。上述的算法均有非常高的收斂效率,但由于共享內存系統的限制,并行效率并不高。Zingg等[5-6]采用近似Schur并行預處理算子實現了塊結構網格上GMRES的分布式并行化,采用無矩陣方法處理矩陣向量積,取得了非常好的收斂效果和并行效果,在一定條件下出現了加速比的超線性,NASA CRM模型千萬量級網格湍流模擬僅需3000步收斂。同一算例半隱式的LU-SGS方法需要數萬步的迭代才能收斂[7]。燕振國等[8]將GMRES應用于高階耗散緊致格式,張健等[9]基于PETsc科學計算工具包建立三維混合網格分布式并行GMRES算法,龔小權等[10]將GMRES應用于間斷Galerkin有限元,均取得了不錯的效果。綜合分析國內外研究的前沿,在混合網格上GMRES的分布式并行化是隱式方法研究的難點。需要恰當的處理矩陣運算并行化、預處理算子并行化、通訊模型、網格重排序和負載均衡,才能提高GMRES方法并行效率。近年來,計算機硬件高速發展,可用內存越來越大,大規模問題中存儲雅可比矩陣成為可能。全隱式方法中計算和存儲雅可比矩陣,并采用并行GMRES方法求解線性方程組是提高CFD收斂性的主要途徑。

本文采用變子空間數GMRES算法,結合Hybrid LU-SGS方法作為并行預處理算子,引入更魯棒的湍流模型,結合不同網格排序方法,在MPI并行環境中,發展出一套計算和存儲雅可比矩陣的分布式并行全隱式計算方法。在南京航空航天大學自主開發的OVERU軟件[11]中應用該方法,經驗證具有非常高的收斂性和并行效率,適用于定常和非定常流問題的數值模擬。

1 控制方程與離散方法

控制方程為任意拉格朗日歐拉(ALE)形式的三維雷諾平均N-S方程(RANS),湍流模型為一方程SA模型,本文采用的是一種新型的Negative-SA模型[12],該模型中ν~允許為負值,以抑制網格質量較差時的非物理解,可以改善整體隱式方法的魯棒性,提高迭代中CFL數的取值,從而間接達到提高全隱式方法收斂性的目的。Negative-SA模型方程為:當ν~≥0時

式中各符號與原始SA模型一致。其思想是當ν~≥0求解原始SA模型,當出現非物理解ν~＜0時求解式(2),式(2)源項中破壞項cw1(ν~/d)2與原SA模型符號相反,通過迭代求解使得ν~重新更新為正值。

OVERU軟件采用非結構混合網格有限體積方法求解控制方程。空間離散采用格點格式,可處理多種單元類型的混合網格。基于多種常用的通量差分格式HLLC、ROE、HLLEM和通量分裂格式AUSMPW＋等,通過U-MUSCL方法進行高階重構,限制器采用Venkatakrishna限制器。非定常計算采用雙時間步長方法進行二階精度的時間離散,內迭代采用本文所發展的全隱式方法。

2 并行全隱式方法

2.1 GMRES方法

全隱式方法是定常、非定常CFD計算中提高收斂性和計算效率有效途徑之一。其中的Krylov子空間迭代法中的GMRES方法[13]是一種高效的解線性方程組方法,結合不同預處理算子在CFD數值計算領域廣泛使用。本文的分布式并行GMRES方法的架構如圖1所示,含全隱式方法的并行化實現及優化。

在重啟型GMRES算法的基礎上,本文采用變子空間數GMRES算法,其初始化過程為v0＝RAΔW0,r0:＝P－1v0,β:＝‖r0‖,v1:＝r0/β。初始化后,計算流程如下:

圖1 分布式并行計算架構Fig.1 Framework of distributed parallel

重啟型GMRES中通過殘差收斂判據確定是否需要GMRES重啟,每個重啟過程的子空間數不變均為m。通過增加殘差極小值問題的求解次數(式(3)中增加第10步),提前判斷是否滿足殘差收斂判據。當滿足收斂判據時,即使子空間數小于m也退出GMRES內循環,從而減少子空間迭代數,提高計算效率。存儲空間方面,GMRES方法需要(m＋2)×nnodes×neqn×8字節內存空間。其中m為Krylov子空間數,nnodes為網格點數,neqn為方程數,如三維Euler方程為5。

變子空間數GMRES算法,可以在基本不影響收斂性的前提下,減小了每次迭代計算時間,從而提高整體計算效率。變子空間數GMRES算法相比于原重啟型GMRES算法在迭代啟動若干步的魯棒性上有一定程度的提高,改善了全隱式方法的啟動問題。

2.2 雅可比矩陣計算方法

隱式方法在處理雅可比矩陣時可分為無矩陣方法和計算矩陣方法。其中無矩陣方法一般采用數值差分方法方法直接得到矩陣向量積,不需要計算和存儲雅可比矩陣,節約內存空間。有矩陣方法中雅可比矩陣需要提前計算并存儲,雖然增加了內存需求和計算量,但是矩陣方法比無矩陣方法具有更好的收斂性。

全隱式方法中雅可比矩陣以雙精度實數存儲需要(nnodes＋2nedges)×neqn2×8字節內存空間,處理大規模問題時受到內存的限制往往采用無矩陣方法。近幾年計算機硬件的發展,內存空間不斷增大,價格也不斷降低。內存空間已經逐漸不再制約有矩陣方法的應用。雅可比矩陣的計算方法可通過人工編程或自動微分來實現。本文采用人工編程方法實現雅可比矩陣的計算,人工編程方法可以帶來較優化的代碼和較高的效率。雅可比矩陣由無黏項和黏性項組成,無黏項采用Van-Leer通量分裂雅可比矩陣,黏性項則基于薄層N-S假設。為提高計算效率,實際計算中雅可比矩陣每隔若干次迭代重新計算一次,在這幾次迭代中凍結。采用這個策略在不影響收斂性前提下,減少整體計算時間。

無矩陣方法中的矩陣向量積本文采用Luo等[3]的近似方法,如式(4)所示,簡化通量雅可比計算,可節約差分近似中兩次殘差的計算量:

2.3 GMRES并行化

基于MPI非阻塞通訊模型,采用分區并行策略實現GMRES的分布式并行化。雅可比矩陣的對角項按點存儲,非對角項按邊存儲,均存儲在局部內存中。通過分區邊界上的虛擬點進行數據通信,GMRES中矩陣向量積(式(3)中第2步)僅需相鄰分區之間的通信,減小通信時間。GMRES中向量內積(式(3)中第4、7步)在各分區計算完成后進行全局歸約操作,得到全局Hessenberg矩陣,通信量極小。采用Metis庫進行分區實現并行負載均衡,分區后可以保證各分區的網格點數基本一樣,同時分區交界面分割的邊數盡量少。可擴展性是并行計算方法的重要指標[14],非結構網格分區并行實現方法的可擴展性與分區方法的選取有直接關系,問題規模和分區數等比例增大后,如果分區交界面分割的邊數基本不變,則并行方法的可擴展性會很好。

Wissink等[15]在DP-LUR基礎上進行改進,發展出并行化的Hybrid LU-SGS方法。該方法的本質是在各處理器的網格分區內部使用LU-SGS方法,進行Gauss-Seidel迭代;在網格分區的邊界上,需要數據通信的網格點上使用DP-LUR方法,進行Jacobi迭代。Hybrid LU-SGS方法具有較高的并行效率,且計算時間比DP-LUR方法的計算時間減少約45%。相比其它預處理算子Hybrid LU-SGS效率更高,因此本文選擇Hybrid LU-SGS作為GMRES方法中的預處理算子P(式(3)中第2步),計算流程如下:

交換分區邊界數據,然后邊界點上有:

分區內部進行LU-SGS向前掃和向后掃:

2.4 網格重排序

本文所發展的數值方法中的網格重排序包含兩個方面:

(1)非結構網格LU-SGS算法中的網格重排序[3-4]。經過分層和染色兩個步驟將網格重排之后,保證每一層中的點與本層其它點不相連、與任一點相連的上下對角點的個數相當,目的是改善上三角矩陣U和下三角矩陣L的平衡性,構造非結構網格超平面,提高LU-SGS方法的效率。本文采用全局排序方法將分層信息發送到各計算分區。

(2)網格節點序號的重排序[16]。對網格節點序號進行重排序,其目的是減小稀疏矩陣的帶寬,使矩陣中非零項靠近主對角項,減小矩陣存儲空間,提高解線性方程組的效率。CFD計算中網格節點序號的重排序還可以減小計算過程中緩存未命中幾率,提高計算效率。本文采用RCM(Reverse-Cuthill-Mc Kee)方法進行網格節點序號的重排序。數值試驗表明經過RCM排序后可提高約10%的計算效率。

3 收斂特性驗證

3.1 二維定常流動

本節將對比不同時間離散方法在二維平面無黏流和NACA0012翼型黏性繞擾流算例中的收斂性。二維平面無黏流計算網格為三角形網格,如圖2,共1426個網格點。圖3給出了殘差收斂到機器精度時顯式方法(RK)、LU-SGS方法、無矩陣GMRES方法和矩陣GMRES方法(MGMRES)的對比。算例中最大子空間數m＝10,GMRES重啟次數為2,GMRES內殘差收斂判據ε＜0.1,本文后續算例均采用此設置。三種隱式方法均在1000步內收斂到機器精度,收斂性大大優于顯示方法,收斂性最好的是MGMRES方法。

圖2 平面流計算網格Fig.2 Computational grid for plane flow

圖3 平面流殘差收斂曲線Fig.3 Convergence history for plane flow

NACA0012翼型黏性擾流計算網格采用四邊形單元,共23萬網格點,如圖4。計算條件為:來流馬赫數Ma∞＝0.15,迎角α＝10°,雷諾數為Re＝6×106。該算例中數值方法的收斂性受到密網格和高雷數嚴苛的考驗。LU-SGS方法CFL數取1×106,GMRES和MGMRES方法均沒有啟動問題,CFL數可直接取300。圖5分別給出迭代收斂曲線和CPU時間收斂曲線。可見殘差收斂6個量級,LU-SGS方法需要約80 000步,GMRES方法約10 000步,MGMRES方法僅需2500步。計算時間上由于MGMRES每迭代步的計算時間較長,迭代初期其計算效率接近或低于GMRES方法。迭代后期MGMRES的收斂速度優勢明顯,僅需GMRES方法不到二分之一的計算時間。無矩陣GMRES方法在計算效率較LU-SGS方法提高約20%。需要特別說明的是收斂6個量級主要是為了對比不同方法的收斂性,實際計算中并不需要如此大量的迭代。表1給出了本文計算得到的升阻力系數與文獻[17]中不同程序計算結果的對比,符合得很好。

圖4 NACA0012計算網格Fig.4 Computational grid for NACA0012 case

表1 NACA0012結果與其他代碼比較Table 1 Present NACA0012 results compared with those of other codes

圖5 NACA0012擾流收斂曲線Fig.5 Convergence history for NACA0012 case

3.2 翼身組合體擾流

對DLR-F6翼身組合體[18-19]有黏和無黏繞流進行計算分析,其中無黏計算網格采用四面體單元約90萬網格點503萬單元,黏性網格采用六面體單元,約580萬網格點,如圖6。計算條件為:來流馬赫數M∞＝0.75,迎角α＝0.52°。黏性計算中雷諾數Re＝3×106。LU-SGS方法CFL數取1×106,GMRES方法CFL數取300,MGMRES方法無黏算例CFL數取300。MGMRES方法在黏性算例遇到了啟動問題,本文種采用變CFL數取值的策略解決啟動問題,如下式所示:

首先進行網格節點序號重排序方法的驗證,圖7給出了采用網格節點序號的重排序方法(RCM)和不采用網格節點序號的重排序方法(NO-RCM)的對比,其中TET表示無黏擾流算例,HEX為黏性擾流算例。可見RCM方法減小了計算過程中緩存未命中幾率,經過RCM排序后可提高約10%的計算效率。

圖6 DLR-F6計算網格Fig.6 Computational grid for DLR-F6 case

圖7 RCM排序對計算時間的影響Fig.7 Effect of RCM reorder on CPU time

圖8 重啟型和改進GMRES方法收斂曲線對比Fig.8 Comparison of convergence history between original and modified GMRES

圖8給出了黏性擾流算例中重啟型GMRES方法和本文改進的變子空間數GMRES方法殘差收斂曲線的對比,兩者收斂效率一致。為達到GMRES內殘差收斂判據,當重啟型GMRES方法中一般需要1次重啟,共10＋10次子空間運算,而變子空間數GMRES方法在重啟后子空間數會減少,共10＋n次子空間運算(n≤10)。變子空間數GMRES方法較重啟型GMRES方法,計算時間減少約15%。

圖9給出了三種隱式方法在進行黏性擾流計算時的迭代收斂曲線和CPU時間收斂曲線。其中MGMRES方法的收斂性最好,僅需800步左右殘差下降3個量級。殘差下降3個量級MGMRES方法較GMRES方法計算時間減少約37%,整體計算時間是LU-SGS方法的四分之一不到。并且MGMRES方法在殘差下降3個量級后收斂性依然較好,而GMRES方法收斂曲線出現一定程度波動。圖10給出了升力系數的迭代和CPU時間收斂曲線。從升力系數上判斷MGMRES方法收斂效率也是最高,略高于GMRES方法,計算效率約為LU-SGS方法的4倍多。

圖9 DLR-F6擾流收斂曲線Fig.9 Convergence history for DLR-F6 case

針對第六屆AIAA阻力測試會議中的NASA CRM翼身組合體構型,采用NASA提供的混合網格約3000萬網格點1.2億網格單元,驗證本文的方法在數千萬量級網格上的收斂性。計算條件為:來流馬赫數Ma∞＝0.85,迎角α＝2.75°,雷諾數Re＝5×106。圖11給出了三種隱式方法在進行黏性擾流計算時的升力系數收斂曲線。其中MGMRES方法升力系數收斂僅需1000步左右,而GMRES方法需要約1600步。相同計算狀態及網格條件下,文獻[20]中NASA的FUN3D軟件采用全隱式方法需要約3000步收斂。本文的方法在千萬量級網格上的也具有很好的收斂性。

圖10 DLR-F6擾流升力系數收斂曲線Fig.10 Lift coefficient convergence history for DLR-F6 case

圖11 NASA CRM擾流升力系數收斂曲線Fig.11 Lift coefficient convergence history for NASA CRM case

3.3 戰斗機非定常氣動彈性模擬

本節采用F16戰斗機外形的非定常氣動彈性數值模擬算例,考察所發展的隱式方法在三維非定常黏性擾流計算中的效果。計算網格采用混和網格,圖12為機身物面網格和對稱面網格,約214萬網格點。計算條件為:來流馬赫數Ma∞＝0.8,迎角α＝4.0°,CFL數取值與3.2節一致。結構有限元模型基于ECERTA項目[21]提供的結構模態數據。在時域內耦合求解非定常N-S方程和結構模態方程,實現飛機氣動彈性的數值模擬。其中流固界面插值方法為薄平板樣條法(TPS),動網格方法采用高效的Delaunay圖映射方法。雙時間步長方法中內迭代殘差下降3個量級為判斷其收斂的依據。

圖12 F16計算網格Fig.12 Computational grid for F16 case

圖13 F16非定常殘差收斂曲線Fig.13 Convergence history for F16 aeroelastics flow

圖13給出了三組雙時間步長內迭代殘差曲線,MGMRES方法僅需約7步迭代殘差下降3個量級,GMRES方法需要約8步,LU-SGS方法則需要約20步。收斂性仍然是MGMRES方法最好。計算時間上GMRES方法略快于MGMRES方法,較LU-SGS方法減少約23%。需要注意的是,此時間僅統計了流場計算的時間,未考慮氣動彈性計算中網格變形、流固界面插值等時間。

圖14給出了采用MGMRES方法計算得到的渦量等值面及壓強云圖,符合擾流基本規律。圖15給出了前四階模態廣義位移隨時間的變化曲線,各階廣義位移振幅逐漸減小,在此狀態飛機未發生顫振。

圖14 F16壓強分布和渦量等值面Fig.14 F16 surface pressure and vorticity iso-surface

圖15 廣義位移時間歷程Fig.15 Generalized displacement history

3.4 旋翼前飛流場

將本文方法發展到非結構重疊網格上,通過旋翼非定常流模擬研究方法的收斂特性。采用2.5m級單旋翼模型,包括4片槳葉。旋翼實度0.077,翼型采用NACA0012,弦長0.072 m,幾何扭轉0°,前飛馬赫數0.127。計算網格方面,如圖16所示采用的非結構重疊網格系統,共包含6個子網格,其中兩塊為背景網格點數分別為601萬和139萬,每片旋翼葉片的網格隨旋翼運動,網格點數為107萬,總網格數約1200萬。

雙時間步長方法中內迭代殘差下降4個量級為判斷其收斂的依據。在內迭代中MGMRES方法遇到了魯棒性問題,非定常CFL最大僅可取到80,而LU-SGS方法可取1×106。圖17給出了采用MGMRES方法計算得到的Q等值面,其中槳渦干擾現象明顯。圖18給出了兩組內迭代殘差曲線,MGMRES方法僅需約35步迭代殘差下降4個量級,LU-SGS方法則需要約107步。收斂性仍然是MGMRES方法最優。計算時間上MGMRES方法LU-SGS方法減少約40%。但是,內迭代殘差下降3個量級可滿足一般的非定常計算需要,此時LU-SGS方法在計算時間上還要少于MGMRES方法。主要原因是此算例中MGMRES方法的CFL數較小,當達到GMRES內殘差收斂判據ε＜0.1時,子空間數僅有1~3個,降低了GMRES方法的收斂速度。因此提高GMRES方法的魯棒性的研究在未來仍然是很有必要的。

圖16 旋翼重疊網格系統Fig.16 Overset grids system for rotors

圖17 Q等值面Fig.17 Q criterion iso-surface

圖18 旋翼非定常殘差收斂曲線Fig.18 Convergence history for unsteady rotor flow

4 并行特性驗證

本節采用3.2節的翼身組合體繞流算例進行本文發展方法的并行特性驗證。圖19為各隱式方法內存空間需求,無矩陣GMRES方法比LU-SGS方法內存需求增加有限。MGMRES方法需要的內存最多,黏性算例在580萬網格上需要約13GB內存。MGMRES方法平均每千萬網格需要約22GB內存空間。在目前硬件條件下,采用MGMRES方法單臺工作站就能夠滿足千萬量級網格的內存需求,高性能分布式集群則可進行更大規模網格的計算。

圖19 內存需求Fig.19 Memory requirements

在國家超級計算天津中心的天河1A高性能計算機上,對所發展方法的并行特性展開研究。圖20為采用MGMRES方法、1和240個計算核心(NP)計算得到的機翼上四個站位(η)壓強系數分布與實驗值的對比,兩者符合得較好,驗證了并行方法的正確性。

圖21給出了Hybrid LU-SGS方法和LU-SGS方法殘差收斂曲線。Hybrid LU-SGS方法在分區邊界網格點上使用DP-LUR方法而LU-SGS方法對分區邊界不做處理(圖中w/o表示)。可見在分區較少時兩種方法殘差收斂曲線非常接近。當240個分區時,分區邊界網格點不做處理的并行方法在4000步左右計算發散。因此,在進行較多分區并行計算時,采用Hybrid LU-SGS方法是必要的。圖22給出了Hybrid LU-SGS預處理算子和MGMRES方法并行加速比。兩種方法的加速比非常接近理想值,加速效率非常高,MGMRES方法的加速比略低于Hybrid LU-SGS方法。特別的,在240個計算核心的并行計算中均出現了加速比超線性(并行效率超過100%)。一般認為,在負載均衡和極端通信時間條件下,加速比超線性主要是因為計算機CPU緩存在更多分區時利用率更高。

圖20 壓強系數分布計算值與實驗值對比Fig.20 Comparison between experimental and computed pressure coef f icient distributions

圖21 Hybrid LU-SGS方法收斂曲線Fig.21 Convergence history using hybrid LU-SGS

圖22 并行加速比Fig.22 Parallel speedup

5 結 論

本文基于變子空間數GMRES算法,結合Hybrid LU-SGS方法作為并行預處理算子,引入更魯棒的湍流模型,結合不同網格排序方法,發展出一套計算和存儲雅可比矩陣的分布式并行全隱式計算方法。利用南京航空航天大學自主開發的OVERU軟件,在天河1A高性能計算機上對若干二維和三維擾流問題進行了數值模擬,經驗證具有非常高的收斂性和并行效率,結論如下:

(1)變子空間GMRES方法與重啟GMRES方法具有相同的收斂特性,計算時間可較原方法減少約15%。計算和存儲雅可比矩陣的全隱式方法收斂性優于無矩陣方法,較LU-SGS方法提高20%~200%計算效率。采用RCM網格重排序方法優化緩存命中可將計算效率再提高10%。

(2)矩陣GMRES方法每千萬網格約需要22GB的內存空間,當前硬件水平下可以滿足大規模網格計算的內存需求。

(3)分布式并行GMRES方法結合Hybrid LUSGS預處理算子,所發展的全隱式方法并行效率非常高。特別的,在240個計算核心出現了加速比的超線性。本文研究僅限于CPU并行計算,后續研究應面向未來E級超算的CPU＋GPU(或加速卡)的異構架構,開展異構架構下全隱式并行算法研究。