基于非結構網格流場超大規(guī)模并行計算

, , ,*,

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

0 引 言

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

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

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

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

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

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

1 數(shù)值方法

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

(1)

對流通量可以表示為：

(2)

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

V=nxu+nyv+nzw

(3)

H為焓值，定義為：

(4)

式中，e為流動內能。

黏性通量相對比較復雜，本文略去，具體內容可以參考文獻[17]。

在網格單元內積分，得到離散形式的控制方程：

(5)

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

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

格心格式的有限體積方法，流場未知變量存儲在網格單元的形心上，如圖1所示的三角形形心點J1、J2和J3。然而，對流通量和黏性通量計算卻需要網格單元對應的網格面心上的流場值，如K1、K2和K3。只有采用一階格式時，面上的流場值等于形心值。當采用緊致型的二階格式時，網格面心上左右兩側的流場值由同側的網格單元體心值構造。計算黏性通量時，面心流場值采用中心格式，以K1面為例，表示為:

(6)

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

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

Hv,K1=f(QK1)

(7)

黏性通量計算的具體表達式參考文獻[21]。

采用Roe格式時，對流通量可以表示為

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

(8)

式中，QK1,L和QK1,R分別為面K1左右兩側的流場值，通過各自同側的體心值構造

(9)

采用Roe格式計算對流通量的具體表達式參考文獻[20]。

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

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

(10)

則式(5)表示為:

(11)

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

(12)

式中：

(13)

通過對殘差項的線化，可以得到

(14)

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

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

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

2.1 網格分區(qū)技術

非結構網格分區(qū)方法采用基于圖論算法的Metis程序庫，該程序庫使用多級二分法或多級多支分區(qū)算法，如圖2所示。

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

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

(15)

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

圖3 非結構網格連接關系CSR表示Fig.3 Graph data in CSR format of unstructured grid

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

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

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

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

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

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

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

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

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

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

3 并行測試結果

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

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

(16)

式中：n為測試并行核數(shù)，T為每步計算平均時間；n0為并行測試基數(shù)核數(shù)，T0為基數(shù)核數(shù)對應的每步計算平均時間。

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

3.1 并行正確性及收斂性測試

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

圖7 旋成體外形計算網格Fig.7 Computing grid of revolving body

圖8～圖10分別給出了升力、阻力和俯仰力矩系數(shù)的收斂歷程。由于并行計算和串行計算迭代過程中，網格單元通量計算的順序并不完全相同，因此收斂過程有小量的差別，但收斂趨勢一致且收斂后的值在誤差范圍內相同。不同核數(shù)對應收斂后的數(shù)值列于表1，可以看出，并行計算收斂數(shù)值和串行計算偏差小于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分別給出了加速比和并行效率隨計算核數(shù)增加的變化，可以看出，512核是加速比快速偏離理論值的臨界核數(shù)，大于該值后效率快速降低，主要原因為網格單元很少，并行通信占用較多的時間。

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

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

圖12 并行效率隨計算核數(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ī)模并行效率測試

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

圖15 運輸機構型并行分區(qū)后的表面網格Fig.15 Partitioned surface grid for transport aircraft

為了增加可靠性，效率測試每步計算時間是取迭代200步的平均時間。由于網格量大，單核計算時間太長，因此并行加速比S和效率E測試都是以64核為基準。

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

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

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

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

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

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

4 結 論

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

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

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