CCFD重疊網(wǎng)格并行算法設計和優(yōu)化

劉夏真，袁 武，馬文鵬，胡曉東，陸忠華，張 鑒

(1.中國科學院計算機網(wǎng)絡信息中心，北京 100190;2.中國科學院大學，北京 100049;3.信陽師范學院計算機與信息科技學院,河南 信陽 464000)

1 引言

重疊網(wǎng)格(Overlapping Grid)[1]的基本思想是將計算區(qū)域劃分為多個相對簡單的子區(qū)域獨立生成網(wǎng)格，子區(qū)域之間存在相互覆蓋關系，重疊區(qū)域的流場信息通過插值進行匹配和耦合。重疊網(wǎng)格放寬了網(wǎng)格拓撲要求，能在剛性運動中保留初始網(wǎng)格品質(zhì)，故而被廣泛應用于研究諸如飛行器拋殼拋彈、級間分離、無人機旋翼、高鐵會車、船舶操縱運動等[2 - 5]具有相對運動為主要特征的多體流體動力學問題，顯示出獨到優(yōu)勢。

迄今為止，國外發(fā)展的重疊網(wǎng)格程序有Pegasus[6]、SUGGAR[7]、DCF3D[8]、OVERGRID[9]、BEGGAR[10]等，在工程領域廣泛使用的一些商業(yè)軟件(如CFD++、FASTRAN等)也具備網(wǎng)格重疊功能。國內(nèi)方面，中國空氣動力研究與發(fā)展中心[11]、北京航空航天大學[12,13]、西北工業(yè)大學[14,15]、南京航空航天大學[16]等單位均發(fā)展了各自的算法和軟件。早期的程序主要是串行計算代碼，隨著高性能計算在計算流體力學CFD(Computational Fluid Dynamics)領域的普遍應用，適應分布式環(huán)境的重疊網(wǎng)格裝配算法受到了廣泛重視[17,18]。

中國科學院計算機網(wǎng)絡信息中心在國家和院所項目的支持下,在重疊網(wǎng)格并行算法和優(yōu)化上開展了一些工作[19 - 21]，開發(fā)了計算流體力學并行計算軟件CCFD(China CFD)。本文將簡要介紹CCFD并行重疊網(wǎng)格方法，包括并行算法設計、負載平衡和通信優(yōu)化等。

2 重疊網(wǎng)格的基本概念

首先，簡要介紹重疊網(wǎng)格的基本概念，以空間中五球體為例，說明建立結構重疊網(wǎng)格的主要步驟[12]。

(1)對每個球體分別生成彼此覆蓋的初始網(wǎng)格，如圖1a所示。

(2)落入物面內(nèi)部的網(wǎng)格點在流場計算中無實際意義，因此需要將這部分網(wǎng)格點屏蔽掉，一般記為“洞內(nèi)點”，以區(qū)別于參與流場計算的網(wǎng)格點即“洞外點”。人們將標識洞內(nèi)點的過程形象地稱為“挖洞”，挖洞的結果即產(chǎn)生緊密圍繞洞內(nèi)點的初始洞邊界，如圖1b所示。

(3)挖洞產(chǎn)生的初始洞邊界一般貼近物面，使網(wǎng)格重疊區(qū)域龐大，同時插值區(qū)流場梯度很大，一般在挖洞結束后需要對洞面進行一些優(yōu)化，以改善洞面質(zhì)量。洞面優(yōu)化結束后，洞邊界和緊鄰的一層洞內(nèi)點轉(zhuǎn)換為插值邊界，用于各子區(qū)網(wǎng)格的信息交換，如圖1c所示。

Figure 1 Schematic diagram of the overlapping process of five spheres in space圖1 空間五球體重疊過程示意[12]

