唐靜，張健，李彬,*，崔鵬程，周乃春

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

計(jì)算流體力學(xué)(CFD)已成為空氣動(dòng)力學(xué)研究和應(yīng)用的主要手段之一。影響CFD精度的兩個(gè)主要因素是數(shù)值方法和計(jì)算網(wǎng)格。數(shù)值方法研究較多且應(yīng)用經(jīng)驗(yàn)豐富，但由于復(fù)雜工程外形的網(wǎng)格生成強(qiáng)烈依賴用戶經(jīng)驗(yàn)，計(jì)算網(wǎng)格已成為影響CFD精準(zhǔn)度的關(guān)鍵因素[1]。近年來(lái)，耦合流場(chǎng)信息、自動(dòng)調(diào)整網(wǎng)格分布的自適應(yīng)技術(shù)得到了廣泛的研究，且在NASA報(bào)告 “2030年CFD展望研究”中將其列為亟需發(fā)展的一種網(wǎng)格技術(shù)[2]。

結(jié)構(gòu)網(wǎng)格通常只能采用移動(dòng)網(wǎng)格點(diǎn)的r型自適應(yīng)技術(shù)[3-4]，優(yōu)化網(wǎng)格的能力有限。h型自適應(yīng)技術(shù)可動(dòng)態(tài)增減網(wǎng)格單元數(shù)量，能有效地優(yōu)化網(wǎng)格分布，包括笛卡爾網(wǎng)格、四面體網(wǎng)格和混合網(wǎng)格3種類型。笛卡爾網(wǎng)格可基于八叉樹(shù)的網(wǎng)格生成算法進(jìn)一步開(kāi)展網(wǎng)格自適應(yīng)[5]，但模擬附面層難度很大。四面體自適應(yīng)技術(shù)發(fā)展較為完備[6]，然而也存在附面層模擬精度不高的問(wèn)題。

非結(jié)構(gòu)混合網(wǎng)格兼顧了黏性模擬精度和易于生成的要求，廣泛應(yīng)用于黏性流動(dòng)的模擬，但臨近物面的層結(jié)構(gòu)網(wǎng)格是自適應(yīng)技術(shù)的難點(diǎn)。基于多面體單元類型的自適應(yīng)方法[7-9]具備同時(shí)改變物面法向和切向的網(wǎng)格分布能力，并且有效消除了相鄰網(wǎng)格單元間的相互影響，魯棒性較高。

為快速解決實(shí)際飛行器的氣動(dòng)問(wèn)題，CFD流場(chǎng)求解通常都采用大規(guī)模的并行技術(shù)。網(wǎng)格自適應(yīng)也須相應(yīng)地采用匹配的并行技術(shù)，否則將成為自適應(yīng)CFD模擬流程中的性能瓶頸。

CFD流場(chǎng)求解的并行算法較為成熟，較多研究機(jī)構(gòu)[10-11]和超算中心[12-13]都已實(shí)現(xiàn)了CFD流場(chǎng)求解的大規(guī)模并行。流場(chǎng)求解并行實(shí)現(xiàn)通常采用區(qū)域分解策略，即在負(fù)載平衡的約束下進(jìn)行網(wǎng)格分區(qū)。為了減少內(nèi)存的開(kāi)銷，流場(chǎng)求解時(shí)每個(gè)并行進(jìn)程只存儲(chǔ)局部區(qū)域的網(wǎng)格，各并行進(jìn)程在分區(qū)交界面附近進(jìn)行流場(chǎng)數(shù)據(jù)的通信[10]。

