, , ,*,

(1. 西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072; 2. 中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽 621000)

0 引 言

大規(guī)模分布式并行計(jì)算快速發(fā)展，使計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)具備了解決復(fù)雜外形工程問題的能力。大規(guī)模CFD并行計(jì)算成為飛行器氣動(dòng)性能計(jì)算[1]、多體分離安全性評(píng)估[2]和飛行器外形優(yōu)化[3-5]的有效手段。

CFD并行計(jì)算屬于數(shù)據(jù)通信密集型，采用的數(shù)值格式和軟件并行設(shè)計(jì)對(duì)并行效率影響很大[6]。Knight首先系統(tǒng)性地闡述了并行計(jì)算機(jī)系統(tǒng)架構(gòu)和CFD并行程序設(shè)計(jì)模式[7]。 國內(nèi)陸林生和董超群等[8]研究了并行程序概念設(shè)計(jì)方法，陳剛和王磊等[9]提出了面向大規(guī)模并行的軟件實(shí)現(xiàn)方法。近年來，CFD并行計(jì)算廣泛應(yīng)用于航空航天等領(lǐng)域[10-12]，但并行規(guī)模都在萬核以下。

對(duì)于分布式并行系統(tǒng)，消息傳遞接口(Message Passing Interface, MPI)是最成熟和有效的并行程序?qū)崿F(xiàn)規(guī)范，廣泛應(yīng)用于計(jì)算力學(xué)領(lǐng)域。OpenMPI和MPICH是兩種廣泛使用的MPI編程程序包，都成功用于CFD并行程序編碼[13-14]。

在使用并行解算器模擬流動(dòng)前，通常需要采用合適的方法對(duì)生成的網(wǎng)格進(jìn)行多核并行分區(qū)。分區(qū)算法極大地影響著并行計(jì)算時(shí)的負(fù)載平衡和并行效率。對(duì)于非結(jié)構(gòu)網(wǎng)格，采用多級(jí)不規(guī)則圖算法的Metis程序庫[15]，是目前應(yīng)用最廣、最快速和高質(zhì)量的并行分區(qū)工具。由于沒有拓?fù)浣Y(jié)構(gòu)的限制，分區(qū)后的非結(jié)構(gòu)網(wǎng)格通常比結(jié)構(gòu)網(wǎng)格有更高的負(fù)載平衡性能。

非結(jié)構(gòu)網(wǎng)格的有限體積方法通常采用積分形式的雷諾平均Navier-Stokes方程和湍流模型方程[16]。二階MUSCL (Monotonic Upwind Scheme for Conservation Laws) 類數(shù)值格式具有類似緊致性質(zhì)，即某個(gè)網(wǎng)格單元通量的計(jì)算只需要臨近共面單元的流場(chǎng)信息，因而廣泛應(yīng)用于CFD計(jì)算軟件且發(fā)展較為成熟。并行計(jì)算時(shí)，僅需要通過MPI交換并行交界面兩側(cè)單元的流場(chǎng)數(shù)據(jù)，有利于并行計(jì)算的實(shí)現(xiàn)和減少并行傳遞的數(shù)據(jù)量。

本文通過基于非結(jié)構(gòu)混合網(wǎng)格和有限體積方法，采用緊致型二階數(shù)值格式，發(fā)展了適用于復(fù)雜工程問題的萬核級(jí)并行實(shí)現(xiàn)方法，測(cè)試和評(píng)估了并行效率。

1 數(shù)值方法

雷諾平均Navier-Stokes方程(RANS)守恒形式的積分方程為：

(1)

對(duì)流通量可以表示為：

(2)

式中，p是靜壓，ρ是密度，V是逆變速度:

V=nxu+nyv+nzw

(3)

H為焓值，定義為：

(4)

式中，e為流動(dòng)內(nèi)能。

黏性通量相對(duì)比較復(fù)雜，本文略去，具體內(nèi)容可以參考文獻(xiàn)[17]。

在網(wǎng)格單元內(nèi)積分，得到離散形式的控制方程：

(5)