由上述內(nèi)容可知，重疊網(wǎng)格研究的主要內(nèi)容就是建立網(wǎng)格重疊關系，因此在分布式環(huán)境下建立高效、魯棒的挖洞和洞面優(yōu)化算法，是重疊網(wǎng)格并行算法設計關注的重點。此外，基于結構重疊網(wǎng)格體系的CFD求解器，在計算負載和數(shù)據(jù)通信上，也和一般結構網(wǎng)格求解器有所區(qū)別，結合重疊網(wǎng)格算法特點，設計良好的負載平衡模型和通信模式，是程序并行優(yōu)化需要關注的重點。

3 網(wǎng)格并行裝配算法

3.1 局部洞映射模型

常見的挖洞方法[22,23]有點矢法、射線求交法、洞映射方法、Object X-Ray、叉樹挖洞方法等，其中，洞映射方法通過構建輔助的直角笛卡爾網(wǎng)格來近似物體的挖洞曲面，將點與曲面之間的關系轉(zhuǎn)化為點與笛卡爾單元之間的簡單關系，因此效率和自動化程度都很高，在Pegasus 5[24]中得到了應用。

傳統(tǒng)的洞映射方法以整體物面的包圍盒作為笛卡爾邊界，這是比較直接的串行思維。在并行環(huán)境下，傳統(tǒng)方法需要每個處理器搜集和處理所有幾何體的全局物面信息，這種做法導致數(shù)據(jù)局部性較差、數(shù)據(jù)擴展性不強。CCFD設計了一種具有局部數(shù)據(jù)特性的洞映射方法，如圖2所示，對幾何體與其它網(wǎng)格塊相交的局部物面構建笛卡爾近似，局部物面是分散于各個進程的、不封閉的物面，該局部模型可以在單獨進程上進行填充和標識屬性，而不依賴于其它進程，非常適合在分布式環(huán)境下計算。

Figure 2 Local hole mapping model圖2 局部洞映射模型

局部洞映射模型是否合理，需要考慮物面不封閉引起的笛卡爾單元屬性標識問題。傳統(tǒng)方法在區(qū)分洞映射模型單元屬性時，首先要確定“種子”單元屬性，例如標記4個角區(qū)位置的單元為外部單元，再由角區(qū)向網(wǎng)格域內(nèi)推進和填充。局部模型則由于物面的不連續(xù)、不封閉，缺少類似的先驗信息。

考慮到洞映射模型作為自身計算域的映射，其內(nèi)部單元對應計算域的物體內(nèi)部，外部單元對應計算域的網(wǎng)格部分。因此，對于各個被分割的互不連通的洞外區(qū)域，如圖3所示，笛卡爾網(wǎng)格的外部單元必然落在計算網(wǎng)格中，反之，計算網(wǎng)格也必然落在非內(nèi)部單元(外部單元或邊緣單元)，這是一個充分必要條件。

Figure 3 Relationship between local hole mapping model and computing grids圖3 局部洞映射模型與計算網(wǎng)格關系

根據(jù)上述推斷，CCFD設計了一種適用于局部洞映射模型的標識方法，如圖4所示，首先，標記笛卡爾網(wǎng)格與物面相交的部分為邊緣單元，見圖4a。其次，對落入洞映射邊界的計算網(wǎng)格進行下標計算，獲得所有包含自身網(wǎng)格點的笛卡爾單元，這些笛卡爾單元必然是外部單元或邊緣單元，且分布在各個分割區(qū)域，見圖4b。最后，再依據(jù)外部單元的相鄰單元一定是非內(nèi)部單元進行遞歸掃描，即可識別所有外部單元。由于物面或局部連續(xù)或終止于邊界，外部單元的掃描過程不會穿越邊緣單元進入洞內(nèi)區(qū)域，因此外部單元的掃描結束后剩余的未識別單元就是內(nèi)部單元，見圖4c和圖4d。

Figure 4 Mesh attribute identification method for local hole mapping model圖4 局部洞映射模型的單元屬性標識方法

3.2 隱式洞面優(yōu)化

