大規(guī)模四面體網(wǎng)格并行生成方法

王小慶，金先龍

(1.上海交通大學(xué) 機(jī)械系統(tǒng)與振動國家重點(diǎn)實(shí)驗(yàn)室，上海 200240;2.上海交通大學(xué) 機(jī)械與動力工程學(xué)院，上海 200240)

網(wǎng)格生成是整個有限元計算的重要一環(huán)。在大型工程領(lǐng)域，比如抗震分析，抗爆炸沖擊分析等，數(shù)值計算方法的應(yīng)用越來越多，對求解精度的要求也越來越高，而這些對單元網(wǎng)格的需求規(guī)模也越來越大［1－4］，傳統(tǒng)的串行網(wǎng)格生成限于時間及內(nèi)存的限制［5－6］，很難滿足超大規(guī)模數(shù)值計算的精度需要，而超級計算機(jī)的快速發(fā)展，也為網(wǎng)格生成的并行化提供了必需的硬件環(huán)境，并行網(wǎng)格劃分已經(jīng)成為網(wǎng)格生成領(lǐng)域的研究熱點(diǎn)。

并行網(wǎng)格生成的研究方法主要集中在算法并行和問題并行兩個方面。采用問題并行方式的研究較多，Laemmer等［7－8］提出了直接對幾何體進(jìn)行區(qū)域分解的二維三角形和三維四面體網(wǎng)格的并行生成方法，其對幾何體的分區(qū)分別采用最大和最小極慣性軸遞歸對分的方法，提供了幾何體劃分的一種新思路，但是較難達(dá)到負(fù)載的均衡。Ivanov等［9－10］提出了對初始面網(wǎng)格分區(qū)，然后再并行生成體網(wǎng)格的方法。Ito等［11］提出了對初始體網(wǎng)格進(jìn)行分區(qū)，然后利用前沿推進(jìn)法進(jìn)行并行網(wǎng)格生成的方法。在算法并行方面，Blelloch等［12－13］闡述了并行 Delaunay網(wǎng)格生成方法。Lhner等［14］對前沿推進(jìn)法的并行化進(jìn)行了研究。國內(nèi)方面，關(guān)振群等［15］總結(jié)了有限元網(wǎng)格生成方法的研究進(jìn)展，指出網(wǎng)格并行生成是未來有限元網(wǎng)格的發(fā)展方向。李水鄉(xiāng)等［16］采用問題并行方式，實(shí)現(xiàn)了基于局域網(wǎng)的小規(guī)模有限元網(wǎng)格的并行生成。齊龍等［17］對GPU在并行網(wǎng)格生成中的應(yīng)用進(jìn)行了研究。王磊等［18］總結(jié)了各類Delaunay四面體網(wǎng)格并行生成算法的優(yōu)缺點(diǎn)，并探討了其發(fā)展趨勢。

本文在上述研究基礎(chǔ)上，提出了一種大規(guī)模非結(jié)構(gòu)化四面體網(wǎng)格并行生成方法，采用問題并行的方式，直接讀取CAD幾何體文件，然后進(jìn)行初始網(wǎng)格劃分，并通過相對體積比及最優(yōu)分區(qū)來控制初始網(wǎng)格數(shù)量，接著采用圖劃分理論對有限元網(wǎng)格進(jìn)行區(qū)域分解，然后各個分區(qū)采用分裂法進(jìn)行并行網(wǎng)格劃分。利用初始體網(wǎng)格進(jìn)行分區(qū)克服了利用初始面網(wǎng)格進(jìn)行區(qū)域分解所帶來的負(fù)載不均衡問題。提出的初始網(wǎng)格數(shù)量控制方法，最大限度的降低了此串行部分對并行效率的影響。針對分區(qū)間節(jié)點(diǎn)匹配的難題，提出了一種基于共享單元的邊界處理方法，有效地解決了分區(qū)間邊界協(xié)調(diào)的問題。最后通過算例證明了整個四面體并行生成方法可以獲得良好的加速比和并行效率，所產(chǎn)生的網(wǎng)格具有較高質(zhì)量可以直接用于后續(xù)并行計算。