當(dāng)前，網(wǎng)格自適應(yīng)并行技術(shù)的研究主要集中于h型自適應(yīng)技術(shù)需要的動(dòng)態(tài)負(fù)載平衡技術(shù)，用以解決因不同進(jìn)程網(wǎng)格單元增減數(shù)量差異大造成的后續(xù)流場(chǎng)計(jì)算負(fù)載不平衡的問(wèn)題[14-15]。動(dòng)態(tài)負(fù)載平衡算法可大致分為擴(kuò)散法[16-19]和兩步法[20-26]兩種。擴(kuò)散方法采用迭代思想，每步只針對(duì)并行邊界兩側(cè)的網(wǎng)格單元實(shí)施遷移操作，直到達(dá)到全局負(fù)載平衡的指標(biāo)，一般耗時(shí)較長(zhǎng)。兩步法則是通過(guò)重分區(qū)和數(shù)據(jù)遷移先后兩個(gè)過(guò)程實(shí)現(xiàn)。根據(jù)重分區(qū)是否并行，可分為管理進(jìn)程重分區(qū)[20-22]和并行重分區(qū)[23-26]。前者須首先匯總所有進(jìn)程的網(wǎng)格數(shù)據(jù)，繼而基于串行算法進(jìn)行分區(qū)，最后將分區(qū)后的網(wǎng)格分發(fā)到各自進(jìn)程。該方法可達(dá)到很高的負(fù)載平衡性能，但由于涉及大量網(wǎng)格數(shù)據(jù)的匯總和分發(fā)，并行效率通常較低且時(shí)間開(kāi)銷較大。并行重分區(qū)則是所有進(jìn)程都采用并行分區(qū)算法(如ParMetis程序[27]等)進(jìn)行協(xié)同分區(qū)，然后將少量的網(wǎng)格數(shù)據(jù)傳遞到目標(biāo)進(jìn)程，雖然負(fù)載平衡性能略差，但通常效率更高和計(jì)算時(shí)間更少。

網(wǎng)格自適應(yīng)技術(shù)一般包括誤差估計(jì)、網(wǎng)格單元分布優(yōu)化、表面網(wǎng)格幾何投影、空間網(wǎng)格協(xié)調(diào)匹配和流場(chǎng)插值等主要步驟[9]。為了提高整個(gè)自適應(yīng)系統(tǒng)的并行效率，需要針對(duì)每個(gè)步驟涉及的自適應(yīng)技術(shù)設(shè)計(jì)相應(yīng)的并行算法。

本文基于分布式并行構(gòu)架及消息傳遞接口(MPI)，針對(duì)混合網(wǎng)格自適應(yīng)系統(tǒng)中的關(guān)鍵技術(shù)開(kāi)展并行算法研究，主要包括網(wǎng)格單元分布優(yōu)化并行技術(shù)、表面網(wǎng)格投影并行技術(shù)和空間網(wǎng)格協(xié)調(diào)并行技術(shù)，以及動(dòng)態(tài)負(fù)載平衡技術(shù)。首先，在網(wǎng)格單元分布優(yōu)化方面，提出了“先唯一后同一”的兩步法策略實(shí)現(xiàn)了并行增添網(wǎng)格對(duì)象，采用預(yù)留編號(hào)區(qū)間方法保證了新增網(wǎng)格對(duì)象編號(hào)的唯一性，并設(shè)計(jì)了網(wǎng)格單元和網(wǎng)格點(diǎn)的特征數(shù)據(jù)恢復(fù)了網(wǎng)格對(duì)象的同一性。其次，采用局部曲面擬合技術(shù)實(shí)現(xiàn)了曲面重構(gòu)和新增物理網(wǎng)格點(diǎn)投影的完全并行，基于數(shù)據(jù)通信實(shí)現(xiàn)了新增虛擬網(wǎng)格點(diǎn)的投影。第三，提出了半并行算法實(shí)現(xiàn)了空間網(wǎng)格匹配過(guò)程的并行化，快速解決了網(wǎng)格投影造成的網(wǎng)格單元交錯(cuò)問(wèn)題。第四，發(fā)展了適用于基于單元剖分自適應(yīng)方法的動(dòng)態(tài)負(fù)載平衡技術(shù)，采用基于根節(jié)點(diǎn)網(wǎng)格單元及子單元加權(quán)的方法實(shí)現(xiàn)了并行重分區(qū)，采用數(shù)據(jù)打包/解包技術(shù)提高了網(wǎng)格數(shù)據(jù)遷移的并行性能，繼而采用圓柱激波流場(chǎng)自適應(yīng)模擬驗(yàn)證了動(dòng)態(tài)負(fù)載平衡技術(shù)。最后，采用三角翼網(wǎng)格自適應(yīng)加密測(cè)試了本文建立的自適應(yīng)系統(tǒng)的并行效率。

1 網(wǎng)格自適應(yīng)系統(tǒng)的并行技術(shù)

