面向JAUMIN的并行AFT四面體網格生成*

鄭　澎，方　維，徐　權，冷玨琳，熊　敏，于長華

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

2.北京應用物理與計算數學研究所，北京 100088

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

1　引言

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

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

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

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

2　網格數據結構

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

Fig.1　Mesh connected relations in memory圖1　內存中的網格鄰接關系

3　AFT并行四面體網格生成方法

3.1　并行策略

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

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

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

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

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

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

Fig.4　Parallel strategy by coarse mesh partition method圖4　基于粗網格分區方法的并行策略

Fig.5　Example of coarse mesh partition圖5　粗網格分區示意

3.2　貼體加密

在本文采用的并行策略中，表面網格的精細化加密采用了貼體加密方法。

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

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

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

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

3.3　表面網格的并行優化

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

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

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

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

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

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

節點優化算法實現了CPT、ODT和CVT共3種移點方法。

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

ODT方法是用Pi鄰接單元的外接圓圓心的平均來替代Pi坐標。設與Pi相鄰節點為Cj，則：

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

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

優化后的表面網格如圖10所示。

Fig.10　Surface mesh and optimized surface mesh圖10　表面網格和優化后的表面網格

表面網格優化方法解決了分區表面網格不夠光滑的問題，但是對碎片化分區問題仍無法有效解決。

4　應用測試

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

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

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

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

Fig.12　Shengguang III target ball device and local details of mesh generation圖12　神光III靶球裝置和生成網格的局部細節

5　結束語

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

[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]莫則堯.高性能數值模擬編程框架研究進展[J].科研信息化技術與應用,2015,6(4):11-19.

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