圖1 并行網(wǎng)格生成流程圖Fig.1 Flow Chart of parallel mesh generation

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

此處所闡述的并行網(wǎng)格劃分方法，作為后續(xù)自主開發(fā)的并行求解器的前處理器的核心算法，克服了傳統(tǒng)串行網(wǎng)格劃分時間及內(nèi)存耗費(fèi)的瓶頸［19］，程序采用C++和Fortran語言編寫，利用MVAPICH并行庫實(shí)現(xiàn)并行通信。本文采用問題并行方式，串行網(wǎng)格劃分部分采用前沿推進(jìn)法。

1.1 并行網(wǎng)格劃分算法

如圖1為并行網(wǎng)格劃分的流程圖。輸入文件為通用CAD幾何文件，采用前沿推進(jìn)法［20］進(jìn)行初始網(wǎng)格劃分，并控制其產(chǎn)生合理的初始網(wǎng)格數(shù);接著利用圖論理論的圖劃分方法，進(jìn)行區(qū)域分解，將各區(qū)域的數(shù)據(jù)發(fā)送到其余進(jìn)程，然后采用分裂法進(jìn)行并行網(wǎng)格生成。后續(xù)通過邊界條件處理提交并行求解器進(jìn)行求解計算，最后形成可視化后處理文件。在整個文件系統(tǒng)中，除了初始幾何文件為單一文件傳遞外，其余文件均采用多文件傳輸處理，每個核負(fù)責(zé)處理自己的文件，并行多文件的利用，在很大程度上提高了程序的執(zhí)行效率。

1.2 最優(yōu)初始網(wǎng)格生成方法

初始網(wǎng)格劃分的數(shù)量決定整個算法中串行部分的時間，劃分?jǐn)?shù)量過多，串行時間會過長，因而降低并行效率，劃分?jǐn)?shù)量過少，會導(dǎo)致模型失真，網(wǎng)格質(zhì)量下降，同時所能進(jìn)行的分區(qū)數(shù)也受到限制，影響并行規(guī)模的擴(kuò)大。為合理解決此問題，此處采用雙重約束對初始網(wǎng)格數(shù)量進(jìn)行控制和優(yōu)化。首先，保證初始網(wǎng)格的總體積和原幾何體體積相對差值不大于5%，從而保證初始網(wǎng)格的幾何貼切度，使有限元模型不失真，如式(1)為其約束條件。不同相對體積比初始網(wǎng)格效果如圖2所示，可以看到對于單位球，其初始網(wǎng)格的體積和原幾何體體積相對差值控制在4.8%時基本能保證模型不失真。其次，在保證最佳分區(qū)效果即保證分區(qū)均衡度的情況下，確定每個分區(qū)最少單元數(shù)，由于圖剖分和圖本身的拓?fù)溆泻艽笙嚓P(guān)性，因此只能通過多次試驗(yàn)的方法確定其在一定分區(qū)數(shù)時對應(yīng)的最佳單元數(shù)。

式中 a，b，c是共頂點(diǎn)的三條棱邊的長度，α，β，γ 為相鄰棱邊組成的面角，ω為這三個面角和的一半。

圖2 不同體積比初始網(wǎng)格效果示意圖Fig.2 Effect of different volume ration of initial mesh

1.3 區(qū)域分解技術(shù)

有限單元網(wǎng)格的區(qū)域分解，采用圖論圖劃分方法。首先將單元信息轉(zhuǎn)換為其對應(yīng)的對偶圖，然后將對偶圖利用圖論的方法進(jìn)行區(qū)域分解，其過程分為粗化圖階段，圖分區(qū)階段及圖還原階段，最后再將分區(qū)結(jié)果映射還原為有限元網(wǎng)格的區(qū)域分解信息。圖3所示為有限元網(wǎng)格與其相應(yīng)對偶圖的相互對應(yīng)關(guān)系，每個有限元單元對應(yīng)對偶圖中的一個頂點(diǎn)，具有共享邊(面)的單元，在對偶圖中其相應(yīng)的頂點(diǎn)之間，形成一條邊。圖中頂點(diǎn)和邊上的數(shù)字代表權(quán)重，此處所有節(jié)點(diǎn)和邊按照等權(quán)重處理，大小均置為1。