式中，V是網(wǎng)格單元體積，N是與該網(wǎng)格單元共面網(wǎng)格單元總數(shù)。

控制方程的離散采用格心格式的有限體積法，時(shí)間項(xiàng)采用LU-SGS迭代[18]，湍流模型為一方程的SA模型[19]。對(duì)流項(xiàng)采用Roe格式[20]，黏性項(xiàng)采用中心格式[21]。

格心格式的有限體積方法，流場(chǎng)未知變量存儲(chǔ)在網(wǎng)格單元的形心上，如圖1所示的三角形形心點(diǎn)J1、J2和J3。然而，對(duì)流通量和黏性通量計(jì)算卻需要網(wǎng)格單元對(duì)應(yīng)的網(wǎng)格面心上的流場(chǎng)值，如K1、K2和K3。只有采用一階格式時(shí)，面上的流場(chǎng)值等于形心值。當(dāng)采用緊致型的二階格式時(shí)，網(wǎng)格面心上左右兩側(cè)的流場(chǎng)值由同側(cè)的網(wǎng)格單元體心值構(gòu)造。計(jì)算黏性通量時(shí)，面心流場(chǎng)值采用中心格式，以K1面為例，表示為:

(6)

圖1 通量計(jì)算模板示意圖Fig.1 Stencil for flux calculation

K2和K3面心流場(chǎng)值可以類似得到。黏性通量可以表示為:

Hv,K1=f(QK1)

(7)

黏性通量計(jì)算的具體表達(dá)式參考文獻(xiàn)[21]。

采用Roe格式時(shí)，對(duì)流通量可以表示為

Hc,K1=f(QK1,L+QK1,R)

(8)

式中，QK1,L和QK1,R分別為面K1左右兩側(cè)的流場(chǎng)值，通過各自同側(cè)的體心值構(gòu)造

(9)

采用Roe格式計(jì)算對(duì)流通量的具體表達(dá)式參考文獻(xiàn)[20]。

需要指出的是，如果采用更高階的數(shù)值格式，面心流場(chǎng)值的構(gòu)造將會(huì)用到更多臨近網(wǎng)格單元上的體心值。對(duì)于本文所述的二階格式，只需要緊鄰的網(wǎng)格單元體心值，因而并行計(jì)算時(shí)，只有并行交界面兩側(cè)的網(wǎng)格單元需要進(jìn)行流場(chǎng)數(shù)據(jù)的交換。

式(5)中的第二項(xiàng)通常稱為殘差項(xiàng)，記為:

(10)

則式(5)表示為:

(11)

為了增加計(jì)算的數(shù)值穩(wěn)定性和提高收斂速度，時(shí)間項(xiàng)通常采用隱式格式[17]，即

(12)

式中：

(13)

通過對(duì)殘差項(xiàng)的線化，可以得到

(14)

上式求解采用LU-SGS迭代方法，具體內(nèi)容可參考文獻(xiàn)[18]。定常流場(chǎng)由邊界條件決定，迭代過程中，邊界信息由邊界單元逐漸向流場(chǎng)內(nèi)部傳遞，直到收斂。通常情況下，可能需要迭代上千步甚至更多流場(chǎng)才能收斂，當(dāng)網(wǎng)格量較大時(shí)，計(jì)算量巨大，單核計(jì)算需要很長(zhǎng)的時(shí)間。對(duì)于工程復(fù)雜外形，并行計(jì)算是減少計(jì)算時(shí)間的有效方法。

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

采用LU-SGS迭代方法，當(dāng)前步流場(chǎng)計(jì)算依賴于前一步的流場(chǎng)，當(dāng)流場(chǎng)為非定常時(shí)依賴更多。但對(duì)于某一迭代步，采用緊致二階格式，網(wǎng)格單元的流場(chǎng)值只與相鄰一層網(wǎng)格有關(guān)。因此對(duì)于分布式并行計(jì)算系統(tǒng)，采用對(duì)空間進(jìn)行分區(qū)更適合CFD并行計(jì)算，并采用MPI進(jìn)行多個(gè)分區(qū)間流場(chǎng)數(shù)據(jù)的交換。并行實(shí)現(xiàn)算法主要包括網(wǎng)格分區(qū)算法和流場(chǎng)數(shù)據(jù)交換算法。