從幾何圖形學的角度看，洞邊界的建立具有一定的隨意性，只要能將落入其它物體內(nèi)部的網(wǎng)格點“挖掉”就能達到屏蔽流場信息的目的。但是，對流場計算而言，如果直接將挖洞面作為洞邊界，將使插值邊界十分靠近物面的高梯度流動區(qū)域，嚴重影響流場解算精度。因此，需要在挖洞結束后對初始洞邊界進行優(yōu)化，常見的洞面優(yōu)化技術[22,25]有陣面推進方法、割補法和隱式切割技術等。其中，陣面推進方法和割補法都是建立在迭代優(yōu)化思想的基礎上，適合于串行計算，算法并發(fā)度低，不利于并行擴展和應用。

隱式切割技術基于網(wǎng)格密度準則，對重疊區(qū)域內(nèi)所有網(wǎng)格單元進行對比和判斷，在一個輪次中完成網(wǎng)格重疊關系的建立，沒有反復的迭代過程，因此十分適合并行計算。常用的網(wǎng)格密度判斷準則包括網(wǎng)格單元體積、網(wǎng)格到物面的距離[25]等。CCFD設計了一種結合網(wǎng)格體積和壁面距離組合加權模型，如式(1)所示：

Q=1/(va·db)

(1)

其中，v為網(wǎng)格單元體積；d為網(wǎng)格到物面的距離；參數(shù)a和b可調(diào)節(jié)，比如選取a=1,b=0時即僅考慮網(wǎng)格尺度效應，選取a=0，b=1時即僅考慮物面距離的影響，不同參數(shù)的組合為設計人員調(diào)節(jié)網(wǎng)格品質(zhì)提供了一個干預手段。

傳統(tǒng)的隱式方法沒有顯式地區(qū)分洞內(nèi)點和洞外點，認為無效的插值會被洞邊界隔離，但必須要求洞邊界連續(xù)、封閉。這一要求在處理縫隙、尖角等復雜形狀或者網(wǎng)格質(zhì)量較差時難以得到保證，有可能出現(xiàn)“挖進物面”的問題，即由于洞邊界不連續(xù)，導致洞內(nèi)區(qū)域和洞外區(qū)域錯誤地出現(xiàn)信息傳遞。CCFD采取的策略是，使用魯棒性很高且并行性良好的局部洞映射方法進行初始挖洞，再使用隱式洞面優(yōu)化方法對重疊關系進行調(diào)整和優(yōu)化。CCFD將洞內(nèi)點區(qū)分為2類，如圖5所示，Ⅰ型洞內(nèi)點是由挖洞產(chǎn)生的，落入物面內(nèi)部的點；Ⅱ型洞內(nèi)點是洞面優(yōu)化過程中產(chǎn)生的被覆蓋區(qū)域。前者被物面屏蔽，不在流場中，故不允許提取流場信息，這樣就避免了出現(xiàn)無效插值點，完全克服了“挖進物面”的問題。

Figure 5 Two types of overlapping hole generated by flaps digging holes on the main wing grid圖5 襟翼對主翼網(wǎng)格挖洞產(chǎn)生的2類洞內(nèi)點

4 負載平衡和通信優(yōu)化

4.1 兩級負載平衡模型

計算網(wǎng)格的區(qū)域分解，是CFD并行計算的前提條件。常用的區(qū)域分解方法有遞歸二分法、平衡切割樹法[26]、貪婪算法[27]、KL(Kerninghan-Lin)算法、遺傳算法[28]等。結構網(wǎng)格劃分的算法，主要是考慮了計算負載在各個處理器上的平衡，當通信時間的占比隨并行規(guī)模增長而增加時，這種忽略通信開銷的做法存在不足。特別是引入重疊網(wǎng)格后，跨結點的尋點、插值等操作，將加劇通信延遲的問題。