圖3 有限元網(wǎng)格及其對偶圖Fig.3 Finite element model and its dual graph

有限元模型轉(zhuǎn)換形成的對偶圖，利用CSR(compressed storage)格式進(jìn)行存儲。CSR格式在存儲稀疏圖領(lǐng)域有廣泛應(yīng)用。這種格式利用兩個數(shù)組來表達(dá)圖的鄰接信息，第一個數(shù)組存儲第i個頂點(diǎn)相鄰頂點(diǎn)的開始和終了地址索引，第二個數(shù)組存儲著每個頂點(diǎn)的鄰接頂點(diǎn)，根據(jù)第一個數(shù)組中第i個頂點(diǎn)的索引可以在第二個數(shù)組中找到與其鄰接的頂點(diǎn)。程序中具體實(shí)現(xiàn)采用metis的k路圖劃分方法［21］。基于圖論理論的區(qū)域分解能夠保證負(fù)載的均衡，提高并行效率。

1.4 基于分裂法的并行網(wǎng)格生成

初始網(wǎng)格進(jìn)行區(qū)域分解后，由主進(jìn)程發(fā)往各個進(jìn)程，各從進(jìn)程接收到數(shù)據(jù)后，利用分裂法對網(wǎng)格進(jìn)行并行加密，實(shí)現(xiàn)網(wǎng)格的并行生成，直到滿足網(wǎng)格數(shù)量及精度要求［20，22］。分裂法是在原單元各個邊的中點(diǎn)生成新的節(jié)點(diǎn)，作為新增網(wǎng)格的頂點(diǎn)，然后分裂產(chǎn)生新的單元，一個三角形單元分裂為四個三角形單元，一個四面體單元分裂為8個四面體單元，其過程如圖4所示。

圖4 二維、三維單元分裂示意圖Fig.4 Schematic of 2D and 3D element split

采用分裂法很大程度上減少了各個分區(qū)之間的通信量，但是為了保證網(wǎng)格的整體性，各個分區(qū)邊界上的節(jié)點(diǎn)應(yīng)該保持一致的編號，否則無法進(jìn)行后續(xù)計算或者網(wǎng)格的合并。為此，各個分區(qū)邊界之間需要保持一定的通信，這里提出了一種基于共享單元的邊界判定方法，首先根據(jù)對偶圖區(qū)域分解結(jié)果，確定邊界四面體單元;通過遍歷所有邊界單元的邊，并在其所有共享進(jìn)程之間進(jìn)行通信，確定邊界共享邊;采用隨機(jī)算法在共享進(jìn)程間分配共享邊;劃分每個分區(qū)的邊為內(nèi)部邊、屬于自己的邊界邊、不屬于自己的邊界邊三種類型，分別插入新的節(jié)點(diǎn)，并在相應(yīng)邊界之間交換節(jié)點(diǎn)信息;產(chǎn)生新的單元;進(jìn)入下一個循環(huán)。基于共享單元的方法，借鑒了文獻(xiàn)［10］基于共享節(jié)點(diǎn)的方法，同時克服了其無法處理圖5情況的一個缺陷，圖示該邊的兩個節(jié)點(diǎn)96和113被兩個進(jìn)程共享，采用基于共享節(jié)點(diǎn)的方法會判定其為邊界邊，但是事實(shí)上它是內(nèi)部邊;而采用基于單元的邊界判定方法，通過對該邊界單元四條邊的遍歷，及與其共享進(jìn)程的通信，發(fā)現(xiàn)其共享進(jìn)程Core1并沒有該邊，這樣可以判定該邊為內(nèi)部邊，而非邊界邊。基于共享單元分裂法的具體實(shí)現(xiàn)流程如圖6所示。

