陳 龍,夏 健,田書玲

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

0 引 言

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

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

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

1 控制方程與離散方法

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

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

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

2 并行全隱式方法

2.1 GMRES方法

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

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

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

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

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

2.2 雅可比矩陣計算方法

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

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

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

2.3 GMRES并行化

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

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

交換分區(qū)邊界數(shù)據(jù),然后邊界點上有:

分區(qū)內(nèi)部進行LU-SGS向前掃和向后掃:

2.4 網(wǎng)格重排序

本文所發(fā)展的數(shù)值方法中的網(wǎng)格重排序包含兩個方面:

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

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

3 收斂特性驗證

3.1 二維定常流動

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

圖2 平面流計算網(wǎng)格Fig.2 Computational grid for plane flow

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

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

圖4 NACA0012計算網(wǎng)格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]有黏和無黏繞流進行計算分析,其中無黏計算網(wǎng)格采用四面體單元約90萬網(wǎng)格點503萬單元,黏性網(wǎng)格采用六面體單元,約580萬網(wǎng)格點,如圖6。計算條件為:來流馬赫數(shù)M∞＝0.75,迎角α＝0.52°。黏性計算中雷諾數(shù)Re＝3×106。LU-SGS方法CFL數(shù)取1×106,GMRES方法CFL數(shù)取300,MGMRES方法無黏算例CFL數(shù)取300。MGMRES方法在黏性算例遇到了啟動問題,本文種采用變CFL數(shù)取值的策略解決啟動問題,如下式所示:

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

圖6 DLR-F6計算網(wǎng)格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方法和本文改進的變子空間數(shù)GMRES方法殘差收斂曲線的對比,兩者收斂效率一致。為達到GMRES內(nèi)殘差收斂判據(jù),當重啟型GMRES方法中一般需要1次重啟,共10＋10次子空間運算,而變子空間數(shù)GMRES方法在重啟后子空間數(shù)會減少,共10＋n次子空間運算(n≤10)。變子空間數(shù)GMRES方法較重啟型GMRES方法,計算時間減少約15%。

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

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

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

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

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

3.3 戰(zhàn)斗機非定常氣動彈性模擬

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

圖12 F16計算網(wǎng)格Fig.12 Computational grid for F16 case

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

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

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

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

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

3.4 旋翼前飛流場

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

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

圖16 旋翼重疊網(wǎng)格系統(tǒng)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 并行特性驗證

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

圖19 內(nèi)存需求Fig.19 Memory requirements

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

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

圖20 壓強系數(shù)分布計算值與實驗值對比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 結 論

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

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

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

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