面向復雜幾何模型的多級并行四面體網格生成算法*

徐 權，劉田田，冷玨琳，楊 洋，鄭 澎

(1. 中物院高性能數值模擬軟件中心， 北京 100088； 2. 北京應用物理與計算數學研究所， 北京 100094；3. 中國工程物理研究院計算機應用研究所， 四川 綿陽 621900)

近年來，高性能計算機發展得非常快，預計國內外將會有多臺E級超級計算機研制成功，這些超級計算機的研制為高性能數值模擬的發展提供了必需的硬件環境。隨著高性能計算機的發展，國內外的研究人員已開發出諸多適應于高性能計算機的并行數值模擬軟件，同時也研制出了數值模擬領域的并行編程框架，如JASMIN[1]、JAUMIN[2]、PHG[3]等，這些編程框架均可以支持并行數值模擬程序的快速開發，使得應用程序快速擴展至數萬到數十萬處理器核。目前，在計算流體力學、復雜電磁環境、反應堆、大型工程力學分析等領域進行高精度的數值模擬所需的網格規模已高達數十億量級，而網格生成作為數值模擬中的重要一環，其發展則相對滯后于高性能計算機和數值模擬解法器的發展，針對并行網格生成算法的研究相對較少，國內外尚沒有比較成熟的適用于高性能計算機體系結構的并行網格生成軟件，而傳統的網格生成方法或軟件已經不能滿足高性能數值模擬的需求，因此需要面向高性能數值模擬研究復雜幾何模型的超大規模非結構網格并行生成方法，并與領域編程框架和并行數值模擬求解器實現無縫對接。

前沿推進方法(Advancing Front Technique, AFT)和Delaunay三角化方法是目前非結構四面體網格生成最常用的兩種方法，國內外研究人員基于這兩種方法已經開發出一批比較成熟的串行四面體網格生成軟件，如TetGen[4]、Gmsh[5]、Netgen[6]等，并得到了廣泛應用。當前并行四面體網格生成算法的研究主要基于AFT和Delaunay三角化開展，如：L?hner等[7-8]發展了并行前沿推進方法，提出了一種新的可擴展的并行前沿推進方法，并對當前AFT方法的發展趨勢進行了展望；Chrisochoides等[9]則基于其團隊多年的研究，提出了多級并行的Delaunay三角化方法，該方法可以適應于當前的高性能計算機體系結構，可擴展性較好，但需要復雜的通信和調度機制，編程難度較大；Ivanov、Chen等[10-11]分別提出了基于區域分解的并行網格生成方法，該方法一般需要對面網格進行區域分解，但會引入人工界面，算法的可擴展性較低；Pebay、Wang等[12-13]提出了基于初始網格的并行網格加密方法，該方法在進程間無須進行通信，具有良好的可擴展性，但網格質量相對較差；文獻[14-16]則分別提出了基于共享內存的多線程并行四面體網格生成方法。

目前的超級計算機普遍采用分布式-共享內存結構，計算節點之間是分布式內存，計算節點內部是共享式內存，每個計算節點上分布多個中央處理器(Central Processing Unit, CPU)，每個CPU內分布多個核，如圖1所示。采用分布式算法和共享內存算法耦合的并行網格生成可以更好地適應當前的高性能計算機體系結構，同時可以有效提高網格生成算法的可擴展性。

圖1 主流超級計算機的體系結構Fig.1 Architecture of current major supercomputers

本文針對高性能數值模擬領域對數十億乃至上百億四面體網格生成的需求，提出了一種適用于復雜幾何模型的多級并行四面體網格生成框架，同時提出了面向幾何實體的區域分解方法和多實體間的并行面網格生成方法，在生成四面體網格時復用了基于三角形網格的區域分解和多線程四面體網格生成方法。該算法可以適應現有的高性能計算機體系結構，實現多實體復雜幾何模型的大規模四面體網格的并行生成。

1 多級并行四面體網格生成

1.1 并行網格生成框架

本文的多級并行四面體網格生成算法是基于中物院高性能數值模擬軟件中心研制的并行數值模擬前處理軟件SuperMesh[17]研發的。SuperMesh是一款面向大規模復雜數值模擬的并行前處理引擎，它以“可計算幾何模型”為核心，建立了模型的統一表達方式，屏蔽了底層各種幾何核心，發展了模型處理和網格生成兩大子系統，具備面向復雜幾何模型的結構網格、非結構網格、組合幾何的并行生成能力，支持應用軟件前處理界面的快速定制，支持與并行領域編程框架JASMIN、JAUMIN和數值模擬應用程序的無縫對接，從而支持大規模的數值模擬。