網(wǎng)格自適應(yīng)技術(shù)通常包括網(wǎng)格單元誤差估計(jì)、自適應(yīng)網(wǎng)格優(yōu)化和網(wǎng)格流場(chǎng)插值3個(gè)部分。圖1給出了結(jié)合網(wǎng)格自適應(yīng)的CFD流場(chǎng)并行模擬的執(zhí)行流程，并行算法涉及各模塊的數(shù)據(jù)交換。其中，網(wǎng)格誤差估計(jì)和流場(chǎng)插值都是基于當(dāng)前網(wǎng)格和進(jìn)程間并行對(duì)應(yīng)關(guān)系開(kāi)展幾何量和流場(chǎng)數(shù)據(jù)的代數(shù)運(yùn)算，其并行實(shí)現(xiàn)可直接使用流場(chǎng)求解的并行算法。自適應(yīng)網(wǎng)格優(yōu)化可細(xì)分為網(wǎng)格單元分布優(yōu)化、表面網(wǎng)格幾何投影和空間網(wǎng)格協(xié)調(diào)匹配3個(gè)部分，是自適應(yīng)系統(tǒng)并行算法設(shè)計(jì)的核心。

圖1 并行網(wǎng)格自適應(yīng)數(shù)值模擬執(zhí)行流程Fig.1 Flowchart of parallel numerical simulation with mesh adaptation

在基于網(wǎng)格分區(qū)的CFD并行實(shí)現(xiàn)中，虛擬網(wǎng)格技術(shù)常用來(lái)存儲(chǔ)并行流場(chǎng)數(shù)據(jù)，以實(shí)現(xiàn)物理邊界條件和并行邊界條件的統(tǒng)一處理[10]。本文虛擬網(wǎng)格單元的模板擴(kuò)展為與網(wǎng)格點(diǎn)相連的網(wǎng)格單元，如圖2所示，實(shí)線圍成的網(wǎng)格單元為真實(shí)物理單元，虛線相關(guān)的網(wǎng)格單元為虛擬單元。

圖2 并行邊界虛擬網(wǎng)格單元Fig.2 Ghost elements on parallel interface

通常情況下，并行算法與自適應(yīng)系統(tǒng)采用的具體方法有關(guān)。本文針對(duì)非結(jié)構(gòu)混合網(wǎng)格，且主要基于文獻(xiàn)[9]中h型混合網(wǎng)格自適應(yīng)實(shí)現(xiàn)方法，分別設(shè)計(jì)了自適應(yīng)網(wǎng)格優(yōu)化3項(xiàng)技術(shù)的并行算法。另外，針對(duì)網(wǎng)格加密或稀疏造成的負(fù)載不平衡，本文在兩步法的基礎(chǔ)上設(shè)計(jì)了適合多級(jí)網(wǎng)格自適應(yīng)的動(dòng)態(tài)負(fù)載平衡算法。

1.1 網(wǎng)格單元分布優(yōu)化并行技術(shù)

采用基于單元剖分的h型自適應(yīng)方法實(shí)現(xiàn)網(wǎng)格單元分布優(yōu)化時(shí)，各進(jìn)程并行地剖分網(wǎng)格單元會(huì)同時(shí)新增網(wǎng)格單元和網(wǎng)格點(diǎn)。為了準(zhǔn)確建立進(jìn)程間的并行連接關(guān)系，新增的網(wǎng)格對(duì)象(包括網(wǎng)格單元和網(wǎng)格點(diǎn))需滿足如下兩條并行相容性準(zhǔn)則：

1) 唯一性準(zhǔn)則：真實(shí)物理網(wǎng)格對(duì)象的編號(hào)和并行進(jìn)程編號(hào)在所有進(jìn)程中具有唯一性。

2) 同一性準(zhǔn)則：虛擬網(wǎng)格對(duì)象的編號(hào)和并行進(jìn)程編號(hào)與存在于其他進(jìn)程上的真實(shí)網(wǎng)格對(duì)象的編號(hào)和并行進(jìn)程編號(hào)相同。

并行進(jìn)程編號(hào)的并行相容性實(shí)現(xiàn)相對(duì)簡(jiǎn)單，網(wǎng)格單元和網(wǎng)格點(diǎn)分別按如下方法設(shè)置：通過(guò)剖分得到的新增子網(wǎng)格單元的進(jìn)程編號(hào)與其父級(jí)網(wǎng)格單元相等；新增網(wǎng)格點(diǎn)的進(jìn)程編號(hào)等于與該網(wǎng)格點(diǎn)相連的所有網(wǎng)格單元的進(jìn)程編號(hào)的最小值。

