非結構混合網格自適應并行技術

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

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

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

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

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

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

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

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

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

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

1 網格自適應系統的并行技術

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

圖1 并行網格自適應數值模擬執行流程Fig.1 Flowchart of parallel numerical simulation with mesh adaptation

在基于網格分區的CFD并行實現中，虛擬網格技術常用來存儲并行流場數據，以實現物理邊界條件和并行邊界條件的統一處理[10]。本文虛擬網格單元的模板擴展為與網格點相連的網格單元，如圖2所示，實線圍成的網格單元為真實物理單元，虛線相關的網格單元為虛擬單元。

圖2 并行邊界虛擬網格單元Fig.2 Ghost elements on parallel interface

通常情況下，并行算法與自適應系統采用的具體方法有關。本文針對非結構混合網格，且主要基于文獻[9]中h型混合網格自適應實現方法，分別設計了自適應網格優化3項技術的并行算法。另外，針對網格加密或稀疏造成的負載不平衡，本文在兩步法的基礎上設計了適合多級網格自適應的動態負載平衡算法。

1.1 網格單元分布優化并行技術

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

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

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

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

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

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

(1)

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

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

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

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

對于新增網格單元的特征數據De，取為父級單元編號Ip和子單元序號Sc組成的數對：

Depair〈Ip,Sc〉

(2)

對于新增網格點的特征數據Dp，取如下四元素：

(3)

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

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

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

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

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

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

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

1.2 表面網格幾何投影并行技術

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

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

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

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

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

1.3 空間網格匹配并行技術

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

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

2 流場求解方法簡介

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

3 動態負載平衡技術

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

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

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

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

圖6 網格單元數據打包/解包示意圖Fig.6 Schematic diagram of data packing/unpacking for mesh elements

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

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

圖8 自適應負載平衡并行網格動態變化過程Fig.8 Dynamical variation of parallel mesh during mesh adaptation for load balance

4 網格自適應并行技術效率測試

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

(4)

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

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

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

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

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

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

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

表1 自適應系統及其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 密度殘差和升力系數收斂曲線Fig.12 Convergence history of density residual and lift coefficient

5 結 論

針對非結構混合網格，設計了網格單元分布優化、表面網格幾何投影和空間網格匹配相關的并行算法。同時，為后續流場模擬建立了基于并行重分區-網格數據遷移的動態負載平衡技術。最后，采用三角翼模型測試了不同并行規模下3項并行算法的并行性能。

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

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

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

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