非結(jié)構(gòu)網(wǎng)格上可壓縮Euler方程的并行算法

馬欣榮，史瑞琪

（咸陽師范學(xué)院 數(shù)學(xué)與信息科學(xué)學(xué)院，陜西 咸陽 712000）

非結(jié)構(gòu)網(wǎng)格具有優(yōu)越的幾何靈活性和普適性，對(duì)復(fù)雜區(qū)域邊界和約束情形有很強(qiáng)的適應(yīng)性，克服了結(jié)構(gòu)網(wǎng)格（直交網(wǎng)格）不能適應(yīng)在任意形狀和任意區(qū)域進(jìn)行剖分的缺陷，更容易處理復(fù)雜幾何外形?；诜墙Y(jié)構(gòu)網(wǎng)格的算法不斷發(fā)展，但往往存在存儲(chǔ)量較大和計(jì)算量大等缺點(diǎn)，存儲(chǔ)格式中需要存儲(chǔ)節(jié)點(diǎn)坐標(biāo)以及標(biāo)識(shí)單元間的相鄰接的關(guān)系，直接影響了這類算法的計(jì)算效率。為了提高計(jì)算效率，各種高效的數(shù)值計(jì)算方法不斷發(fā)展，尤其是適合于線性和非線性系統(tǒng)方程的求解、不受控制方程時(shí)間和空間離散格式的限制，并且滿足最優(yōu)漸進(jìn)的性質(zhì)的多重網(wǎng)格方法[1-5]。

隨著計(jì)算機(jī)技術(shù)和數(shù)值分析方法的發(fā)展，工程中大型計(jì)算問題需要成千上萬億次級(jí)的運(yùn)算能力，串行算法已經(jīng)不能滿足計(jì)算需求，所以集群計(jì)算成了當(dāng)前科學(xué)計(jì)算的必要途徑之一。區(qū)域分解算法便是其中有效的求解方法之一，它具有諸多優(yōu)點(diǎn)：縮小計(jì)算規(guī)模；使用熟知的算法；無需采用整體一致網(wǎng)格，甚至子區(qū)域間可以使用不同的離散方法；子區(qū)域可以選用不同的數(shù)學(xué)模型；算法高度并行。由此在CFD應(yīng)用中得到了研究工作者的青睞。朱國林等[6]人主要研究了分布式并行處理系統(tǒng)的搭建與實(shí)施以及CFD計(jì)算程序的并行化，提出了并行化的關(guān)鍵是隱式算法中邊界上各類通量的信息傳遞的正確性。Hauser等[7]提出了求解Navier-Stokes方程直接數(shù)值模擬并行算法，降低了大型問題的海量計(jì)算成本。文獻(xiàn)[8]提出了三維Euler方程的三分區(qū)、四分區(qū)網(wǎng)格上的并行算法，負(fù)載平衡較好時(shí)能夠保持較高的并行效率。文獻(xiàn)[9]基于MPI中的染色分層通信技術(shù)，適當(dāng)處理塊邊界通信問題，采用隱式時(shí)間離散算法進(jìn)一步提高計(jì)算效率，求解了三維黏性流場。趙偉波等[10]在非結(jié)構(gòu)網(wǎng)格上設(shè)計(jì)了層次化的網(wǎng)格數(shù)據(jù)結(jié)構(gòu)實(shí)現(xiàn)了并行有限元方法，并應(yīng)用于彈性力學(xué)方程組的求解。

本文在文獻(xiàn)[11]基礎(chǔ)上，最大限度保證數(shù)據(jù)信息傳遞最少，各處理器保持負(fù)載均衡。采用LU-SGS隱式時(shí)間離散技術(shù)提高推進(jìn)效率，結(jié)合最大當(dāng)?shù)貢r(shí)間步長技術(shù)以及V循環(huán)-多重網(wǎng)格方法求解歐拉方程組，數(shù)值計(jì)算結(jié)果表明了本文算法的可行性，適合在集群系統(tǒng)下進(jìn)行像群計(jì)算。

1 數(shù)值離散算法

1.1流動(dòng)方程

在笛卡兒坐標(biāo)系下的非定常歐拉方程守恒形式：

其中Q、F分別為守恒變量和通量，定義如下：

