徐 權(quán)，劉田田，冷玨琳，楊 洋，鄭 澎

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

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

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

目前的超級(jí)計(jì)算機(jī)普遍采用分布式-共享內(nèi)存結(jié)構(gòu)，計(jì)算節(jié)點(diǎn)之間是分布式內(nèi)存，計(jì)算節(jié)點(diǎn)內(nèi)部是共享式內(nèi)存，每個(gè)計(jì)算節(jié)點(diǎn)上分布多個(gè)中央處理器(Central Processing Unit, CPU)，每個(gè)CPU內(nèi)分布多個(gè)核，如圖1所示。采用分布式算法和共享內(nèi)存算法耦合的并行網(wǎng)格生成可以更好地適應(yīng)當(dāng)前的高性能計(jì)算機(jī)體系結(jié)構(gòu)，同時(shí)可以有效提高網(wǎng)格生成算法的可擴(kuò)展性。

圖1 主流超級(jí)計(jì)算機(jī)的體系結(jié)構(gòu)Fig.1 Architecture of current major supercomputers

本文針對(duì)高性能數(shù)值模擬領(lǐng)域?qū)?shù)十億乃至上百億四面體網(wǎng)格生成的需求，提出了一種適用于復(fù)雜幾何模型的多級(jí)并行四面體網(wǎng)格生成框架，同時(shí)提出了面向幾何實(shí)體的區(qū)域分解方法和多實(shí)體間的并行面網(wǎng)格生成方法，在生成四面體網(wǎng)格時(shí)復(fù)用了基于三角形網(wǎng)格的區(qū)域分解和多線程四面體網(wǎng)格生成方法。該算法可以適應(yīng)現(xiàn)有的高性能計(jì)算機(jī)體系結(jié)構(gòu)，實(shí)現(xiàn)多實(shí)體復(fù)雜幾何模型的大規(guī)模四面體網(wǎng)格的并行生成。

1 多級(jí)并行四面體網(wǎng)格生成

1.1 并行網(wǎng)格生成框架

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

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

圖2 多級(jí)并行四面體網(wǎng)格生成流程Fig.2 Flow chart of multilevel parallel tetrahedral mesh generation

1.2 基于幾何實(shí)體的區(qū)域分解

在實(shí)際工程應(yīng)用中，CAD模型一般均是裝配體，裝配體包含了CAD模型的建模過程信息，并由多個(gè)幾何實(shí)體組成。為了提高并行網(wǎng)格生成方法的可擴(kuò)展性，本文針對(duì)含多個(gè)幾何實(shí)體的CAD模型提出了基于幾何實(shí)體的區(qū)域分解方法，該方法可以保證網(wǎng)格生成中幾何實(shí)體的完整性和四面體網(wǎng)格的質(zhì)量。在對(duì)CAD模型進(jìn)行分組時(shí)，需要考慮兩方面因素：①減少計(jì)算節(jié)點(diǎn)間的通信；②保證計(jì)算節(jié)點(diǎn)間的負(fù)載平衡。為了減少通信，需要減少不同區(qū)域間的交界面的個(gè)數(shù)，即將相鄰的幾何實(shí)體分配到同一個(gè)計(jì)算節(jié)點(diǎn)，為此在分組時(shí)考慮了幾何實(shí)體間的鄰接關(guān)系。為了保證網(wǎng)格生成的負(fù)載平衡，本方法首先建立了網(wǎng)格的尺寸場(chǎng)，然后基于尺寸場(chǎng)對(duì)每個(gè)幾何實(shí)體的網(wǎng)格規(guī)模進(jìn)行了預(yù)估，并在構(gòu)建無向圖時(shí)作為圖中節(jié)點(diǎn)的權(quán)重，從而達(dá)到良好的負(fù)載平衡?；趲缀螌?shí)體的區(qū)域分解算法步驟如算法1所示，流程如圖3所示。

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

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

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

圖4 加入網(wǎng)格規(guī)模預(yù)估前后網(wǎng)格數(shù)對(duì)比Fig.4 Comparison of number of cells before and after mesh size prediction

1.3 并行三角形面網(wǎng)格生成

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

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

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

1.4 基于三角形網(wǎng)格的區(qū)域分解方法

