面向JAUMIN的并行AFT四面體網(wǎng)格生成*

鄭　澎，方　維，徐　權(quán)，冷玨琳，熊　敏，于長華

1.中國工程物理研究院 計算機應(yīng)用研究所，四川 綿陽 621900

2.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100088

3.中國工程物理研究院 高性能數(shù)值模擬軟件中心，北京 100088

1　引言

并行自適應(yīng)非結(jié)構(gòu)網(wǎng)格支撐軟件框架JAUMIN[1]（J adaptive unstructured mesh applications infrastructure）是北京應(yīng)用物理與計算數(shù)學(xué)研究所研制的面向拉格朗日、歐拉等有限元計算方法的應(yīng)用軟件編程框架，它提升了非結(jié)構(gòu)網(wǎng)格并行應(yīng)用軟件的計算效率，降低了研制難度，并已在重大科學(xué)裝置結(jié)構(gòu)力學(xué)分析與優(yōu)化設(shè)計、裂變能源等重大項目中得到應(yīng)用。在這些應(yīng)用中，網(wǎng)格生成是實現(xiàn)高保真數(shù)值模擬的重要步驟，它需要生成能精確刻畫復(fù)雜幾何區(qū)域的高質(zhì)量和高精度網(wǎng)格，同時還需要網(wǎng)格生成、網(wǎng)格數(shù)據(jù)結(jié)構(gòu)必須與求解程序的并行操作相容，因此并行網(wǎng)格生成是本文研究的重點。

在有限元計算中，四面體網(wǎng)格能很好地逼近模型邊界且自動化生成的程度高，比較有代表性的串行四面體網(wǎng)格生成方法有AFT[2-3]（advancing front technique）方法和Delaunay方法[4-5]。但對于同一幾何模型，非結(jié)構(gòu)六面體網(wǎng)格雖能用更小規(guī)模的網(wǎng)格取得同樣的計算精度，生成過程需要更少的內(nèi)存，但針對任意復(fù)雜幾何模型自動生成六面體網(wǎng)格十分困難，全自動化六面體方法的適用性還在逐步研究中[6-7]。因此，權(quán)衡利弊，對于高性能數(shù)值模擬計算來說，其運行在具有數(shù)千上萬核的、有超大規(guī)模內(nèi)存的高性能計算機上，采用較成熟的四面體網(wǎng)格生成方法，自動化程度高，能夠快速響應(yīng)計算需求，且超大規(guī)模的高質(zhì)量四面體網(wǎng)格也能夠提高計算精度。目前，基于AFT方法的廣為應(yīng)用的開源四面體網(wǎng)格生成工具包有Netgen、商業(yè)工具包VKI等，均是花費10～15年時間研發(fā)的串行工具，離人們的需求所缺乏的是并行化。

迄今為止，網(wǎng)格并行生成是一項緊密耦合計算幾何、網(wǎng)格數(shù)據(jù)結(jié)構(gòu)、負(fù)載平衡和網(wǎng)格分區(qū)的復(fù)雜技術(shù)，主要有兩類實現(xiàn)途徑：第一類是直接對網(wǎng)格生成算法進行并行化[8-10]；另一類是將幾何模型分解為子區(qū)域后再進行網(wǎng)格生成，這種方法可復(fù)用串行網(wǎng)格生成算法。例如：非開源并行網(wǎng)格生成工具ParTGen[11]，其并行核數(shù)最多能到500核，遠遠不能滿足萬億次以上規(guī)模計算的需求；ITAPS（interoperable technologies for advanced petascale simulations）[12]是美能源部SciDAC計劃支持的面向千萬億次高性能計算的研究項目，目標(biāo)在于研發(fā)高可用的網(wǎng)格、幾何、場生成的工具，它聯(lián)合了美國各大學(xué)和三大實驗室的力量，在不斷發(fā)展中。

本文重點討論第二類方法的實現(xiàn)，在分區(qū)時，考慮要盡量減少通訊和達到負(fù)載平衡。在并行計算領(lǐng)域，串行的分區(qū)程序Metis[13]、Zoltan[14]和它們的并行版都被廣泛使用，本文在并行策略中，利用這些分區(qū)工具進行粗網(wǎng)格的并行分區(qū)，為大規(guī)模網(wǎng)格并行生成打下了基礎(chǔ)。

2　網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

