非結(jié)構(gòu)混合網(wǎng)格自適應并行技術(shù)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1) 唯一性準則：真實物理網(wǎng)格對象的編號和并行進程編號在所有進程中具有唯一性。

2) 同一性準則：虛擬網(wǎng)格對象的編號和并行進程編號與存在于其他進程上的真實網(wǎng)格對象的編號和并行進程編號相同。

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

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

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

(1)

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

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

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

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

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

Depair〈Ip,Sc〉

(2)

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

(3)

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

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

2) 面心：若為三角形面，C1、C2和C3分別取為該三角形3個頂點的編號，C4=Big；若為四邊形面，C1、C2、C3和C4分別取為該四邊形4個頂點的編號。

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

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

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

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

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

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

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

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

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

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

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

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

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

2 流場求解方法簡介

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

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

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

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

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

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

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

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

圖7 自適應前后流場對比Fig.7 Comparison of flow fields before and after adaptation

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

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

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

(4)

式中：T為一定迭代步數(shù)計算的平均時間；n為并行核數(shù)；下標0表示相應的基準值。

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

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

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

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

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

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

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

圖11 自適應3個模塊和并行通信耗時比例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é) 論

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

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

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

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

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