基于并行前處理引擎SuperMesh，提出了面向復雜幾何模型的適應高性能計算機體系結構的多級并行四面體網格生成框架，該框架的整體流程如圖2所示。首先是對計算機輔助設計(Computer Aided Design, CAD)模型進行模型處理，使模型轉換為可計算幾何模型，然后自動識別模型的幾何特征，建立網格的尺寸場，基于尺寸場和幾何實體間的鄰接關系對CAD模型的幾何實體進行分組，并將分組后的幾何實體分配到不同的計算節點中，各節點之間并行生成三角形面網格，接著在每個計算節點內對三角形面網格進行二次區域分解，分解為多個子面網格，并將分解后的子面網格分配到各進程中，最后每個進程內采用多線程并行四面體網格生成方法實現大規模非結構四面體網格的生成。下面將對算法流程中的關鍵算法進行介紹。

圖2 多級并行四面體網格生成流程Fig.2 Flow chart of multilevel parallel tetrahedral mesh generation

1.2 基于幾何實體的區域分解

在實際工程應用中，CAD模型一般均是裝配體，裝配體包含了CAD模型的建模過程信息，并由多個幾何實體組成。為了提高并行網格生成方法的可擴展性，本文針對含多個幾何實體的CAD模型提出了基于幾何實體的區域分解方法，該方法可以保證網格生成中幾何實體的完整性和四面體網格的質量。在對CAD模型進行分組時，需要考慮兩方面因素：①減少計算節點間的通信；②保證計算節點間的負載平衡。為了減少通信，需要減少不同區域間的交界面的個數，即將相鄰的幾何實體分配到同一個計算節點，為此在分組時考慮了幾何實體間的鄰接關系。為了保證網格生成的負載平衡，本方法首先建立了網格的尺寸場，然后基于尺寸場對每個幾何實體的網格規模進行了預估，并在構建無向圖時作為圖中節點的權重，從而達到良好的負載平衡。基于幾何實體的區域分解算法步驟如算法1所示，流程如圖3所示。

算法1 基于幾何實體的區域分解算法Alg.1 Algorithm of domain decomposition based on geometric entities

圖3 基于幾何實體的區域分解流程Fig.3 Flow chart of domain decomposition based on geometric entities

針對包含1 876個幾何實體的三峽大壩模型進行了測試，測試使用24個處理器核，分別對無網格規模預估和加入網格規模預估進行了測試，兩次測試生成的四面體網格數分別為3 374萬和3 348萬。圖4給出了加入網格規模預估前后的每個處理器核內的網格數，可以看出加入網格規模預估后，不同處理器核上網格分布更加均勻，相比未進行網格規模預估的情況，整體上網格的負載平衡得到大大改善。

圖4 加入網格規模預估前后網格數對比Fig.4 Comparison of number of cells before and after mesh size prediction

1.3 并行三角形面網格生成

目前的數值模擬求解方法大都需要協調地計算網格，為了保證協調四面體網格的生成，需要保證子區域交界面上網格的一致，為此在并行三角形網格生成過程中，需要對每個節點上的邊界進行處理，從而保證節點間網格的協調一致。在并行前處理引擎SuperMesh中，已經研制了面向復雜幾何模型的“線—面—體”的三步分離式串行四面體網格生成算法，本文對該算法進行了擴展，實現了節點間“點—線—面”三角形網格的并行生成，該算法的步驟如算法2所示，流程如圖5所示。在計算節點內實現“點—線—面”網格后，對于不同的子區域間共享的節點、線和面，通過并行通信模塊使得節點間交界面上網格一致，從而保證協調的三角形網格生成。

算法2 并行三角形網格生成算法Alg.2 Algorithm of parallel triangular mesh generation

圖5 并行三角形網格生成流程Fig.5 Flow chart of parallel triangular mesh generation

1.4 基于三角形網格的區域分解方法

單個幾何實體比較復雜時，其需要生成的網格規模會非常大，生成四面體網格時間比較長，此時需要對單個幾何實體的三角形面網格做進一步分解，分解為多個較小的子面網格，然后再并行生成四面體網格。本文采用文獻[11]提出的基于三角形面網格的區域分解算法，該算法對三角形面網格采用遞歸的方式依次分解為子區域，每次分解時會在待分解區域內插入一個交界面網格，將區域一分為二，當所有子區域的網格規模都小于給定的閾值時，遞歸分解過程結束，最終實現面網格的區域分解，并將每個子面網格分發到不同的處理器中。

1.5 多線程四面體網格生成

Delaunay三角化方法由于具備數學基礎好、網格生成時間快、網格質量好等優勢，得到了廣泛的應用。多線程并行Delaunay算法關鍵是減少網格生成過程中多個點同時插入引起的空腔干涉現象，目前處理的思路基本是先將點集進行分組，然后各線程再并行插點，當兩個點同時插入引起空腔干涉時，放棄其中一個線程的操作，直到網格生成結束。但上述方法無法保證負載平衡，且同步等待的時間較長，因此本文采用了文獻[16]提出的改進的基于無鎖原子操作的并行Delaunay算法，該算法中令每一個四面體單元都附一個原子變量及一個占位標志。在計算Delaunay空腔的過程中，線程將其空腔中所有單元及邊界外一層單元的原子變量置為1。若有多個線程同時試圖將某個變量置位，則使用無鎖原子操作的機制保證有且僅有一個線程成功置位占位變量，在更新四面體單元的相鄰關系過程中，此操作僅僅涉及修改局部區域的單元，可以有效減少由于負載不平衡引起的同步等待時間，算法的并行效率較高。