其中ρ,e,p分別為密度、單位體積總內(nèi)能和壓強(qiáng)，ui為速度各方向的分量，δij為克羅內(nèi)克爾符號(hào)?？紤]理想氣體的熱力學(xué)特性：

其中γ為比熱比，空氣取為γ=1.4。

1.2流動(dòng)方程離散方法

采用

有限體積方法可得

在控制體i上對(duì)方程(2)進(jìn)行離散得到

式（3）為關(guān)于時(shí)間方向的一階常微分半離散方程組。相對(duì)于顯式時(shí)間離散格式而言，隱式時(shí)間離散格式可以放寬穩(wěn)定性條件，所以本文采用LU-SGS隱式時(shí)間離散格式。

運(yùn)用Euler向后差分離散以及線化處理后可以得到

其中L為嚴(yán)格下三角矩陣，D為純對(duì)角矩陣，U為嚴(yán)格上三角矩陣。

將式(4)變形并略去高階小量可以得到

在實(shí)際計(jì)算過程中，為了提高計(jì)算效率，本文采用了SOR內(nèi)迭代預(yù)處理技術(shù)。采用Roe[12]格式處理通量函數(shù)，限制器采用Venkatakrishnan[13]限制器。

2并行實(shí)現(xiàn)

2.1網(wǎng)格劃分

采用多級(jí)圖分區(qū)方法劃分網(wǎng)格，保證各子區(qū)域交界面數(shù)量最少且網(wǎng)格獨(dú)立保存，僅在相鄰子區(qū)域間數(shù)據(jù)傳遞，實(shí)現(xiàn)同步計(jì)算，較好地減少通信等待時(shí)間以提高計(jì)算效率。邊界分為物面邊界、遠(yuǎn)場邊界和分區(qū)內(nèi)邊界。物面滿足無穿透條件，遠(yuǎn)場滿足無反射條件。人工邊界引入虛擬單元進(jìn)行數(shù)據(jù)傳遞，見圖1。

圖1子分區(qū)邊界

2.2數(shù)據(jù)傳遞

3 數(shù)值實(shí)驗(yàn)

3.1繞NACA0012翼型跨音速無黏流動(dòng)

計(jì)算狀態(tài)：來流馬赫數(shù) Ma=0.8，攻角α=1.25°，收斂誤差為ε=10-10。圖2給出了八分區(qū)網(wǎng)格，網(wǎng)格單元共10 216個(gè)，網(wǎng)格節(jié)點(diǎn)共5 233個(gè)。包括虛擬單元的子區(qū)域單元網(wǎng)格數(shù)為：1 360、1 366、1 317、1 319、1 389、1 364、1 308、1 422。圖3給出了單個(gè)處理器、多處理器以及翼型實(shí)驗(yàn)值的壓力系數(shù)對(duì)比，可以看出數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)值吻合。實(shí)驗(yàn)數(shù)據(jù)見表1，可以看出本文算法10臺(tái)處理器并行效率保持在70%以上，優(yōu)于文獻(xiàn)[11]的計(jì)算結(jié)果。此外，隨著網(wǎng)格數(shù)量的增加，并行效率會(huì)隨之增加，基本呈線性變化。與參考文獻(xiàn)[14]相比較，單處理器與多處理器計(jì)算結(jié)果顯示物面壓力系數(shù)分布和激波位置與其保持一致。

圖2 NACA0012翼型八分區(qū)網(wǎng)格

表1 NACA0012翼型數(shù)值模擬并行性能

3.2 DLR-F6跨音速繞流

來流馬赫數(shù)：Ma=0.75，攻角：α=0°，圖4給出了四分區(qū)全局網(wǎng)格，網(wǎng)格單元共有1 453 849個(gè)，網(wǎng)格節(jié)點(diǎn)共275 561個(gè)。圖5給出了F6翼身組合體物體表面等密度線。單處理器計(jì)算99 100步需要152.056 9 min，10臺(tái)處理器計(jì)算需要27.080 1 min，并行效率為70.92%，優(yōu)于文獻(xiàn)[11]的計(jì)算結(jié)果，殘值下降12個(gè)量級(jí)。表2給出了具體計(jì)算數(shù)據(jù)，每迭代100步輸出一次計(jì)算結(jié)果。同樣通過數(shù)值實(shí)驗(yàn)?zāi)軌蝌?yàn)證隨著網(wǎng)格量的增加，并行效率會(huì)隨之增加，基本呈線性變化。

