Level set函數快速步進重構并行算法的改進

黃筱云, 董國海, 常佳夫, 蔣學煉

(1.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連116024； 2.長沙理工大學 水利工程學院，湖南 長沙 410114; 3.天津城建大學 天津市軟土特性與工程環境重點實驗室，天津 300384)

?

Level set函數快速步進重構并行算法的改進

黃筱云1,2, 董國海1, 常佳夫2, 蔣學煉3

(1.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連116024； 2.長沙理工大學 水利工程學院，湖南 長沙 410114; 3.天津城建大學 天津市軟土特性與工程環境重點實驗室，天津 300384)

為提高level set函數快速步進重構過程的并行計算效率，本文提出一種改進的分區并行重構算法。與原有分區并行算法相比，優化了子區域間的同步方案，縮短了level set函數并行重構的計算時間。運用OpenMP多線程技術，建立了相應的并行計算模型，實現了圓球、圓環管和啞鈴等值面并行重構。并行重構數值結果表明:只要子區域均分初始表面邊界，level set函數全局或局部并行重構均具有良好加速比，8線程的最大加速比可接近6。

level set函數；快速步進法；重構；并行算法; 多線程技術；OpenMP多線程技術

分區并行計算是將整個計算區域分割成多個子區域，并將各子區域的計算分配給不同線程執行。作為粒子level set法的重要環節，level set函數并行重構的實現是粒子level set法并行化的關鍵。在這方面，筆者曾提出過一種level set函數分區并行重構算法[13]，有效地縮短了計算時間。本文的主要工作是在該算法的基礎上優化并行同步方案，以進一步提高重構效率，建立基于OpenMP多線程技術的并行計算新模型，并通過圓球、圓環管和啞鈴三個等值面重構算例，分析和討論并行計算新模型的精度和效率。

1 快速步進法

(1)

上述格式的特點是網格節點的選擇依賴于節點φ值大小，φ值較小的節點決定φ值較大的節點。即當確定緊鄰交界面網格節點的φ值后，便可從邊界開始逐步向外構造出φ>0區域內其他網格節點的φ值。若改變區域內網格節點上φ值的符號，亦可構造出φ<0區域內網格節點上的φ

值。具體算法如下：

1)標記交界面周圍的已知φ值的網格節點為Alive節點，與Alive節點相鄰的其他節點為Close節點，剩余節點為Far節點，如圖1(a)所示；

2)按式(1)計算Close節點上的φ值；

3)在Close節點集中找到φ值最小的節點A，標記該節點為Alive節點，如圖1(b)所示；

4)重新計算與節點A相鄰的所有節點的φ值，并將其中的Far節點B和C標記為Close節點，如圖1(c)所示；

5)從Close節點集中刪除節點A；

6)若Close節點集為空，計算結束，否則返回步驟3。

圖1 快速步進法構造φ值的過程Fig.1 Construction of φ value by fast marching method

2 并行快速步進法的改進

2.1 計算區域劃分與網格布置

要發揮并行計算的優勢，通常將計算區域劃分成多個大小一致的子區域，以保證每個線程分配的計算量接近。為了方便相鄰區域間交換數據，各子區域外需要布置虛擬網格。例如圖2，整個計算區域平均分成2個子區域，子區域1在右邊界外布置一排虛擬網格，用于備份子區域2左邊界網格上的函數值，同理，子區域2在左邊界外布置一排虛擬網格，用于備份子區域1右邊界網格上的函數值；節點A、C和虛擬網格節點C′的坐標位置一致，子區域1的虛擬節點C′用于接收子區域2上節點A的信息，同樣，子區域2的虛擬節點也將接收子區域1右邊界節點的信息。

圖2 計算區域劃分與虛擬網格布置Fig.2 Subdivision of computation domain and arrangement of ghost mesh

2.2 原有并行算法

原有的分區并行的基本算法[13]如下：

1)對于任一子區域，運用快速步進法構造交界面外所有非Alive節點的φ值；

2)同步子區域；

3)滿足收斂條件，計算結束，否則，返回步驟1。

在子區域同步過程中，本區域的虛擬網格節點首先接收從相鄰區域邊界節點傳遞來的φ值，再檢查虛擬網格節點的φ值是否小于邊界節點的φ值，若成立，將虛擬網格節點的φ值賦給邊界節點，并定義所有賦予值中最小者為φmin，最后，重新計算本區域內所有大于φmin的網格節點上φ值。由于要重新計算所有大于φmin的網格節點，這使得步驟1中部分非Alive節點上的φ值構造可能是無用的。因此，可將同步過程嵌入到快速步進法循環中，通過及時傳遞邊界節點信息，減少網格節點φ值的無效計算次數，提高計算效率。

2.3 改進的并行算法