2 數值實驗與結果分析

2.1 并行性能與分析

在并行算法性能測試中，并行加速比是較為常用的并行性能測試指標。一般地，并行加速比的定義為

S=T1/TP

其中，T1為單進程網格生成時間，TP為P個進程的并行網格生成時間。

由于測試模型網格生成規模較大，無法在單進程下運行，因此對并行加速比的計算公式進行修正，修正后的并行加速比為

S=(TN/TP)×N

其中，TN為N個進程時網格生成時間。

測試環境為曙光集群系統，該系統共有424個計算節點，系統采用Intel Omni-path高速互聯網格、曙光Parastor300全局并行存儲系統，每個計算節點均配置2顆Intel Xeon Gold 6132處理器，每個處理器14核，96 GB DDR4 ECC 2400 MHz內存。

測試模型選取了具有代表性的三峽大壩模型對算法進行了數值驗證(如圖6所示)。三峽大壩模型包含左廠房壩段、泄洪壩段、右廠房壩段、非溢流壩段、升船機、船閘、連接壩段，共計1 876個幾何實體。壩頂長度2 309 m，最大壩高181 m，整個模型包含了豐富的幾何細節，幾何尺寸跨度大，模型比較復雜。

圖6 三峽大壩CAD模型圖Fig.6 CAD model of Three Gorges Dam

在進行并行測試時，由于集群使用限制，本算法測試主要測試了千核下的并行加速比。測試時每個節點使用28個處理器核，從228個處理器核開始進行測試，測試網格尺寸為0.5 m，生成9.36億四面體網格單元，并行網格生成時間和加速比分別如表1和圖7所示。

表1 三峽大壩模型并行四面體網格生成加速比Tab.1 Speed-up ratio of parallel tetrahedral mesh generation for Three Gorges Dam model

圖7 三峽大壩模型四面體網格并行生成加速比Fig.7 Speed-up ratio of parallel tetrahedral mesh generation for Three Gorges Dam model

從表1和圖7可以看出，本算法針對三峽大壩模型，可以獲得良好的加速比，具備在數千處理器核上并行生成數十億規模四面體網格單元的能力。

與文獻[10-11]等提出的并行四面體網格生成算法相比，本文的算法在448處理器核下，9.36億四面體網格單元生成時間在5 min以內，整體上網格生成時間相對較少，同時本文的算法可擴展到千核以上，相對文獻[10-11]提出的算法具有更好的并行可擴展性。

2.2 四面體網格質量分析

圖8 三峽大壩模型四面體網格質量統計結果Fig.8 Statistical results of tetrahedral mesh quality of Three Gorges Dam model

從統計結果中可以看出，并行生成的四面體單元的體積邊長比均集中在0.5～0.9的范圍，基本不存在扁平的單元(<0.2)，生成的網格質量較高。

2.3 實際應用測試

三峽工程計算規模龐大，千萬級單元的模擬只能達到米級，考慮基礎、細節構造的高分辨率模擬，需要數十億自由度級的計算能力，而現有的自主行業軟件和商業軟件難以滿足計算需求。通過本文提出的并行網格生成算法，結合文獻[18]提出的并行加密算法，在“天河二號”上，采用49 152處理器核，并行生成101.1億四面體網格單元，并采用結構靜力與振動大規模并行分析軟件PANDA-StaVib首次實現了三峽大壩精細模型的靜力計算，為三峽大壩工作性態評估提供了能力支撐。圖9給出了三峽大壩的網格結果。

圖9 三峽大壩模型四面體網格生成結果Fig.9 Tetrahedral mesh of Three Gorges Dam model

3 結論

面向高性能數值模擬對大規模非結構網格生成的需求，提出了一種適應復雜幾何模型的多級并行四面體網格生成框架，并針對裝配體模型提出了基于幾何實體的區域分解算法和并行面網格生成算法，能適用具有數千上萬幾何實體的復雜幾何模型。在四面體網格并行生成時，該框架可以復用面網格區域分解和多線程四面體網格生成算法，相比單純的分布式并行算法或共享內存式并行算法，具備良好的并行效率和可擴展性，更好地適應了當前的高性能計算機體系結構，可以在數千處理器核上并行生成數十億規模四面體網格單元。

本文提出的并行非結構網格生成框架，可以推廣到含有多個幾何實體的模型的其他類型網格生成算法，如六面體網格、混合網格等，僅需將框架中的三角形網格生成和四面體網格生成的單元算法替換為相應的網格生成算法，具有較好的推廣應用價值。