CCFD在結構網(wǎng)格區(qū)域分解中使用了多約束條件K-way圖剖分算法[29]，設計了一種兩級負載平衡模型[30]，同時將計算量和通信量作為衡量標準。首先，對原始非均勻網(wǎng)格塊進行網(wǎng)格塊分區(qū)，得到一系列規(guī)模更小、更均衡的子網(wǎng)格塊集合，使用的方法包括遞歸邊界二分法、固定步長法等；然后基于塊分區(qū)后的子網(wǎng)格，建立一個以網(wǎng)格塊為頂點、塊內(nèi)網(wǎng)格量和塊間相鄰面為頂點權重和邊權重的有向加權圖,基于此加權圖進行圖剖分操作，得到的每個集合的頂點權值之和近似，集合間的邊權值之和均衡。將網(wǎng)格塊集合分配到不同的計算單元上，實現(xiàn)網(wǎng)格塊在處理器上的計算量和通信量負載均衡。

(2)

其中，N是幾何體數(shù)量，ni是第i個幾何體的網(wǎng)格塊數(shù)量。

Figure 6 Directed graph of multi-block structure overlapping grid圖6 多塊結構重疊網(wǎng)格有向圖

綜合考慮有向圖DG=V,E的通信量：

(3)

基于有向圖DG=V,E的設計，CCFD針對多塊結構重疊網(wǎng)格系統(tǒng)上提出的兩級負載平衡模型如圖7所示。在粗層，對原始非均勻大小網(wǎng)格塊進行塊分區(qū)處理，形成一系列規(guī)模更均勻、塊數(shù)更多的子網(wǎng)格集合；在細層，對塊分區(qū)后的子網(wǎng)格建立一個以計算量和通信量為頂點權重和邊權重的有向加權圖，使用多層圖剖分算法實現(xiàn)加權圖的分區(qū)，得到計算量和通信量均衡的集合，將集合映射到不同的處理器上，最終實現(xiàn)網(wǎng)格塊在處理器上的計算量和通信量負載均衡。

Figure 7 Two-level load balancing model for multi-block overlapping grid圖7 多塊結構重疊網(wǎng)格的兩級負載平衡模型

4.2 基于塊的通信模式

多塊結構重疊網(wǎng)格的通信比一般的結構網(wǎng)格更為復雜，不但有對接型網(wǎng)格塊對的邊界通信，還有在不同進程、相互重疊的網(wǎng)格塊對的插值數(shù)據(jù)交換，因此需要結合算法特點考慮通信優(yōu)化問題。Djomehri等[31]在對重疊網(wǎng)格求解器OVERFLOW的通信優(yōu)化中，證實了采用非阻塞通信代替阻塞通信，可以有效降低插值數(shù)據(jù)交換的頻次。在一個處理器負責多個網(wǎng)格塊的情況下，以處理器為單位進行通信就必須等到一個處理器上所有網(wǎng)格塊都計算完成后才能發(fā)起通信，這樣通信過程獨立于流場計算，通信時間無法得到隱藏，可能存在并行擴展性瓶頸。

CCFD結合應用特點對不同類型的進程間通信進行了優(yōu)化，在計算過程中，存在諸如數(shù)據(jù)量小而頻繁的數(shù)據(jù)通信、數(shù)據(jù)量大而次數(shù)固定的網(wǎng)格邊界更新、一對多的廣播數(shù)據(jù)、多對一的主從通信等，需要針對性地采用靈活多變的通信手段，包括集合通信、非阻塞通信等，以減少通信頻次和數(shù)據(jù)量。此外，由于負載不均衡引起的進程等待問題也將產(chǎn)生很大的通信時間開銷，CCFD采用了計算和通信重疊的方法，如圖8所示，數(shù)據(jù)的發(fā)送和接收由后臺處理，只在使用的時候才進行通信完成檢測，這種方法可以有效隱藏通信時間。

Figure 8 Calculation and communication overlap mode圖8 計算與通信重疊模式