本文虛擬網(wǎng)格單元也一并開(kāi)展剖分或粗化操作，以減少子網(wǎng)格信息的并行通信。由于新增網(wǎng)格對(duì)象的同一性準(zhǔn)則只能通過(guò)并行通信實(shí)現(xiàn)，同時(shí)考慮唯一性和同一性準(zhǔn)則會(huì)由于并行通信將并行操作退化為串行。因此，本文采用兩步法策略先后實(shí)現(xiàn)網(wǎng)格對(duì)象編號(hào)的并行相容性，即在新增網(wǎng)格對(duì)象時(shí)先保證唯一性準(zhǔn)則，完成后再統(tǒng)一恢復(fù)新增網(wǎng)格對(duì)象的同一性。

唯一性準(zhǔn)則要求為新增的網(wǎng)格對(duì)象分配全局唯一的編號(hào)。串行環(huán)境下，可以采用計(jì)數(shù)器為每個(gè)新增網(wǎng)格對(duì)象依次分配唯一編號(hào)。然而，并行環(huán)境下計(jì)數(shù)器分配編號(hào)將會(huì)出現(xiàn)沖突。本文采用預(yù)留編號(hào)區(qū)間方法實(shí)現(xiàn)并行編號(hào)。主要的思想是每個(gè)并行進(jìn)程計(jì)數(shù)器的取值限定在互不重疊的編號(hào)區(qū)間內(nèi)，只在某個(gè)并行進(jìn)程的計(jì)數(shù)器取值超出了可用區(qū)間時(shí)才需要通過(guò)并行通信確定新的可用區(qū)間。并行進(jìn)程p的可用編號(hào)區(qū)間[Imin,Imax]可確定為

(1)

式中：Nmax為當(dāng)前所有進(jìn)程的計(jì)數(shù)器的最大值；np為并行進(jìn)程總數(shù)，p=1,2,…,np；NT為計(jì)數(shù)器區(qū)間跨度，本文取512。新增網(wǎng)格點(diǎn)也可按類似方法并行地進(jìn)行編號(hào)分配。

各并行進(jìn)程采用預(yù)留編號(hào)區(qū)間方法獨(dú)立編號(hào)新增網(wǎng)格對(duì)象后，每個(gè)網(wǎng)格單元或網(wǎng)格點(diǎn)的編號(hào)都滿足全并行域內(nèi)的唯一性，如圖3所示，包括虛擬網(wǎng)格單元(圖中用字母c和d表示)和虛擬網(wǎng)格點(diǎn)(圖中用字母p表示)。

圖3 新增網(wǎng)格單元和網(wǎng)格點(diǎn)實(shí)現(xiàn)唯一性編號(hào)示意圖Fig.3 Schematic diagram of uniquely numbering for new mesh elements and points

本文同一性的恢復(fù)采用基于網(wǎng)格對(duì)象特征數(shù)據(jù)的并行通信來(lái)實(shí)現(xiàn)，包括“詢問(wèn)”和“返回”兩個(gè)主要步驟。網(wǎng)格特征數(shù)據(jù)是區(qū)分不同網(wǎng)格對(duì)象的唯一性標(biāo)識(shí)，在“詢問(wèn)”步驟通過(guò)特征數(shù)據(jù)在真實(shí)物理網(wǎng)格對(duì)象所在的進(jìn)程查詢對(duì)象編號(hào)，然后在“返回”步驟設(shè)置對(duì)應(yīng)的對(duì)象編號(hào)。

對(duì)于新增網(wǎng)格單元的特征數(shù)據(jù)De，取為父級(jí)單元編號(hào)Ip和子單元序號(hào)Sc組成的數(shù)對(duì)：

Depair〈Ip,Sc〉

(2)

對(duì)于新增網(wǎng)格點(diǎn)的特征數(shù)據(jù)Dp，取如下四元素：

(3)

新增網(wǎng)格點(diǎn)位于原網(wǎng)格單元的相對(duì)位置包括體心、面心和網(wǎng)格邊3種情況，分別設(shè)置為

1) 體心：C1為原網(wǎng)格單元編號(hào)，C2=C3=C4=Big，Big為程序語(yǔ)言允許的最大正整數(shù)，下同。

2) 面心：若為三角形面，C1、C2和C3分別取為該三角形3個(gè)頂點(diǎn)的編號(hào)，C4=Big；若為四邊形面，C1、C2、C3和C4分別取為該四邊形4個(gè)頂點(diǎn)的編號(hào)。

3) 網(wǎng)格邊：C1和C2分別取為網(wǎng)格邊兩個(gè)端點(diǎn)的編號(hào)，C3=C4=Big。