2.1 網(wǎng)格分區(qū)技術(shù)

非結(jié)構(gòu)網(wǎng)格分區(qū)方法采用基于圖論算法的Metis程序庫，該程序庫使用多級(jí)二分法或多級(jí)多支分區(qū)算法，如圖2所示。

圖2 Metis多級(jí)分區(qū)算法示意圖[15]Fig.2 Multi-level partition used by Meits[15]

分區(qū)算法以網(wǎng)格數(shù)據(jù)構(gòu)成的圖作為輸入，由網(wǎng)格面和網(wǎng)格單元組成。包含n個(gè)節(jié)點(diǎn)V和m個(gè)邊E的圖G定義為:

(15)

需要指出的是，對(duì)于格心格式的有限體積法，圖定義中的節(jié)點(diǎn)和邊分別指網(wǎng)格中的網(wǎng)格單元和網(wǎng)格面。輸入的圖信息通常采用壓縮存儲(chǔ)格式(Compressed Storage Format， CSR)，該格式廣泛應(yīng)用于稀疏圖的存儲(chǔ)。網(wǎng)格連接關(guān)系表示成圖的信息可以用兩個(gè)整數(shù)數(shù)組來指定：xadj和adjncy，數(shù)組長(zhǎng)度分別為n+1和2m。xadj[i+1]-xadj[i]的值表示與第i個(gè)網(wǎng)格單元共面的網(wǎng)格單元個(gè)數(shù)，數(shù)組adjncy中第xadj[i] 到第xadj[i+1] 位置存儲(chǔ)了相應(yīng)的共面網(wǎng)格單元的編號(hào)。以圖1網(wǎng)格數(shù)據(jù)為例，數(shù)組xadj和adjncy存儲(chǔ)數(shù)據(jù)如圖3所示。

圖3 非結(jié)構(gòu)網(wǎng)格連接關(guān)系CSR表示Fig.3 Graph data in CSR format of unstructured grid

圖4給出了M6機(jī)翼網(wǎng)格并行分區(qū)表面網(wǎng)格的示意圖，共3個(gè)并行分區(qū)，其中相同顏色的網(wǎng)格單元屬于相同的計(jì)算進(jìn)程。

2.2 數(shù)據(jù)交換并行實(shí)現(xiàn)

本文并行算法基于分布式大規(guī)模并行計(jì)算機(jī)系統(tǒng)，采用所有CPU都參與迭代計(jì)算的對(duì)等并行模式，MPI數(shù)據(jù)傳遞采用非阻塞通信，盡量加大計(jì)算和通信的重疊，只在必要的地方進(jìn)行數(shù)據(jù)同步以保證計(jì)算的正確性。

圖4 M6機(jī)翼網(wǎng)格分區(qū)圖Fig.4 Partitioned grid of M6 wing with Metis

采用對(duì)等并行模式，首先需要對(duì)網(wǎng)格進(jìn)行前處理，按指定的并行規(guī)模采用Metis劃分網(wǎng)格。為了減少每個(gè)計(jì)算進(jìn)程內(nèi)存消耗以提高程序的運(yùn)行效率，每個(gè)分區(qū)網(wǎng)格只保存必需的局部網(wǎng)格數(shù)據(jù)和并行分區(qū)邊界上網(wǎng)格單元及頂點(diǎn)的對(duì)應(yīng)關(guān)系。網(wǎng)格分區(qū)完成后，每個(gè)計(jì)算進(jìn)程讀入相應(yīng)分區(qū)的網(wǎng)格，獨(dú)立地逐一計(jì)算各自分區(qū)中每個(gè)網(wǎng)格單元的通量。