圖5 特殊邊界棱邊示意圖Fig.5 Schematic of special edge

2 算例分析

本節(jié)通過算例分析該算法的有效性和實(shí)用性。算例測試環(huán)境均為上海超級計算機(jī)中心“魔方”曙光5000A超級計算機(jī)，整機(jī)理論計算峰值為200 Tflops，具體計算區(qū)域?yàn)镃區(qū)，該區(qū)由800個刀片式服務(wù)器組成，每個服務(wù)器包含四顆完全相同的AMD Barcelona 1.9 GHz四核處理器，每個刀片共享內(nèi)存為64 G。所有服務(wù)器操作系統(tǒng)均為64位SuSE Linux Enterprise Server 10 SP2，通過LSF作業(yè)系統(tǒng)進(jìn)行作業(yè)管理。

2.1 算例1

算例1選擇在電力、冶金、石油、化工、水利等行業(yè)有廣泛用途的滑動式馬鞍支座，通過該算例驗(yàn)證算法的并行性能及可擴(kuò)展性。首先通過圖7直觀的展現(xiàn)采用八個計算核心時，并行網(wǎng)格劃分的具體實(shí)現(xiàn)過程，即包括幾何體讀入、初始網(wǎng)格劃分、區(qū)域分解以及并行網(wǎng)格生成等。

圖6 基于分裂法的并行網(wǎng)格生成流程圖Fig.6 Flow chart of parallel mesh generation based on the split method

通過不同規(guī)模的網(wǎng)格及不同規(guī)模的核數(shù)，對算例進(jìn)行分析，評估其時間性能，基于初始網(wǎng)格數(shù)量的控制方法，對于128核及其以下核數(shù)時，串行部分初始網(wǎng)格取877個，256核及其以上核數(shù)時，取初始網(wǎng)格取7 016個。如表1所示為不同規(guī)模不同核數(shù)時，網(wǎng)格劃總分時間。可以看到，對于千萬數(shù)量網(wǎng)格，可以在幾秒內(nèi)完成網(wǎng)格的劃分;而對串行很難完成的億以及十億量級的網(wǎng)格劃分，利用該并行網(wǎng)格生成方法可以很快完成，其中對將近二十億的網(wǎng)格采用1024個處理核時，可以在大約0.28 min內(nèi)完成網(wǎng)格劃分。

圖7 八分區(qū)并行網(wǎng)格生成示意圖Fig.7 Diagram of parallel mesh generation with 8 partitions

表1 生成不同數(shù)量網(wǎng)格總時間(s)Tab.1 Time of different amounts of element generation

并行算法的可擴(kuò)展性是指隨問題規(guī)模及計算資源的同比擴(kuò)大，并行算法效率的穩(wěn)定性。加速比是衡量可擴(kuò)展性的常用方法。如圖8為在問題規(guī)模固定為億量級網(wǎng)格時，加速比隨核數(shù)的變化。由圖可以看到，在不考慮串行部分即初始網(wǎng)格劃分和分區(qū)時間的情況下，該算法可以獲得相當(dāng)好的加速比。在考慮串行部分后，整體的加速比有所下降，128核、256核和512核均可以達(dá)到70%的并行效率，在1 024核情況下仍然達(dá)到大約56%的并行效率。

2.2 算例2

算例2選擇具有實(shí)際工程應(yīng)用背景的某核電站防浪堤抗震分析模型進(jìn)行大規(guī)模非結(jié)構(gòu)化網(wǎng)格的并行生成，主要考察生成網(wǎng)格的質(zhì)量。整個模型長度2 km、寬度950 m、高度210 m。該算例采用1024個核進(jìn)行網(wǎng)格的并行劃分，總運(yùn)行時間為100.3 s，最終的網(wǎng)格單元數(shù)2 750 414 848，節(jié)點(diǎn)數(shù)480 617 472，如圖9所示為防浪堤抗震分析模型1024分區(qū)整體網(wǎng)格生成結(jié)果及局部網(wǎng)格詳圖。