不同進(jìn)程上定義同一網(wǎng)格面或網(wǎng)格邊所采用的網(wǎng)格點(diǎn)的順序可能不同，因而同一網(wǎng)格點(diǎn)的特征數(shù)據(jù)包含的4個(gè)元素的順序也可能不同。對(duì)4個(gè)元素進(jìn)行升序排序，則兩個(gè)特征數(shù)據(jù)相等判斷只需逐一判斷每個(gè)元素。另外，雖然特征數(shù)據(jù)分為3種情況分別設(shè)置，但得到的特征數(shù)據(jù)互不相同，因而在并行算法設(shè)計(jì)時(shí)可采用同一數(shù)據(jù)結(jié)構(gòu)。圖3中新增的網(wǎng)格單元和網(wǎng)格點(diǎn)恢復(fù)編號(hào)的同一性后，如圖4所示。

圖4 新增網(wǎng)格單元和網(wǎng)格點(diǎn)編號(hào)同一性恢復(fù)示意圖Fig.4 Schematic diagram of identity recovery for new mesh elements and points

以上并行相容性算法可能會(huì)造成網(wǎng)格對(duì)象的編號(hào)不連續(xù)，包括同一進(jìn)程內(nèi)編號(hào)不連續(xù)和全并行域內(nèi)某些編號(hào)未使用。因此，在實(shí)現(xiàn)了新增網(wǎng)格對(duì)象的并行相容后，可基于并行對(duì)應(yīng)關(guān)系重新排序所有物理網(wǎng)格對(duì)象，并通過(guò)并行通信同步虛擬網(wǎng)格對(duì)象在重排序后的新編號(hào)。

1.2 表面網(wǎng)格幾何投影并行技術(shù)

當(dāng)表面網(wǎng)格投影基于CAD系統(tǒng)時(shí)，每個(gè)進(jìn)程可以并行地調(diào)用CAD系統(tǒng)進(jìn)行投影。當(dāng)采用基于曲面擬合方法重構(gòu)幾何信息時(shí)，曲面擬合方法的可并行性直接影響并行算法的設(shè)計(jì)。通常非局部的曲面擬合方法的可并行性較低，并行通信量大，并行效率低。曲面擬合采用模板少且緊致的方法，如文獻(xiàn)[9]中局部Coons曲面擬合方法，有利于降低并行算法實(shí)現(xiàn)難度和提高并行性能。

對(duì)于包含新增網(wǎng)格點(diǎn)的表面網(wǎng)格單元，局部Coons曲面擬合只需要表面網(wǎng)格單元的頂點(diǎn)和頂點(diǎn)處的法向信息，后者可采用點(diǎn)相連的網(wǎng)格面法向信息加權(quán)得到。因此，該擬合方法只涉及當(dāng)前表面網(wǎng)格單元和頂點(diǎn)相連的臨近表面網(wǎng)格單元。

根據(jù)本文虛擬網(wǎng)格單元的設(shè)置方法可知，真實(shí)的物理表面網(wǎng)格單元，如圖5中的網(wǎng)格面F，需要的所有表面網(wǎng)格都已存在于當(dāng)前的并行進(jìn)程中，因而曲面擬合不需要并行通信。同時(shí)，由于重構(gòu)的曲面為參數(shù)化曲面，新增表面網(wǎng)格點(diǎn)的投影目標(biāo)點(diǎn)也不需要并行通信。因此，物理表面網(wǎng)格單元上新增網(wǎng)格點(diǎn)投影可以完全獨(dú)立執(zhí)行。

圖5 曲面擬合需要的臨近表面網(wǎng)格單元示意圖Fig.5 Schematic diagram of linked surface elements for surface fitting

對(duì)于虛擬表面網(wǎng)格上的新增網(wǎng)格點(diǎn)，曲面擬合需要其他并行進(jìn)程上的表面網(wǎng)格信息，必須通過(guò)并行通信獲取。然而，考慮到該虛擬表面網(wǎng)格單元同時(shí)作為物理網(wǎng)格單元存在于其他并行進(jìn)程上，且新增網(wǎng)格點(diǎn)的目標(biāo)點(diǎn)也已經(jīng)求出，因此只需通過(guò)并行通信直接同步投影的目標(biāo)點(diǎn)即可。

1.3 空間網(wǎng)格匹配并行技術(shù)