重疊網(wǎng)格的并行插值計算具有局部性，即一個網(wǎng)格塊需要的插值信息僅依賴于數(shù)個網(wǎng)格塊，在同一時間步內(nèi)，任何一個網(wǎng)格塊的數(shù)據(jù)更新后都可以將更新信息發(fā)送出去，以供插值依賴塊在下一個時間步來使用。因此，CCFD設計了以塊為單位的重疊網(wǎng)格插值策略：每個處理器依次計算、更新若干網(wǎng)格塊，在當前塊的流場信息計算完畢后立即以非阻塞的通信發(fā)送給所插值依賴塊所在的處理器，此部分通信的開銷可以被該處理器上其他網(wǎng)格塊的流場計算所隱藏。與以處理器為單位的通信模式不同，以塊為單位的通信模式可能存在一個處理器與另一個處理器重復多次通信，但是可以實現(xiàn)通信和計算的完全異步，在大規(guī)模并行計算時能夠有效隱藏通信開銷。圖9展示了Processor #0中的A和J網(wǎng)格塊，Processor #1中的C和G網(wǎng)格塊向Processor #2中的B、F、I網(wǎng)格塊發(fā)送插值信息的過程。Processor #2按照信息來源的網(wǎng)格塊號升序分配、組織接收緩存區(qū)，Processor #0和Processor #1中提供插值的網(wǎng)格塊，按照目標塊所在的處理器升序組織發(fā)送緩存區(qū)，保證同一網(wǎng)格塊提供給同一處理器的數(shù)據(jù)連續(xù)存放。

Figure 9 Block communication mode:reorganization,sending,receiving processes圖9 塊通信模式：重組、發(fā)送、接收過程

5 測試與討論

5.1 雙機翼模型測試

通過開展雙機翼模型在不同并行規(guī)模下的重疊網(wǎng)格測試，研究本文方法的計算和通信負載均衡特性。該模型采用2個平行布置的ONERA-M6機翼，圖10所示為串行計算使用的粗網(wǎng)格和重疊結果示意，雙機翼網(wǎng)格對稱，插值邊界居中，重疊區(qū)域整潔。采用CCFD軟件包提供的網(wǎng)格劃分工具SPLITTER，對串行計算的粗網(wǎng)格進行加密和分區(qū)，生成的并行計算網(wǎng)格規(guī)模約6 740萬，網(wǎng)格塊數(shù)量為5 080塊，并行測試規(guī)模為64～1024核。測試環(huán)境是中國科學院計算機網(wǎng)絡信息中心的高性能集群，每個結點配置2顆Intel Xeon E5520 CPU處理器，編譯環(huán)境為Intel Fortran+OpenMPI。

Figure 10 Two-wing model and overlapping grids圖10 雙機翼模型和重疊網(wǎng)格示意

圖11是雙機翼模型測試的計算負載均衡因子CalFactor和通信負載均衡因子ComFactor隨進程數(shù)變化情況。負載均衡因子定義為目標量的最大值/平均值，該比值越大，表示目標量均衡性越差。由圖11可知，在進程數(shù)為64時，CalFactor和ComFactor分別是1.035和1.284，可知CCFD提出的兩級負載平衡模型，通過考慮網(wǎng)格的重疊關系，使存在重疊關系的網(wǎng)格塊在若干結點內(nèi)聚集以減少通信開銷，從而可獲得良好的計算和通信負載均衡。當進程數(shù)從64增長到256時(以64核為基數(shù)擴展4倍)，CalFactor和ComFactor分別增長了1.294%和2.522%；進程數(shù)從256增長到1 024時(以256核為基數(shù)擴展4倍)，CalFactor和ComFactor分別增長了11.182%和7.040%，可見通過本文的方法，計算與通信負載在強擴展條件下仍然能維持較好的均衡水平。

5.2 機翼掛載分離測試

機翼掛載分離標模[32]被廣泛用于考核非定常數(shù)值模擬方法和重疊網(wǎng)格方法。本文使用CCFD并行求解器研究了這一問題，該算例的粗網(wǎng)格規(guī)模約480萬，分為22個子區(qū)域并行。計算條件和數(shù)值方法包括：來流馬赫數(shù)0.95，模擬大氣高度7.92 km，單位雷諾數(shù)7.874×106/m；層流計算，空間離散采用Roe’s FDS方法，非定常時間推進采用雙時間步方法，時間間隔取Δt=0.65E-02。典型時刻流場和計算結果比較分別如圖12和圖13所示。圖13是分離時間0.3 s內(nèi)的3方向位移和姿態(tài)對比，圖中實線是CFD計算數(shù)據(jù)，虛線是風洞實驗數(shù)據(jù)，橫軸是有量綱的時間軸，計算值和實驗值符合良好。