非結(jié)構(gòu)四面體網(wǎng)格并行生成工具的輸入是Brep格式描述的幾何區(qū)域，而輸出則是直接對接JAUMIN的J網(wǎng)格文件描述格式。在程序內(nèi)部，緊湊的網(wǎng)格數(shù)據(jù)結(jié)構(gòu)為網(wǎng)格并行生成提供了支撐，在單個并行計算節(jié)點上，它采用了局部信息存儲方式，四面體-幾何頂點、三角面-四面體、三角面-幾何頂點等多種鄰接關(guān)系映射。其中，四面體-頂點鄰接關(guān)系包含所有網(wǎng)格單元，而其余鄰接關(guān)系只保存在幾何體和分界面上的三角面的鄰接信息（如圖1所示）。在單個并行計算節(jié)點上，還有描述原始幾何區(qū)域的幾何模型。

Fig.1　Mesh connected relations in memory圖1　內(nèi)存中的網(wǎng)格鄰接關(guān)系

3　AFT并行四面體網(wǎng)格生成方法

3.1　并行策略

串行AFT四面體網(wǎng)格生成過程中，首先基于幾何形狀生成表面三角形網(wǎng)格，再從表面三角形網(wǎng)格出發(fā)生成體網(wǎng)格?？紤]AFT網(wǎng)格生成算法并行化問題時，描述并發(fā)性的一種方式是將每個計算節(jié)點上的網(wǎng)格生成定義為一個任務(wù)，任務(wù)間僅有的共享數(shù)據(jù)是輸入的模型，因為要訪問的輸入模型是只讀的，所以任務(wù)間彼此獨立，從而能夠高效地使用任意合理數(shù)目的計算節(jié)點，節(jié)點間通信開銷幾乎為零。

為實現(xiàn)并行生成，本文采用先生成粗網(wǎng)格，再進行并行細(xì)化的方法，先后采用了兩種并行策略。第一種策略如圖2所示，首先直接針對Brep格式的輸入模型進行劃分，劃分模型的方法采用KD-Tree方法，生成若干塊子模型；然后按順序給子模型生成表面網(wǎng)格，且各子模型塊相鄰面上的表面網(wǎng)格保持一致，將子表面網(wǎng)格分配到各計算節(jié)點，進行后續(xù)的體網(wǎng)格生成。

Fig.2　Parallel strategy by geometry partition method圖2　基于幾何分區(qū)方法的并行策略

這種方法并行效率高，可以將模型分成若干子模型，如圖3所示。但幾何計算不能達到負(fù)載均衡，且網(wǎng)格尺寸由初始尺寸決定，可擴展性不好。

Fig.3　Example of geometry partition圖3　幾何分割模型示意

第二種方法如圖4所示，首先針對輸入模型生成粗網(wǎng)格，然后采用串行的Metis或者Zoltan工具，根據(jù)粗網(wǎng)格的幾何關(guān)聯(lián)性進行分區(qū)，得到子區(qū)域網(wǎng)格（如圖5所示），每個子區(qū)域網(wǎng)格具有幾乎相同數(shù)量的節(jié)點和單元；接著將子區(qū)域網(wǎng)格分配到各計算節(jié)點，計算子區(qū)域的表面網(wǎng)格，不斷重復(fù)網(wǎng)格精細(xì)化步驟，直到網(wǎng)格尺寸達到要求，再生成體網(wǎng)格。這種方法并行效率高，幾何計算負(fù)載均衡性更好，能夠不斷精細(xì)化網(wǎng)格，可擴展性好。

Fig.4　Parallel strategy by coarse mesh partition method圖4　基于粗網(wǎng)格分區(qū)方法的并行策略

Fig.5　Example of coarse mesh partition圖5　粗網(wǎng)格分區(qū)示意

3.2　貼體加密

在本文采用的并行策略中，表面網(wǎng)格的精細(xì)化加密采用了貼體加密方法。

一般有代表性的表面網(wǎng)格細(xì)分加密方法有Loop[15]算法、Doo-Sabin算法、Catmull-clark算法和改進型的蝶形算法，這些算法都是通過逼近或者插值的方法對表面網(wǎng)格進行不斷的光滑，并沒有考慮表面網(wǎng)格所對應(yīng)的輸入模型的幾何形狀。

貼體加密方法示意如圖6所示。其特點是對表面三角形網(wǎng)格進行邊中點加密時，采用向分區(qū)模型的外表面進行投影的方法求得邊中點P對應(yīng)的加密點P′。

Fig.6　Surface mesh and after two times refinement圖6　表面網(wǎng)格和一次、二次加密后的表面網(wǎng)格

由于表面三角形網(wǎng)格包括了輸入模型的表面網(wǎng)格和分區(qū)產(chǎn)生的分界面網(wǎng)格，針對一個表面三角形單元的3個頂點的不同位置情況采用不同的加密點位置計算方法：當(dāng)只有兩個頂點位于輸入模型表面網(wǎng)格時，且它們組成的邊不為其他表面三角形單元共享時，用最短距離法求得加密點P′；當(dāng)兩個頂點位于分界面網(wǎng)格上，P′即為P；當(dāng)兩個頂點位于輸入模型表面網(wǎng)格，一個頂點位于分界面網(wǎng)格時，采用法線方向的最短距離法求得加密點P′。由于加密原則一致，位于不同計算節(jié)點的相鄰分界面上的網(wǎng)格保持一致。