空間網(wǎng)格匹配是用來(lái)解決因表面網(wǎng)格投影造成的空間網(wǎng)格單元交叉等問(wèn)題，通常采用動(dòng)網(wǎng)格技術(shù)將表面網(wǎng)格點(diǎn)的位移傳遞到內(nèi)部網(wǎng)格點(diǎn)上?？臻g網(wǎng)格匹配并行技術(shù)的性能與采用的動(dòng)網(wǎng)格技術(shù)的并行算法直接相關(guān)。

針對(duì)文獻(xiàn)[9]中的動(dòng)網(wǎng)格技術(shù)，可以采用半并行算法[28]，其主要思想是每個(gè)并行進(jìn)程都收集完整的表面網(wǎng)格點(diǎn)數(shù)據(jù)，獨(dú)立地處理各自的內(nèi)部空間網(wǎng)格點(diǎn)。

2 流場(chǎng)求解方法簡(jiǎn)介

本文流場(chǎng)求解采用課題組自主研發(fā)的MFlow軟件[29]。該軟件基于積分型式的雷諾平均Navier-Stokes方程和采用格心格式的有限體積方法。文中采用Roe格式[30]計(jì)算對(duì)流通量和采用中心差分格式計(jì)算黏性通量，采用LU-SGS方法[31]進(jìn)行方程的迭代求解，湍流模擬采用標(biāo)準(zhǔn)SA模型[32]。

3 動(dòng)態(tài)負(fù)載平衡技術(shù)

通常情況下，增添或刪除網(wǎng)格單元的操作主要集中在少數(shù)幾個(gè)并行進(jìn)程中，造成后續(xù)流場(chǎng)求解的負(fù)載嚴(yán)重不平衡，將顯著降低并行效率。本文基于并行重分區(qū)-網(wǎng)格數(shù)據(jù)遷移的兩步法思想，發(fā)展了適用于多層級(jí)網(wǎng)格自適應(yīng)[9]的動(dòng)態(tài)負(fù)載平衡技術(shù)。

在完成增添或刪除網(wǎng)格單元后，多層級(jí)網(wǎng)格單元間的連接關(guān)系存在跨層級(jí)情況，構(gòu)建連接關(guān)系的算法較為復(fù)雜。在并行環(huán)境下構(gòu)建單元連接關(guān)系，不同的網(wǎng)格單元需要不同層級(jí)的虛擬網(wǎng)格單元。為了避免虛擬網(wǎng)格單元的遴選過(guò)程帶來(lái)的計(jì)算開(kāi)銷和降低算法復(fù)雜度，本文將所有層級(jí)的虛擬網(wǎng)格單元都通過(guò)并行通信進(jìn)行傳遞，以適應(yīng)所有物理網(wǎng)格單元的連接關(guān)系的構(gòu)建。

并行重分區(qū)采用ParMetis程序庫(kù)[27]，分區(qū)所需的網(wǎng)格連接圖的構(gòu)建基于最初的根節(jié)點(diǎn)網(wǎng)格單元，根節(jié)點(diǎn)網(wǎng)格單元的權(quán)重取為該單元包含活動(dòng)子網(wǎng)格單元(流場(chǎng)求解需要)的總數(shù)。該方法能有效保證后續(xù)流場(chǎng)求解的負(fù)載平衡性能。在完成根節(jié)點(diǎn)網(wǎng)格單元的并行進(jìn)程編號(hào)設(shè)置后，所有子網(wǎng)格單元的并行進(jìn)程編號(hào)設(shè)置為對(duì)應(yīng)的根節(jié)點(diǎn)網(wǎng)格單元的進(jìn)程編號(hào)值。

網(wǎng)格重分區(qū)完成后，不再屬于當(dāng)前進(jìn)程的網(wǎng)格單元，以及臨近的虛擬網(wǎng)格單元，都需要通過(guò)并行通信遷移到目標(biāo)進(jìn)程。網(wǎng)格單元的數(shù)據(jù)包含整數(shù)和實(shí)數(shù)等多種數(shù)據(jù)類型，且通常較為分散地存儲(chǔ)在內(nèi)存中。為了減少并行通信函數(shù)的啟動(dòng)頻率，提高通信效率，本文采用數(shù)據(jù)打包/解包的方法，如圖6所示，將多種數(shù)據(jù)類型、分散的網(wǎng)格數(shù)據(jù)轉(zhuǎn)為連續(xù)存儲(chǔ)，增大并行通信的數(shù)據(jù)粒度。

圖6 網(wǎng)格單元數(shù)據(jù)打包/解包示意圖Fig.6 Schematic diagram of data packing/unpacking for mesh elements