格心格式的有限體積法，并行邊界間數(shù)據(jù)傳遞需要的對(duì)應(yīng)關(guān)系包括網(wǎng)格單元與網(wǎng)格單元和網(wǎng)格頂點(diǎn)與網(wǎng)格頂點(diǎn)間的對(duì)應(yīng)關(guān)系，其中網(wǎng)格單元間對(duì)應(yīng)關(guān)系主要用于流場(chǎng)變量的迭代，網(wǎng)格頂點(diǎn)間的對(duì)應(yīng)關(guān)系用于二階格式中梯度的計(jì)算。并行網(wǎng)格對(duì)應(yīng)關(guān)系建立耗時(shí)長(zhǎng)，然而對(duì)于定常問題只需要計(jì)算一次，程序設(shè)計(jì)時(shí)宜采用存儲(chǔ)該對(duì)應(yīng)關(guān)系，以空間換時(shí)間，提高并行效率。

對(duì)于串行計(jì)算，封閉的邊界條件決定了當(dāng)前流場(chǎng)的特征，邊界信息由邊界網(wǎng)格單元向內(nèi)部網(wǎng)格單元傳播。并行計(jì)算時(shí)，對(duì)于單個(gè)計(jì)算進(jìn)程，只有部分邊界網(wǎng)格存在相應(yīng)的分區(qū)網(wǎng)格中，邊界條件不再封閉。將并行交界面上的網(wǎng)格面視為特殊類型的邊界面時(shí)，則邊界面滿足封閉性質(zhì)。同時(shí)，與真實(shí)邊界面類似地采用虛擬邊界單元，真實(shí)邊界面對(duì)應(yīng)的虛擬單元用于施加邊界條件，并行邊界面對(duì)應(yīng)的虛擬單元?jiǎng)t用來存儲(chǔ)并行交換的數(shù)據(jù)。在計(jì)算網(wǎng)格單元的通量時(shí)，并行邊界面與真實(shí)邊界面采用相同的計(jì)算方法，程序設(shè)計(jì)亦能與串行程序保持一致。并行邊界與虛擬單元設(shè)計(jì)如圖5所示，則對(duì)應(yīng)任意并行計(jì)算進(jìn)程，計(jì)算及存儲(chǔ)模式和串行時(shí)保持一致。

MPI并行接口中并行數(shù)據(jù)必須以連續(xù)內(nèi)存的方式提供，程序設(shè)計(jì)必須考慮減少并行計(jì)算時(shí)額外數(shù)據(jù)拷貝的開銷以提高并行效率。格心格式有限體積法任何分區(qū)內(nèi)的某個(gè)單元只對(duì)應(yīng)唯一其它分區(qū)內(nèi)的唯一單元，可以采用先物理網(wǎng)格單元后虛擬網(wǎng)格單元的網(wǎng)格排序策略，使對(duì)應(yīng)的每一個(gè)并行分區(qū)需要接收的數(shù)據(jù)在內(nèi)存中集中連續(xù)存放，如圖6所示，消除了網(wǎng)格單元并行數(shù)據(jù)的拷貝過程。

圖5 并行邊界及虛擬單元Fig.5 Ghost elements for data exchange between different processors

圖6 并行物理網(wǎng)格單元及虛擬網(wǎng)格單元排序存儲(chǔ)示意圖Fig.6 Memory pattern for physical elements and ghost elements of grid data

3 并行測(cè)試結(jié)果

采用兩個(gè)算例分別對(duì)本文非結(jié)構(gòu)流動(dòng)解算器的并行正確性和并行效率進(jìn)行測(cè)試。算例一使用網(wǎng)格量相對(duì)較少的旋成體外形，迭代計(jì)算直到流場(chǎng)收斂，比較分析串行計(jì)算和并行計(jì)算結(jié)果的一致性；算例二使用運(yùn)輸機(jī)構(gòu)型，其網(wǎng)格量大，用于測(cè)試大規(guī)模并行計(jì)算的效率，其最大的并行規(guī)模達(dá)18816核。

并行加速比S和效率E定義為：

(16)

式中：n為測(cè)試并行核數(shù)，T為每步計(jì)算平均時(shí)間；n0為并行測(cè)試基數(shù)核數(shù)，T0為基數(shù)核數(shù)對(duì)應(yīng)的每步計(jì)算平均時(shí)間。