3.3　表面網(wǎng)格的并行優(yōu)化

采用Metis進行分區(qū)時，由于不會考慮天然幾何外形的特點，會出現(xiàn)分區(qū)表面網(wǎng)格不夠光滑和碎片化分區(qū)，如圖7所示，最終影響生成的體網(wǎng)格的質(zhì)量。

因此，對圖4中的并行策略進行了改進，增加對分區(qū)表面網(wǎng)格進行優(yōu)化的步驟，如圖8所示，不斷重復(fù)對分區(qū)表面網(wǎng)格進行優(yōu)化和精細(xì)化加密的步驟，直到達到網(wǎng)格尺寸要求為止，進一步提高體網(wǎng)格的質(zhì)量。由于對表面網(wǎng)格優(yōu)化的步驟需要在計算節(jié)點間進行少量的通信，本文實際是借助JAUMIN框架建立了通信關(guān)系。

分區(qū)表面網(wǎng)格優(yōu)化工作由本文作者和湘潭大學(xué)魏華祎老師小組共同完成，主要采用了節(jié)點優(yōu)化算法，通過移動網(wǎng)格節(jié)點的坐標(biāo)來改善節(jié)點周圍三角形單元的形狀。給定一節(jié)點Pi所在的單元片，如圖9所示。

Fig.7　Coarse mesh of aircraft and fragmentations caused by partition圖7　飛行器粗網(wǎng)格及分區(qū)造成的碎片區(qū)域

Fig.8　Adding optimization step to parallel strategy圖8　并行策略中增加優(yōu)化步驟

Fig.9　Triangular element圖9　三角形單元片

節(jié)點優(yōu)化算法實現(xiàn)了CPT、ODT和CVT共3種移點方法。

CPT方法采用Pi鄰接單元的重心的平均來替代Pi坐標(biāo)。設(shè)與Pi相鄰單元的重心為bj，則：

ODT方法是用Pi鄰接單元的外接圓圓心的平均來替代Pi坐標(biāo)。設(shè)與Pi相鄰節(jié)點為Cj，則：

把Pi鄰接單元的外接圓圓心依次連接起來，就會形成Pi的Voronoi區(qū)域Vi，用Vi的質(zhì)心代替Pi，就是CVT方法，即：

節(jié)點移動過程中，移動方向可分解為沿曲面切線方向和法線方向的移動。其中，沿切線方向的移動改善單元的形狀和尺寸，而沿法線方向的移動改善表面網(wǎng)格的光滑性。但為了保持輸入模型的幾何特征，對位于輸入模型表面上的網(wǎng)格點要特殊處理，只能沿著輸入模型的表面或特征曲線的切線方向移動，移動之后還必須采用3.2節(jié)提到的投影方法映射回模型表面。綜合評價上述移點方法，CPT方法比較穩(wěn)定，適應(yīng)性更強，因為它是重心的平均，所以不會移出單元片之外，導(dǎo)致無效單元出現(xiàn)；CVT和ODT穩(wěn)定性次之，但ODT的優(yōu)化速度更快。

優(yōu)化后的表面網(wǎng)格如圖10所示。

Fig.10　Surface mesh and optimized surface mesh圖10　表面網(wǎng)格和優(yōu)化后的表面網(wǎng)格

表面網(wǎng)格優(yōu)化方法解決了分區(qū)表面網(wǎng)格不夠光滑的問題，但是對碎片化分區(qū)問題仍無法有效解決。

4　應(yīng)用測試

基于本文的AFT四面體并行網(wǎng)格生成技術(shù)，在浪潮TS10K集群進行了同一基準(zhǔn)模型的強可擴展測試，精細(xì)化加密次數(shù)均為2，該集群有64個計算節(jié)點，每個計算節(jié)點內(nèi)存為64 GB。測試結(jié)果如表1和圖11所示。

Table 1　ParallelAFT tetrahedral mesh generation test表1　并行AFT四面體網(wǎng)格生成測試

Fig.11　Strong extensible test圖11　強可擴展測試

從圖11中可看到幾乎可以達到線性加速，具有良好的可擴展性，在64節(jié)點并行時，甚至出現(xiàn)了輕微的超線性。分析原因在于模型越復(fù)雜，在分區(qū)越多時，抽取表面網(wǎng)格的速度非線性提高，而模型不復(fù)雜時，則不會存在這一現(xiàn)象，說明該方法更適合復(fù)雜模型。本文利用該方法，在千萬億次高性能計算機上進行了實際應(yīng)用，實現(xiàn)了針對神光III靶球裝置的實際工程模型的高精度四面體網(wǎng)格生成，達到在1 min內(nèi)8 192核并行生成2.2億自由度網(wǎng)格，如圖12所示。