圖8 并行網(wǎng)格生成加速比Fig.8 Speedup of parallel mesh generation

圖9 防浪堤抗震模型1024分區(qū)并行網(wǎng)格劃分結(jié)果Fig.9 Parallel mesh generation result of anti-seismic breakwater model with 1024 partitions

并行網(wǎng)格劃分方法克服了串行劃分方法的時間及內(nèi)存的瓶頸，但是其網(wǎng)格質(zhì)量必須滿足后續(xù)串行或并行求解器的需要，為此，針對該算例分析其網(wǎng)格質(zhì)量，這里將網(wǎng)格質(zhì)量Q定義如式:

式中，VolIdeal=8r3/9，為與該四面體同一個外接球的理想四面體的體積，r為該四面體的外接球半徑。VolActual為每個四面體的體積，通過式(2)可求得。理想的Q值為0，即該四面體為正四面體時，網(wǎng)格質(zhì)量最高。如圖10為防浪堤抗震模型網(wǎng)格質(zhì)量的統(tǒng)計結(jié)果圖，可以看到，Q值集中在0－0.3之間，和商業(yè)軟件串行產(chǎn)生的網(wǎng)格質(zhì)量分布情況類似，可以滿足后續(xù)計算的需要。

圖10 防浪堤抗震模型網(wǎng)格質(zhì)量分布圖Fig.10 Mesh quality distribution of anti-seismic breakwater model

3 結(jié)論

(1)所提出的非結(jié)構(gòu)化四面體網(wǎng)格并行生成方法，克服了串行網(wǎng)格劃分存在的瓶頸問題，通過算例證明了其高效性和所生成網(wǎng)格的實(shí)用性。

(2)提出的初始網(wǎng)格數(shù)量控制方法，減少了串行部分影響，有效地提高了并行效率;所述的基于共享單元的邊界處理方式，有效地解決了各個分區(qū)之間邊界節(jié)點(diǎn)的匹配，對類似邊界之間的協(xié)調(diào)問題具有借鑒意義。

［1］鄧榮兵，金先龍，陳峻.爆炸流場與玻璃幕墻動力響應(yīng)的仿真計算方法［J］.振動與沖擊，2011，30(3):14－17.DENG Rong-bing，JIN Xian-long，CHEN Jun.Numerical simulation method for blast flow and dynamic of glass curtain wall［J］.Journal of Vibration and Shock，2011，30(3):14 －17.

［2］ Zhou M，Xie T，Seol S，et al.Tools to support mesh adaptation on massively parallel computers［J］.Engineering with Computers，2012，28(3):287－301.

［3］Jou W.Comments on the feasibility of LES for commercial airplane wings［R］.AIAA －98－2801，1998.

［4］ Baum J D，Luo H，Lhner R.Numerical simulation of blast in the world trade center［C］.AIAA，Aerospace Sciences Meeting and Exhibit，33rd，Reno，NV.1995.

［5］馬新武，王芳，趙國群.具有內(nèi)部特征約束的四邊形網(wǎng)格生成方法［J］.計算力學(xué)學(xué)報，2012，29(6):855－900.MA Xin-wu，WANGFang，ZHAO Guo-qun.An algorithm for quadrilateral mesh generation with internal feature constraints［J］.Chinese Journal of Computational Mechanics，2012，29(6):855－900.

［6］田素壘，張志毅，陳敏，等.四面體網(wǎng)格生成方法的研究與實(shí)現(xiàn)［J］.計算機(jī)工程與設(shè)計，2012，33(11):4416－4421.TIAN Su-lei，ZHANG Zhi-yi，CHEN Min，et al.Research and implementation of methods for tetrahedral mesh generation［J］.Computer Engineering and Design，2012，33(11):4416－4421.