單個(gè)幾何實(shí)體比較復(fù)雜時(shí)，其需要生成的網(wǎng)格規(guī)模會(huì)非常大，生成四面體網(wǎng)格時(shí)間比較長(zhǎng)，此時(shí)需要對(duì)單個(gè)幾何實(shí)體的三角形面網(wǎng)格做進(jìn)一步分解，分解為多個(gè)較小的子面網(wǎng)格，然后再并行生成四面體網(wǎng)格。本文采用文獻(xiàn)[11]提出的基于三角形面網(wǎng)格的區(qū)域分解算法，該算法對(duì)三角形面網(wǎng)格采用遞歸的方式依次分解為子區(qū)域，每次分解時(shí)會(huì)在待分解區(qū)域內(nèi)插入一個(gè)交界面網(wǎng)格，將區(qū)域一分為二，當(dāng)所有子區(qū)域的網(wǎng)格規(guī)模都小于給定的閾值時(shí)，遞歸分解過程結(jié)束，最終實(shí)現(xiàn)面網(wǎng)格的區(qū)域分解，并將每個(gè)子面網(wǎng)格分發(fā)到不同的處理器中。

1.5 多線程四面體網(wǎng)格生成

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

2 數(shù)值實(shí)驗(yàn)與結(jié)果分析

2.1 并行性能與分析

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

S=T1/TP

其中，T1為單進(jìn)程網(wǎng)格生成時(shí)間，TP為P個(gè)進(jìn)程的并行網(wǎng)格生成時(shí)間。

由于測(cè)試模型網(wǎng)格生成規(guī)模較大，無法在單進(jìn)程下運(yùn)行，因此對(duì)并行加速比的計(jì)算公式進(jìn)行修正，修正后的并行加速比為

S=(TN/TP)×N

其中，TN為N個(gè)進(jìn)程時(shí)網(wǎng)格生成時(shí)間。

測(cè)試環(huán)境為曙光集群系統(tǒng)，該系統(tǒng)共有424個(gè)計(jì)算節(jié)點(diǎn)，系統(tǒng)采用Intel Omni-path高速互聯(lián)網(wǎng)格、曙光Parastor300全局并行存儲(chǔ)系統(tǒng)，每個(gè)計(jì)算節(jié)點(diǎn)均配置2顆Intel Xeon Gold 6132處理器，每個(gè)處理器14核，96 GB DDR4 ECC 2400 MHz內(nèi)存。

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

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

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

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

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

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

與文獻(xiàn)[10-11]等提出的并行四面體網(wǎng)格生成算法相比，本文的算法在448處理器核下，9.36億四面體網(wǎng)格單元生成時(shí)間在5 min以內(nèi)，整體上網(wǎng)格生成時(shí)間相對(duì)較少，同時(shí)本文的算法可擴(kuò)展到千核以上，相對(duì)文獻(xiàn)[10-11]提出的算法具有更好的并行可擴(kuò)展性。

2.2 四面體網(wǎng)格質(zhì)量分析

圖8 三峽大壩模型四面體網(wǎng)格質(zhì)量統(tǒng)計(jì)結(jié)果Fig.8 Statistical results of tetrahedral mesh quality of Three Gorges Dam model

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

2.3 實(shí)際應(yīng)用測(cè)試

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

圖9 三峽大壩模型四面體網(wǎng)格生成結(jié)果Fig.9 Tetrahedral mesh of Three Gorges Dam model

3 結(jié)論

面向高性能數(shù)值模擬對(duì)大規(guī)模非結(jié)構(gòu)網(wǎng)格生成的需求，提出了一種適應(yīng)復(fù)雜幾何模型的多級(jí)并行四面體網(wǎng)格生成框架，并針對(duì)裝配體模型提出了基于幾何實(shí)體的區(qū)域分解算法和并行面網(wǎng)格生成算法，能適用具有數(shù)千上萬幾何實(shí)體的復(fù)雜幾何模型。在四面體網(wǎng)格并行生成時(shí)，該框架可以復(fù)用面網(wǎng)格區(qū)域分解和多線程四面體網(wǎng)格生成算法，相比單純的分布式并行算法或共享內(nèi)存式并行算法，具備良好的并行效率和可擴(kuò)展性，更好地適應(yīng)了當(dāng)前的高性能計(jì)算機(jī)體系結(jié)構(gòu)，可以在數(shù)千處理器核上并行生成數(shù)十億規(guī)模四面體網(wǎng)格單元。

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