此處采用二維圓柱超聲速繞流算例，開(kāi)展基于梯度誤差估計(jì)[9]的網(wǎng)格自適應(yīng)流動(dòng)模擬，以初步驗(yàn)證本文建立的動(dòng)態(tài)負(fù)載平衡技術(shù)。來(lái)流條件為：馬赫數(shù)為2.0，靜壓為30 800 Pa，靜溫為230 K。采用4核并行計(jì)算，網(wǎng)格自適應(yīng)共迭代3次。圖7對(duì)比了初始網(wǎng)格和自適應(yīng)3次后計(jì)算的流場(chǎng)，自適應(yīng)后激波和膨脹波都捕捉得更精細(xì)。自適應(yīng)網(wǎng)格負(fù)載平衡的動(dòng)態(tài)變化過(guò)程見(jiàn)圖8，圖中同一顏色的網(wǎng)格歸屬相同的并行進(jìn)程?？梢钥闯觯瑸榱吮３趾罄m(xù)流場(chǎng)求解的負(fù)載平衡，網(wǎng)格并行分區(qū)的交界面隨著自適應(yīng)迭代動(dòng)態(tài)變化。另外，改變歸屬進(jìn)程的網(wǎng)格單元都集中于并行交界面附近的位置，有利于降低網(wǎng)格單元數(shù)據(jù)遷移的通信開(kāi)銷。

圖7 自適應(yīng)前后流場(chǎng)對(duì)比Fig.7 Comparison of flow fields before and after adaptation

圖8 自適應(yīng)負(fù)載平衡并行網(wǎng)格動(dòng)態(tài)變化過(guò)程Fig.8 Dynamical variation of parallel mesh during mesh adaptation for load balance

4 網(wǎng)格自適應(yīng)并行技術(shù)效率測(cè)試

為了測(cè)試本文自適應(yīng)系統(tǒng)的并行技術(shù)在大規(guī)模并行計(jì)算中的并行效率，采用尖前緣三角翼模型。同時(shí)，對(duì)自適應(yīng)系統(tǒng)3個(gè)主要模塊的并行性能也進(jìn)行了統(tǒng)計(jì)，包括網(wǎng)格分布優(yōu)化、物面幾何投影和空間網(wǎng)格匹配。測(cè)試來(lái)流條件為：雷諾數(shù)為1.38×107，靜溫為322 K，馬赫數(shù)為0.4，迎角為13.3°，側(cè)滑角為0°。由于來(lái)流對(duì)稱，采用包含對(duì)稱面的半模網(wǎng)格開(kāi)展流動(dòng)模擬，初始網(wǎng)格約2 500萬(wàn)，網(wǎng)格單元包含四面體、金字塔和三棱柱3種類型。網(wǎng)格自適應(yīng)采用均勻加密所有單元，自適應(yīng)后網(wǎng)格單元總量約為2.0億。測(cè)試的并行核數(shù)分別為：16、32、64、128、256和512。由于計(jì)算機(jī)內(nèi)存限制，此處以16核作為并行效率計(jì)算的基準(zhǔn)，并行加速比S和并行效率E計(jì)算方法分別為

(4)

式中：T為一定迭代步數(shù)計(jì)算的平均時(shí)間；n為并行核數(shù)；下標(biāo)0表示相應(yīng)的基準(zhǔn)值。

在128核以內(nèi)，整個(gè)網(wǎng)格自適應(yīng)系統(tǒng)和3個(gè)主要模塊的并行加速比差別較小，如圖9所示。更大并行規(guī)模時(shí)，網(wǎng)格分布優(yōu)化和表面網(wǎng)格投影兩個(gè)模塊的加速比優(yōu)于整個(gè)自適應(yīng)系統(tǒng)，空間網(wǎng)格匹配模塊加速比最低。該模塊加速比性能最低的原因是本文采用半并行算法，每個(gè)并行進(jìn)程都針對(duì)所有表面網(wǎng)格點(diǎn)開(kāi)展插值運(yùn)算，降低了該并行算法在較大并行規(guī)模時(shí)的并行性能。

圖9 自適應(yīng)系統(tǒng)和3個(gè)模塊的并行加速比Fig.9 Parallel speed-up ratio of adaption system and its three kernel modules