在改進的并行化快速步進法中，一方面將新標記為Alive邊界節點的函數值及時傳遞到相鄰子區域的虛擬網格節點上；另一方面在每次循環中判斷相鄰子區域計算結果是否對本子區域計算產生影響，即每次循環都檢查子區域虛擬網格節點φ值是否被更新，且小于對應的邊界節點φ值。具體算法如下：

1)執行串行快速步進法中步驟1和2；

2)檢查是否存在新標記為Alive的虛擬網格節點，且該節點上φ值小于位置重合的邊界節點上的φ值；若這種虛擬節點存在，則對應邊界節點標記為Close，同時將虛擬網格節點上φ值賦給該邊界節點，然后將本子區域內所有不小于該φ值的Alive節點重新標記為Close，與交界面相鄰的節點除外；

3)若Close節點集非空，執行串行快速步進法中步驟3)～5)。在步驟4中，若重新計算節點為邊界節點，其φ值應取計算值和原值中的較小者；

4)若上一步中刪除的Close節點為邊界節點，將該邊界節點的φ值傳遞給相鄰子區域的虛擬網格節點，并將虛擬網格節點標記為Alive；

5)若Close節點集非空，返回步驟2；

6)若所有子區域Close節點集皆為空，計算結束；若存在新標記為Alive的虛擬網格節點，則返回步驟2。

圖3給出了子區域間的同步過程。圖3(a)中，節點A和B分別是子區域2和1中φ值最小的Close節點；當節點A被標記成Alive節點后，該節點的相關信息被傳遞給子區域1，并備份到虛擬節點C′上，同時節點B被標記為Alive(圖3(b))；圖3(c)中，子區域1的網格節點將根據虛擬節點C′上φ值進行重新標記，節點B又被標記回Close。由于節點C是子區域1中φ值最小的Close節點，C節點將被標記成Alive節點，并重新計算與C相鄰的非Alive節點上的φ值，同時標記相鄰的Far節點為Close(圖3(d))。

圖3 子區域間同步過程Fig.3 Synchronization between sub-domains

3 并行計算實例分析

3.1 圓球

計算區域為100×100×100，圓球位于整個區域中心，球體半徑為25，如圖4(a)所示。計算過程中，整個計算區域均分成8個大小相等的子區域。圖4(b)給出了φ=-10和φ=10等值面的并行計算結果。表1給出了φ=10等值面所包圍體積的數值結果以及整個區域φ值的計算誤差，其中L1誤差按照文獻[14]計算，即

(2)

在粒子level set方法中，level set方程求解和level set函數重構過程只需在交界面周圍一定范圍內執行，即可保證交界面捕捉的準確性。圖5(a)比較了全局和局部(|φ|<3)并行化快速步進法的加速比?？梢钥闯?，無論全局還是局部操作，并行計算模型均獲得良好的加速比，8線程下全局計算的加速比接近6，而局部計算的加速比也大于5。

圖 4 圓球等值面重構Fig.4 Isosurface reconstruction of sphere

圖5 圓球等值面重構的加速比Fig.5 Speedup of isosurface reconstruction of sphere

Table 1 Computational loss and error of the parallel isosurface reconstruction of sphere

分辨率體積損失/%L1誤差階數真值179594.4———50174782.92.76.24×10-1—100177135.41.43.13×10-10.99200178319.10.71.59×10-10.97

在并行計算過程中，分區并行所需子區域間信息交換和節點重新計算等額外開銷是影響并行計算效率的關鍵。定義并行信息傳遞參數Fc和并行節點重演參數Fr分別為

(3)

(4)

式中：Nt為子區域間信息傳遞總次數，M為計算區域內總節點數,Mr為節點重新計算總次數。

圖6比較了不同線程數下的Fc和Fr值(球心位置[50,50,50])?？梢钥闯觯篎c和Fr都隨著線程數的增加而變大，說明加速比的增長幅度不是恒定的，而是隨線程數增加呈減小的趨勢；Fc和Fr值小于0.2，表明該并行狀態的額外開銷對加速比的影響較??；相比全局并行重構，局部并行重構的Fc和Fr值更大，表明其額外開銷所占的比例更高，使得加速比的增長幅度低于全局并行重構。

圖5(b)給出了位于[25,25,25]的圓球在不同線程下全局和局部level set函數重構的加速比，此時，8線程下全局計算的加速比只達到2，而局部操作甚至不到1。圖7比較了此狀態下的Fc和Fr值?？梢钥闯?，該狀態下的額外計算開銷遠遠高于球心位于[50,50,50]的情況。尤其是Fr值，8線程下局部重構的Fr值甚至超過4，這是導致其加速比隨線程數增加反而減小的主要原因。