圖4 DLR-F6四分區(qū)網(wǎng)格

圖5 DLR-F6翼身組合體等密度線

表2 DLR-F6翼身組合體數(shù)值模擬并行性能

4結(jié)論

本文并行計(jì)算了雙曲守恒率，構(gòu)造虛擬單元，僅在相鄰處理機(jī)間進(jìn)行數(shù)據(jù)通信，使算法具有較好的并行性。結(jié)合當(dāng)?shù)貢r(shí)間步長、隱式LU-SGS時(shí)間離散方法以及多重網(wǎng)格技術(shù)加速收斂，在實(shí)際計(jì)算中，LU-SGS算法采用了預(yù)處理技術(shù)，相對(duì)于GMRES隱式格式計(jì)算效率相當(dāng)，而且能進(jìn)一步減小存儲(chǔ)要求，格式構(gòu)造簡單、易于編程實(shí)現(xiàn)，10臺(tái)處理器并行效率優(yōu)于文獻(xiàn)[11]的計(jì)算結(jié)果，保持在70%以上。為進(jìn)一步求解復(fù)雜外形黏性流場提供了依據(jù)。

參考文獻(xiàn)：

[1]SOLCHENBACH K,TROTTENBERG U.On the multigrid acceleration approach in computational fluid dynamics[C]//Proc of the 4thInt’l DFVLRSeminar on Foundations of Engineering Sciences：Parallel Computing in Science and Engineering.Bonn：Springer-Verlag，1988：145-158.

[2]YANG UM.Parallel algebraic multigrid methods-High performance preconditioners[C]//Proc of the Numerical Solution of Partial Differential Equations on Parallel Computers.Heidelberg，Berlin：Springer-Verlag，2006：209-236.

[3]HULSEMANN F，KOWARSCHIK M，MOHR M，et al.Parallel geometric multigrid[C]//Proc of the Numerical Solution of Partial Differential Equations on Parallel Computers.Heidelberg，Berlin：Springer-Verlag,2005：165-208.

[4]李宗哲.非結(jié)構(gòu)網(wǎng)格的并行多重網(wǎng)格算法研究[D].長沙：國防科技大學(xué)，2012.

[5]劉冰，陸中華，李新亮，等.基于GPU的多重網(wǎng)格Navier-Stokes解算器并行優(yōu)化算法研究[J].科研信息化技術(shù)與應(yīng)用，2013，4（3）：56-67.

[6]朱國林，徐慶新.計(jì)算流體力學(xué)并行計(jì)算技術(shù)研究綜述[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2002，20（S1）：1-6.

[7]HAUSER T，MATTOX T I，LEBEAU R P，et al.High-cost CFD on a low-cost cluster[C]//Proceedings of the 2000 ACM/IEEEConference on Supercomputing,Computer Society.New York：IEEE,2000：55-91.

[8]蘭黔章，呂曉斌.三維Euler方程的分區(qū)和并行計(jì)算[J].力學(xué)季刊，2001，22（4）：433-438.

[9]鄭冠男，鄧守春，韓同來，等.基于混合網(wǎng)格Navier-Stokes方程的并行隱式計(jì)算方法研究[J].應(yīng)用力學(xué)學(xué)報(bào)，2011，28（3）：211-219.

[10]趙偉波，劉青凱，楊揚(yáng).非結(jié)構(gòu)網(wǎng)格上彈性力學(xué)數(shù)值模擬的并行實(shí)現(xiàn)[J].計(jì)算機(jī)研究與發(fā)展，2015，52（5）：1153-1159.

[11]馬欣榮,劉三陽.非結(jié)構(gòu)網(wǎng)格上雙曲守恒率的并行算法研究[J].四川大學(xué)學(xué)報(bào)（工程科學(xué)版），2015，47（2）：123-128.

[12]ROE P L.Approximate Riemann solvers,parameter vectors，and difference schemes[J].Journal of Computational Physics，1981，43（2）：357-372.

[13]VENKATAKRISHNAN V.On the accuracy of limiters and convergence to steady state solutions[R].AIAA Paper 93-0880，1993.

[14]BARTH T J，JESPERSEN D C.The design and application of upwind schemes on unstructured meshes[C]//27thAerospace Sciences Meeting,January 9-12，Reno，Nevada：AIAA Paper 89-0366，1989：1-12.