［7］Laemmer I L.Parallel mesh generation［M］.Solving Irregularly Structured Problems in Parallel.Springer Berlin Heidelberg，1997:1 －12.

［8］ AndrH，Gluchshenko O N，Ivanov E G，et al.Automatic parallel generation of tetrahedral grids by using a domain decomposition approach［J］.Computational Mathematics and Mathematical Physics，2008，48(8):1367 －1375.

［9］ Ivanov E G， Andr H， Kudryavtsev A N. Domain decomposition approach for automatic parallel generation of tetrahedral grids［J］.Comput.Methods Appl.Math.，2006，6(2):178－193.

［10］ Ylmaz Y，zturan C，Tosun O，et al. Parallel mesh generation、migration and partitioning for the elmer application［R］.PRACE，2010:1－12.

［11］ Ito Y，Shih A M，Erukala A K，et al.Parallel unstructured mesh generation by an advancing front method［J］.Mathematics and Computers in Simulation，2007，75(5 －6):200－209.

［12］ Blelloch G E，Miller G L，Hardwick J C，et al.Design and implementation of a practical parallel Delaunay algorithm［J］.Algorithmica，1999，24(3－4):243－269.

［13］Chrisochoides N，Nave D.Parallel Delaunay mesh generation kernel［J］.International Journal for Numerical Methods in Engineering，2003，58(2):161－176.

［14］ Lhner R.A parallel advancing front grid generation scheme［J］. International Journal for Numerical Methods in Engineering，2001，51(6):663－678.

［15］關(guān)振群，宋超，顧元憲，等.有限元網(wǎng)格生成方法研究的新進(jìn)展［J］.計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，2003，15(1):1－14.GUAN Zhen-qun， SONG Chao， GU Yuan-xian， et al.Recent advances of research on finite element mesh generation methods［J］.Journal of Computer-Aided Design ＆ Computer Graphics，2003，15(1):1 －14.

［16］李水鄉(xiāng)，王云鵬，陳永強(qiáng).基于局域網(wǎng)的有限元網(wǎng)格分布式并行生成［J］.計算機(jī)工程與設(shè)計，2005，26(12):3165－3166.LI Shui-xiang， WANG Yun-peng， CHEN Yong-qiang.Distributed parallel FEM mesh generation on LAN［J］.Computer Engineering and Design，2005，26(12):3165－3166.

［17］齊龍，肖素梅，劉云楚，等.基于GPU的并行非結(jié)構(gòu)網(wǎng)格生成技術(shù)研究［J］.機(jī)械設(shè)計與制造，2013(2):184－186.QI Long，XIAO Su-mei，LIU Yun-chu，et al.Research on parallel unstructured mesh generation technology based on GPU［J］.Machinery Design ＆ Manufacturer，2013(2):184－186.

［18］王磊，聶玉峰，李義強(qiáng).Delaunay四面體網(wǎng)格并行生成算法研究進(jìn)展［J］.計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，2011，23(6):923－932.WANGLei，NIE Yu-feng，LI Yi-qiang，Advances of research on parallel Delaunay tetrahedral mesh generation ［J］.Journal of Computer-Aided Design ＆ Computer Graphics，2011，23(6):923 －932.

［19］ Larwood B G，Weatherill N P，Hassan O，et al.Domain decomposition approach for parallel unstructured meshgeneration［J］. International journal for numerical methods in engineering，2003，58(2):177－188.

［20］ Schberl J. NETGEN:an advancing front 2D/3D-mesh generator based on abstract rules［J］. Computing and visualization in science，1997，1(1):41－52.

［21］ Karypis G， Kumar V. Metis:a software package for partitioning unstructured graphs，partitioning meshes，and computing fill-reducing orderings of sparse matrices，version 5.1.0［M］.Minneapolis:Department of Computer Science and Engineering University of Minnesota，2013:26－30.

［22］Garatani K，Nakamura H，Okuda H，et al.GeoFEM:High performance parallel FEM for solid earth［C］.High-Performance Computing and Networking.Springer Berlin Heidelberg，1999:133－140.