圖8比較了8線程下圓球在不同位置的加速比。當球心位置位于[25,25,25]時，整個球體完全落入一個子分區中，過量的信息交換與節點重演帶來巨大的額外開銷，所獲得的加速比最??；當圓球向區域中心逐漸靠近時，各子區域所獲得初始表面面積趨于相等，所獲得的加速比隨之增加，直至最大。

圖6 圓球等值面重構的Fc和Fr(球心位置 [50,50,50])Fig.6 Fc and Fr for isosurface reconstruction of sphere(centre located at [50,50,50])

圖7 圓球等值面重構的Fc和Fr值(球心位置[25,25,25])Fig.7 Fc and Fr for isosurface reconstruction of sphere (center located at [25,25,25])

圖8 8線程下不同位置圓球等值面重構的加速比Fig.8 Speedup of isosurface reconstruction of sphere located at different positions under eight threads

3.2 圓環管

計算區域為100×100×100，圓環管放置在計算區域中心(圖9(a))，其表達式如下

(5)

其中，外徑R=20，內徑r=5。計算過程中，整個計算區域可分成8個大小相等的區域。圖9(b)給出了φ=5和φ=10等值面的計算結果。表2則比較了不同網格分辨率下φ=10等值面所包圍體積的數值結果及φ值重構的計算誤差。

圖9 圓環管等值面重構Fig.9 Isosurface reconstruction of circular cube

Table 2 Computational loss and error of parallel isosurface reconstruction of circular cube

分辨率體積損失/%L1誤差階數真值88826.4———5084721.34.68.26×10-1—10086429.52.74.46×10-10.8920087535.81.42.31×10-10.95

圖10比較了全局和局部(|φ|<3)重構下并行模型的加速比，其中2線程和4線程的局部并行重構加速比取三種計算區域均分方案的平均。這是因為采用2線程和4線程進行局部并行重構時，邊界上需要信息傳遞的邊界節點總數隨計算區域均分方式不同而有所差異。例如，在線程數為2的條件下，若按y方向均分計算區域，所要計算的邊界節點總數要大于其他兩種情況，從表3可以看出，按y方向均分計算區域的Fc和Fr值大于其他兩種均分方式，這也使得其加速比略小于其他兩種方式。

3.3 啞鈴曲面

計算區域為100×100×100，啞鈴曲面的具體表達式如下

(6)

其中，a=25(見圖11(a))。整個計算區域被平均分成8個子區域。圖11(b)為φ=10和φ=20等值面的計算結果。表4比較了不同網格分辨率下φ=10等值面所包圍體積的數值結果及φ值重構的計算誤差。

圖10 圓環管等值面重構的加速比Fig.10 Speedup of isosurface reconstruction of circular cube

表3 2線程下不同區域劃分對的等值面局部并行重構效率的影響

Table 3 Effect of domain division on the efficiency of the local parallel isosurface reconstructions under two threads

加速比FcFrx1.873.39×10-21.60×10-2圓環管y1.601.18×10-12.80×10-2z1.873.29×10-21.55×10-2x1.900.97×10-21.02×10-2啞鈴曲面y1.779.27×10-22.37×10-2z1.769.42×10-22.11×10-2

圖 11 啞鈴等值面重構Fig.11 Isosurface reconstruction of dumbbell

Table 4 Computational loss and error of parallel isosurface reconstruction of dumbbell

分辨率體積損失/%L1誤差階數真值56767.5———5053077.86.51.12—10054665.43.75.99×10-10.9020055642.82.03.31×10-10.86

在本算例中，按x方向均分計算區域時需要信息傳遞的邊界節點數最少，因此，其2線程下的加速比要略高于其他兩種情況(表3)。圖12為全局和局部(|f|<3)并行計算的加速比,同樣，2線程和4線程下局部并行重構的加速比取三種計算區域均分方案的平均?？梢钥闯?，8線程的加速比仍大于4。

圖12 啞鈴等值面重構的加速比Fig.12 Speedup of isosurface reconstruction of dumbbell

4 結論

1)改進的快速步進并行重構方法能夠顯著縮短計算時間，且基本保持1階精度，在8線程的條件下，最大加速比能夠接近6；

2)相同線程數下，子區域對初始表面邊界的劃分將影響并行重構的加速比，均分表面可獲得較高的加速比；

3)局部分區并行重構應盡量減少子區域邊界上需要計算的節點數目，這樣子區域間信息傳遞次數和節點重新計算的次數將相應地減少，而計算效率也隨之提高。

[1]OSHER S, SETHIAN J A. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[J]. Journal of computational physics, 1988, 79(1): 12-49.

[2]SUSSMAN M, SMEREKA P, OSHER S. A level set approach for computing solutions to incompressible two-phase flow [J]. Journal of computational physics, 1994, 114(2): 146-159.