并行計(jì)算使用中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所分布式并行計(jì)算機(jī)系統(tǒng)。該系統(tǒng)包含300個(gè)計(jì)算節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)64核，CPU主頻約為2.8 GHz。

3.1 并行正確性及收斂性測(cè)試

正確性測(cè)試模型為旋成體外形繞流模擬，旨在分析并行和串行計(jì)算結(jié)果的一致性。來流條件為：速度V= 1200 m/s，高度H=24 km，側(cè)滑角β= 0°，迎角α=2.0°。由于流動(dòng)對(duì)稱，宜采用半模計(jì)算，網(wǎng)格如圖7所示，網(wǎng)格單元總數(shù)約為200萬，單元類型包括四面體、棱柱和金字塔。計(jì)算核數(shù)分別為：1、 32、64、128、256、512、1024、2048、4096和8192。

圖7 旋成體外形計(jì)算網(wǎng)格Fig.7 Computing grid of revolving body

圖8～圖10分別給出了升力、阻力和俯仰力矩系數(shù)的收斂歷程。由于并行計(jì)算和串行計(jì)算迭代過程中，網(wǎng)格單元通量計(jì)算的順序并不完全相同，因此收斂過程有小量的差別，但收斂趨勢(shì)一致且收斂后的值在誤差范圍內(nèi)相同。不同核數(shù)對(duì)應(yīng)收斂后的數(shù)值列于表1，可以看出，并行計(jì)算收斂數(shù)值和串行計(jì)算偏差小于10-5。

圖8 升力系數(shù)收斂曲線Fig.8 Convergence history of lift

圖10 俯仰力矩系數(shù)收斂曲線Fig.10 Convergence history of pitch moment

CoresCoefficient of liftDeviationCoefficient of dragDeviationCoefficient of pitch momentDeviation10.0710874--0.137338---0.0433570--320.0710849-2.50×10-60.137318-2.00×10-5-0.04335691.00×10-7640.07108972.30×10-60.137325-1.30×10-5-0.0433607-3.70×10-61280.07108972.30×10-60.137325-1.30×10-5-0.0433607-3.70×10-62560.0710849-2.50×10-60.137318-2.00×10-5-0.04335691.00×10-75120.07109083.40×10-60.1373435.00×10-6-0.0433632-6.20×10-610240.07108972.30×10-60.137325-1.30×10-5-0.0433607-3.70×10-620480.0710868-6.00×10-70.137324-1.40×10-5-0.0433588-1.80×10-640960.07108982.40×10-60.137325-1.30×10-5-0.0433608-3.80×10-681920.07108982.40×10-60.137325-1.30×10-5-0.0433608-3.80×10-6Max.0.07109083.40×10-60.1373435.00×10-6-0.04335691.00×10-7Min.0.0710849-2.50×10-60.137318-2.00×10-5-0.0433632-6.20×10-6

圖11和圖12分別給出了加速比和并行效率隨計(jì)算核數(shù)增加的變化，可以看出，512核是加速比快速偏離理論值的臨界核數(shù)，大于該值后效率快速降低，主要原因?yàn)榫W(wǎng)格單元很少，并行通信占用較多的時(shí)間。

圖13和圖14分別給出了并行核數(shù)為64、256、1024和4096四種情況下密度和殘差的收斂曲線，隨著并行核數(shù)的增加，除256核外其它并行規(guī)模下殘差收斂速度基本一致。

圖11 加速比隨計(jì)算核數(shù)變化Fig.11 Speed up versus processors

圖12 并行效率隨計(jì)算核數(shù)變化Fig.12 Efficiency versus processors

圖13 密度殘差收斂曲線Fig.13 Convergence history of density residual

圖14 能量殘差收斂曲線Fig.14 Convergence history of energy residual

3.2 大規(guī)模并行效率測(cè)試