創(chuàng)新理念和技術的缺乏，直接影響到了獸醫(yī)行業(yè)的發(fā)展。對獸醫(yī)技術進行創(chuàng)新，既能夠提高企業(yè)效益，實現(xiàn)企業(yè)可持續(xù)發(fā)展的目標，又可以增強企業(yè)的綜合競爭力，為新技術、新產(chǎn)品的研發(fā)提供動力，使企業(yè)與市場需求高度契合。另外，創(chuàng)新獸醫(yī)技術，還與全球經(jīng)濟變革所提出的要求間存在著一定的聯(lián)系，只有對技術、產(chǎn)品進行創(chuàng)新，才能緊跟經(jīng)濟發(fā)展的腳步，使相關產(chǎn)業(yè)的積極作用得以充分呈現(xiàn)。

Figure 12 Wing airborne separation model and typical moment flow fields圖12 機翼掛載分離模型和典型時刻流場

Figure 13 Comparison of calculation and experiment value圖13 計算和實驗值比較

采用SPLITTER對在粗網(wǎng)格進行加密和分區(qū)，生成的并行計算網(wǎng)格規(guī)模約5 900萬，網(wǎng)格塊數(shù)量為6 940塊(機翼網(wǎng)格為5 610塊，掛載網(wǎng)格為1 330塊)。擴展性測試規(guī)模為64～1024核，計算并統(tǒng)計了55個非定常時間步的結果，重疊網(wǎng)格裝配過程在每個非定常時間步內(nèi)重復進行。

圖14所示是統(tǒng)計的重疊網(wǎng)格裝配平均時間以及占非定常時間步的平均比值。由圖14可知，重疊網(wǎng)格方法應用在多體相對運動問題仿真時，網(wǎng)格處理時間遠少于流場求解時間，負載平衡模型應主要考慮流場計算部分。圖15所示是并行計算加速比隨進程數(shù)變化情況，以64核為基準的千核級并行效率達到了54.4%，說明CCFD在處理這類涉及非定常模擬和重疊網(wǎng)格計算的復雜多體相對運動仿真時，具有良好的并行效率和可擴展性。

Figure 14 Ratio of assembly time of overlapping mesh to unsteady calculation圖14 重疊網(wǎng)格裝配時間占非定常計算的比值

Figure 15 Parallel speedup ratio change with the number of processes圖15 并行加速比隨進程數(shù)變化

6 結束語

本文介紹了并行計算流體力學軟件CCFD在重疊網(wǎng)格方法上的高效并行實現(xiàn)，包括：結合算法特點發(fā)展了具有局部數(shù)據(jù)特性的新型洞映射模型和適用于分布式環(huán)境的隱式洞面優(yōu)化方法，設計了兼顧計算量和通信量的兩級負載平衡模型，建立了基于網(wǎng)格塊的通信模式并進行了調(diào)優(yōu)。數(shù)值模擬結果表明，該方法對求解有相對運動的多體氣動問題，具有良好的并行效率和可擴展性。

目前，國內(nèi)高性能計算硬件平臺發(fā)展迅猛，而相應的并行軟件稍顯不足。對于國產(chǎn)平臺，不同硬件架構對同一模型算法的并行實現(xiàn)效率差別很大，并行軟件只有協(xié)同國產(chǎn)平臺的硬件環(huán)境開展設計和優(yōu)化，才能充分發(fā)揮硬件性能。在下一步工作中，將瞄準國產(chǎn)E級異構平臺，結合硬件特點開展算法的并行設計、代碼移植和深層優(yōu)化。