Fig.12　Shengguang III target ball device and local details of mesh generation圖12　神光III靶球裝置和生成網(wǎng)格的局部細(xì)節(jié)

5　結(jié)束語

本文的無縫對接JAUMIN的AFT四面體網(wǎng)格生成并行方法，有效地復(fù)用了串行AFT四面體網(wǎng)格生成算法，實現(xiàn)了針對復(fù)雜模型的數(shù)億量級規(guī)模以上的高精度四面體網(wǎng)格的快速生成，可擴展性良好。本文方法已集成到中國工程物理研究院計算機應(yīng)用研究所和高性能數(shù)值模擬軟件中心聯(lián)合研制的大規(guī)模復(fù)雜數(shù)值模擬的前處理軟件SuperMesh V1.5[16]版中。下一步，還需要對碎片化分區(qū)問題和分區(qū)表面網(wǎng)格的拓?fù)鋬?yōu)化算法開展研究。

[1]Mo Zeyao.Progress on high performance programming framework for numerical simulation[J].E-Science Technology&Application,2005,6(4):11-19.

[2]Lo S H.A new mesh generation schema for arbitrary planar domains[J].International Journal for Numercial Methods in Engineering,1985,21(8):1403-1426.

[3]L?hner R.Progress in grid generation via the advancing front technique[J].Engineering with Computers,1996,12(3):186-210.

[4]Alleaume A,Francez L,Loriot M,et al.Large out-of-core tetrahedral meshing[C]//Proceedings of the 16th International Meshing Roundtable,Seattle,Oct 14-17,2007.Berlin:Springer,2008:461-476.

[5]Baker T J.Developments and trends in three-dimensional mesh generation[J].Applied Numerical Mathematics,1989,5(4):275-304.

[6]Owen S J,Shelton T R.Evaluation of grid-based hex meshes for solid mechanics[J].Engineering with Computers,2015,31(3):529-543.

[7]Chen Jianjun,Xiao Zhoufang,Cao Jian,et al.Automatic hexahedral mesh generation for many-to-one sweep volumes[J].Journal of Zhejiang University:Engineering Science,2012,46(2):274-279.

[8]Antonopoulos C D,Blagojevic F,Chernikov A N,et al.A multigrain Delaunay mesh generation method for multicore SMT-based architectures[J].Journal of Parallel and Distributed Computing,2009,69(7):589-600.

[9]D3D Mesh Generator[EB/OL].(2001-05-20)[2016-08-12].http://mesh.fsv.cvut.cz/～dr/d3d.html.

[10]Tendulkar S,Beall M,Shephard M S,et al.Parallel mesh generation and adaptation for CAD geometries[EB/OL].(2012-10-12)[2016-09-05].http://www.scorec.rpi.edu/REPORTS/2011-2.pdf.

[11]Ivanov E,Gluchshenko O,Andr? H,et al.Parallel software tool for decomposing and meshing of 3D structures[EB/OL].(2008-12-21)[2016-08-12].http://www.itwm.fraunhofer.de/fileadmin/ITWM-Media/Zentral/Pdf/Berichte_ITWM/2007/bericht110.pdf.

[12]Shephard M S.Interoperable technologies for advanced petascale simulations[EB/OL].(2010-02-05)[2016-07-15].http://www.osti.gov/scitech/biblio/971531#.

[13]Karypis G,Kumar V.Multilevelk-wayhypergraph partitioning[C]//Proceedings of the 36th Conference on Design Automation,New Orleans,Jun 21-25,1999.New York:ACM,1998:343-348.

[14]Zoltan:paralllel partioning,load balancing and data-management services[EB/OL].(2007-05)[2016-09].http://www.cs.sandia.gov/Zoltan.

[15]Loop C.Smooth subdivision surfaces based on triangles[D].Salt Lake City:University of Utah,1987.

[16]SuperMesh[EB/OL].(2016-07-12)[2016-09-10].http://www.caep-scns.ac.cn/SuperMesh.php.

附中文參考文獻：

[1]莫則堯.高性能數(shù)值模擬編程框架研究進展[J].科研信息化技術(shù)與應(yīng)用,2015,6(4):11-19.

[7]陳建軍,肖周芳,曹建,等.多源掃掠體全六面體網(wǎng)格自動生成算法[J].浙江大學(xué)學(xué)報:工學(xué)版,2012,46(2):274-279.