用于并行效率測(cè)試的模型為運(yùn)輸機(jī)繞流模擬，網(wǎng)格單元總數(shù)約為1.2億。64個(gè)并行分區(qū)后的表面網(wǎng)格如圖15所示，由飛機(jī)頭部和翼梢局部放大圖可以看出，計(jì)算網(wǎng)格很密，其有利于準(zhǔn)確模擬飛行阻力。來流條件為：飛行速度V=300 m/s，飛行高度H=11 km，側(cè)滑角β= 0°，迎角α= 2.0°。

圖15 運(yùn)輸機(jī)構(gòu)型并行分區(qū)后的表面網(wǎng)格Fig.15 Partitioned surface grid for transport aircraft

為了增加可靠性，效率測(cè)試每步計(jì)算時(shí)間是取迭代200步的平均時(shí)間。由于網(wǎng)格量大，單核計(jì)算時(shí)間太長(zhǎng)，因此并行加速比S和效率E測(cè)試都是以64核為基準(zhǔn)。

并行加速比隨并行核數(shù)增加變化曲線如圖16所示，可以看出，該解算器并行性能很好，加速比為線性甚至超線性。出現(xiàn)超線性的原因主要有：a) 非結(jié)構(gòu)網(wǎng)格沒有拓?fù)湎拗?，采用Metis并行分區(qū)后容易達(dá)到很高的負(fù)載平衡性能；b) 在分布式并行系統(tǒng)中，每個(gè)計(jì)算進(jìn)程擁有其單獨(dú)的緩存。隨著并行規(guī)模增加，緩存總量增加，因此每個(gè)核擁有的緩存增加，提高了緩存命中率；c) 當(dāng)并行規(guī)模增加后，每個(gè)核網(wǎng)格量減小，需要的內(nèi)存等資源減少，有利于提高并行系統(tǒng)的性能。

圖16 加速比隨計(jì)算核數(shù)變化Fig.16 Speed up versus processors

并行效率隨并行規(guī)模增加變化曲線如圖17所示，可以看出，由于加速比性能優(yōu)越，在并行核數(shù)小于14336時(shí)，并行效率都在100%以上，直到18816核，并行效率都保持在80%以上。

圖17 并行效率隨計(jì)算核數(shù)變化Fig.17 Efficiency versus processors

圖18給出了并行效率與每核網(wǎng)格單元數(shù)的關(guān)系，當(dāng)網(wǎng)格單元數(shù)太少，并行數(shù)據(jù)交換耗費(fèi)的時(shí)間相比迭代計(jì)算時(shí)間不再是小量時(shí)，并行效率隨網(wǎng)格單元數(shù)減少快速下降。不同的數(shù)值格式和并行實(shí)現(xiàn)下降的拐點(diǎn)不同，本文在每核網(wǎng)格單元數(shù)大于1萬時(shí)，都能保持很高的并行效率。

圖18 并行效率隨每核網(wǎng)格單元數(shù)變化Fig.18 Efficiency versus numbers of grid elements

4 結(jié) 論

本文發(fā)展了基于分布式并行系統(tǒng)的大規(guī)模流動(dòng)模擬方法，包括離散控制方程的緊致型數(shù)值格式、基于Metis的網(wǎng)格分區(qū)技術(shù)和并行邊界處理及虛擬單元技術(shù)等。數(shù)值結(jié)果表明，軟件并行計(jì)算與串行計(jì)算結(jié)果一致，并行效率高。

與串行計(jì)算相比，并行計(jì)算流動(dòng)收斂過程有小量的差別，但收斂數(shù)值在合理的誤差范圍內(nèi)一致。在一定的并行核數(shù)范圍內(nèi)，表現(xiàn)出超線性的加速比特性。由于加速比性能好，所有的測(cè)試算例都具有很高的并行效率。另外，要保證有較高的并行效率，每核網(wǎng)格單元數(shù)不能太少。

總之，本文方法能有效地計(jì)算復(fù)雜飛行器構(gòu)型的氣動(dòng)性能，并可充分利用大規(guī)模并行計(jì)算減少流動(dòng)模擬時(shí)間和提高計(jì)算效率。