圖10給出了整個(gè)網(wǎng)格自適應(yīng)系統(tǒng)和3個(gè)主要模塊的并行效率，整體來(lái)看，128核以內(nèi)自適應(yīng)系統(tǒng)的并行效率大于80%，表明自適應(yīng)系統(tǒng)有較高的可擴(kuò)展性。更大并行規(guī)模時(shí)，空間網(wǎng)格匹配的并行效率下降最快，而表面網(wǎng)格投影模塊并行效率最高。由于網(wǎng)格投影模塊中曲面擬合和物理網(wǎng)格點(diǎn)投影都不需要并行通信，僅需要并行同步虛擬網(wǎng)格點(diǎn)，因而該模塊的并行效率很高。

圖10 自適應(yīng)系統(tǒng)和3個(gè)模塊的并行效率Fig.10 Parallel efficiency of adaption system and its three kernel modules

不同并行核數(shù)時(shí)，網(wǎng)格自適應(yīng)系統(tǒng)、3個(gè)主要模塊和并行通信耗費(fèi)的總時(shí)間及比例見(jiàn)表1，3個(gè)主要模塊中，表面網(wǎng)格投影模塊計(jì)算時(shí)間最少，小于總時(shí)間的1%，網(wǎng)格分布優(yōu)化模塊計(jì)算時(shí)間最長(zhǎng)。圖11給出了耗時(shí)比例隨并行規(guī)模增加的變化，表面網(wǎng)格投影模塊耗時(shí)比例略有減小，網(wǎng)格分布優(yōu)化模塊耗時(shí)比例逐漸減小，空間網(wǎng)格匹配模塊略有增加。另外，并行通信時(shí)間所占比例隨并行規(guī)模較為快速地增加，表明該網(wǎng)格自適應(yīng)系統(tǒng)屬于通信密集型并行系統(tǒng)。

圖12給出了采用256核求解流場(chǎng)的密度殘差和升力系數(shù)的收斂曲線，迭代約11 000步后殘差下降約3個(gè)量級(jí)，氣動(dòng)力獲得收斂解。此時(shí)，流場(chǎng)求解共需約3.2×104s(約8.9 h)。使用256核完成網(wǎng)格自適應(yīng)共需200 s，約為流場(chǎng)求解的0.6%。因此，盡管256核時(shí)網(wǎng)格自適應(yīng)并行效率不高，但從總時(shí)間開(kāi)銷看是完全可接受的。

表1 自適應(yīng)系統(tǒng)及其3個(gè)主要模塊和并行通信的時(shí)間耗費(fèi)Table 1 Time cost for adaptation system, its three kernel modules, and parallel communication

圖11 自適應(yīng)3個(gè)模塊和并行通信耗時(shí)比例Fig.11 Time cost ratio of three kernel modules and parallel communication

圖12 密度殘差和升力系數(shù)收斂曲線Fig.12 Convergence history of density residual and lift coefficient

5 結(jié) 論

針對(duì)非結(jié)構(gòu)混合網(wǎng)格，設(shè)計(jì)了網(wǎng)格單元分布優(yōu)化、表面網(wǎng)格幾何投影和空間網(wǎng)格匹配相關(guān)的并行算法。同時(shí)，為后續(xù)流場(chǎng)模擬建立了基于并行重分區(qū)-網(wǎng)格數(shù)據(jù)遷移的動(dòng)態(tài)負(fù)載平衡技術(shù)。最后，采用三角翼模型測(cè)試了不同并行規(guī)模下3項(xiàng)并行算法的并行性能。

1) 提出的“先唯一后同一”兩步策略有效實(shí)現(xiàn)了網(wǎng)格單元分布優(yōu)化模塊中網(wǎng)格單元和網(wǎng)格點(diǎn)的并行增添或刪除，確保了并行環(huán)境下網(wǎng)格的并行相容性。

2) 采用的局部曲面擬合技術(shù)僅需同步虛擬網(wǎng)格點(diǎn)的信息，物理網(wǎng)格點(diǎn)投影每個(gè)進(jìn)程完全獨(dú)立，該方法具備很高可并行性。

3) 建立的并行自適應(yīng)系統(tǒng)在大規(guī)模并行時(shí)并行效率有所降低，但對(duì)于定常流動(dòng)模擬，與流場(chǎng)求解相比，總時(shí)間消耗很少，該系統(tǒng)的并行性能完全可以接受。

本文提出的兩步法策略和采用局部曲面擬合技術(shù)的思想，對(duì)其他的網(wǎng)格自適應(yīng)方法也具有一定的借鑒意義。