[3]JIANG G S, SHU C W. Efficient implement of weighted ENO schemes [J]. Journal of computational physics, 1996, 126(1): 202-228.

[4]JIANG G S, PENG D P. Weighted ENO schemes for Ham-

ilton-Jacobi equations [J]. SIAM journal on scientific computing, 2000, 21(6): 2126-2143.

[5]ENRIGHT D, FEDKIW R. A hybrid particle level set method for improved interface capturing [J]. Journal of computational physics, 2002, 183(1): 83-116.

[6]ENRIGHT D, LOSASSOM F. A fast and accurate semi-Lagrangian particle level set method [J]. Computers and structure, 2005, 83(6): 479-490.

[7]JIANG L, LIU F B, CHEN D R. A fast particle level set method with optimized particle correction procedure for interface capturing[J]. Journal of computational physics, 2015, 299: 804-819

[8]SETHAIN J. Fast marching methods[J]. SIAM review, 1999, 41(2): 199-235.

[9]STRAIN J. Semi-Lagrangian methods for level set equations[J]. Journal of computational physics, 1999, 151(2):498-533.

[10]黃筱云, 李紹武, 夏波. 一種新型三維水流數值模型 [J]. 海洋學報, 2010, 32(6): 167-173.

HUANG Xiaoyun, LI Shaowu, XIA Bo. A new three dimensional numerical model for stream[J]. Acta oceanologica sinica, 2010, 32(6): 167-173.

[11]HUANG X Y, LI S W. A two-dimensional numerical wave flume based on SA-MPLS method[J]. Acta oceanologica sinica, 2012, 31(3): 18-30.

[12]LI S W, ZHUANG Q, HUANG X Y, et al. 3D simulation of flow with free surface based on adaptive octree mesh system[J]. Transactions of Tianjin University, 2015, 21(1): 32-40.

[13]黃筱云, 董國海, 趙利平，等. Level set函數重新初始化的并行快速步進方法[J]. 哈爾濱工程大學學報, 2016, 37(5): 666-671.

HUANG Xiaoyun, DONG Guohai, ZHAO Liping, et al. A parallelized fast marching method for reinitialization of level set function[J]. Journal of Harbin Engineering University, 2016, 37(5): 666-671.

[14]SUSSMAN M, FATEMI E. An efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow[J]. SIAM journal on scientific computing. 1999, 20(4): 1165-1191.

本文引用格式：

黃筱云, 董國海, 常佳夫, 等. Level set函數快速步進重構并行算法的改進[J]. 哈爾濱工程大學學報， 2017， 38(6)： 836-842.

HUANG Xiaoyun， DONG Guohai， CHANG Jiafu， et al. Improvement of parallel fast marching method for reconstruction of level set function[J]. Journal of Harbin Engineering University， 2017， 38(6)： 836-842.

Improvement of parallel fast marching method for reconstruction of level set function

HUANG Xiaoyun1,2， DONG Guohai1， CHANG Jiafu2， JIANG Xuelian3

(1.State Key Laboratory of Coastal and offshore Engineering, Dalian University of Technology, Dalian 116024, China; 2.School of Hydraulic Engineering, Changsha University of Science and Technology, Changsha 410114, China; 3.Tianjin Key Laboratory of Soft Soil Characteristics & Engineering Environment, Tianjin Chengjian University, Tianjin 300384, China)

To raise the parallel computation efficiency of the fast marching method for reconstruction of the level set function, an improved parallel algorithm based on domain decomposition was developed. Compared with the previous parallel algorithm, the present algorithm optimizes the sub-domain synchronization scheme and reduces the calculation time for parallel reconstruction of the level set function. Using the OpenMP multi-thread technique, a parallel computational model was established, and the parallel isosurface reconstructions of a sphere, circular ring pipe, and dumbbell were enforced. The numerical calculation results on parallel reconstruction show that a good speed-up ratio was obtained in either global or local parallel reconstruction of the level set function, and the maximum speed-up ratio approached six times under eight threads on the condition that initial surface boundary was equally divided by sub-domains.

level set function; fast marching method; reconstruction; parallel algorithm；multi-thread technique；OpenMP multi-thread technique

2016-04-18. 網絡出版日期：2017-04-05.

國家自然科學基金項目(51109018, 51309036)；中國博士后科學基金項目(2014M561230)；湖南省自然科學基金項目(2015JJ2006)；天津市自然科學基金項目(14JCYBJC22100).

黃筱云(1980-), 男, 講師, 博士.

黃筱云，E-mail: Xiaoyun.huang@csust.edu.cn.

10.11990/jheu.201604048

http://www.cnki.net/kcms/detail/23.1390.u.20170405.1558.004.html

TV131.2

文章編號：1006-7043(2017)06-0836-07