Level set函數重新初始化的并行快速步進法

黃筱云, 董國海, 趙利平, 程永舟

(1. 大連理工大學 海岸和近海工程國家重點實驗室, 遼寧 大連 116024；2. 長沙理工大學 水利工程學院, 湖南 長沙 410004；

Level set函數重新初始化的并行快速步進法

黃筱云1,2,3, 董國海1, 趙利平2, 程永舟2

(1. 大連理工大學 海岸和近海工程國家重點實驗室, 遼寧 大連 116024；2. 長沙理工大學 水利工程學院, 湖南 長沙 410004；

3. 河海大學 水文水資源與水利工程科學國家重點實驗室, 江蘇 南京 210098)

摘要：為提高level set函數重新初始化的計算效率，基于分區并行思想，提出一種快速步進法的并行策略，實現level set函數的快速并行重新初始化。通過對圓球、五葉管和圓環管等算例的level set函數重新初始化，討論了新并行算法的準確性和效率。結果表明，與串行快速步進法相比，并行算法保留了串行算法的精度，仍基本保持在1階左右，同時顯著減少了重新初始化的計算時間，特別在8線程條件下，所獲的最佳加速比能夠達到5。

關鍵詞：level set函數；重新初始化；快速步進法；并行；分區；并行算法；加速比

網絡出版地址：http://www.cnki.net/kcms/detail/23.1390.U.20160411.1031.038.html

基于歐拉方法模擬波浪傳播過程，要求精確捕捉或追蹤自由表面位置。與其他界面捕捉或追蹤方法不同，level set方法是通過計算level set函數φ的分布變化隱式地表達自由表面位置變化過程，在模擬波浪破碎、融合等復雜拓撲變化過程方面具有無法比擬的優勢[1-3]。在采用level set方法追蹤自由表面的具體實施步驟中，必須保證每一時刻的level set函數滿足距離函數條件，需要對level set函數進行重新初始化。目前，實現level set函數重構的方法有兩種：迭代法和快速步進法。前者是通過計算特定方程穩定解的方式來完成重新初始化[4-5]；后者則順風地從界面向外構造出level set函數值，相比前者，后者不會造成界面非物理移動，且計算效率較高[6-9]。隨著波浪運動研究的深入，特別是三維波浪與建筑物作用數值仿真開展，計算時間通常都比較可觀，引入并行計算技術則可以有效地減少數值模擬時間。在當前level set方法的并行化處理中，重新初始化過程通常采用迭代方式執行[10-12]，而有關采用快速步進法重新初始化的level set并行算法則很少，因此，有必要開展基于快速步進法的level set函數并行重新初始化研究。本文在文獻[3]基礎上，提出一種快速步進法的并行化策略，以提高level set函數重新初始化的效率，并通過圓球、五葉管和圓環管等算例的level set函數重新初始化結果分析該并行化快速步進法的計算精度，探討多線程并行計算效率。

1Level set方法概述

Level set方法采用的控制方程為

(1)

為保證式(1)計算的level set函數具有符號距離函數的特征，即

(2)

需要對式(1)的計算結果進行重新初始化。迭代重新初始化通過計算下面方程穩定解來完成：

(3)

式中：S(φ)為光滑函數，τ為虛擬時間；而快速步進法則是運用一種迎風格式對方程(2)的邊值問題直接求解。

2快速步進法的基本原理

快速步進方法通常采用一階迎風格式來對式(2)進行求解和推進，其數值離散格式具體如下所示：

(4)

式中△h為單元網格尺寸。可以看出，該格式對網格節點的選擇類似于根據流的信息域選擇節點，這種節點選擇方式可以保證信息始終向一個方向傳遞，即從較小值的節點向較大值的節點傳遞。當獲得交界面周圍網格節點的φ值后，就可以向外構造出φ> 0區域內網格節點上的φ值，若改變區域內網格節點上φ值的符號，就可構造出φ< 0區域內網格節點上的φ值。快速步進法的具體實現步驟可以歸納如下：

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

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

3)在Close節點集中找到φ值最小的節點A，標記該節點為Alive節點，而與之相鄰Far節點則標記為Close節點；

4)重新計算與φ值最小節點相鄰所有Close節點上的φ值；

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

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

如圖1所示，交界面周圍節點標記為Alive節點，與之相鄰的為Close節點，其余均為Far節點；計算所有Close節點上的φ值，找到φ值最小的節點A；標記其周圍Far節點為Close節點；將節點A作為Alive節點，重新計算其周圍Close節點上的φ值。

圖1　快速步進法重新初始化φ值的過程Fig.1　Reinitialization of φ by the fast marching method

3快速步進法并行策略

作為一種主流的并行策略，分區并行是將整個計算區域劃分成多個子區域，并分配給不同線程計算；與單區計算相比，分區并行計算需要在邊界上交換數據，以使計算結果盡量和單區計算的結果一致，要實現快速步進法的分區并行就要考慮步進過程中相鄰區域的影響。并行快速步進算法首先需要將整個計算區域分成多個子區域，然后利用不同的線程來分別執行不同子區域的計算任務，再在子區域邊界上交換數據，保證子區域邊界處節點選擇的迎風性。

3.1子區域劃分

要實現節點選擇的迎風性，需要在計算過程中保持子區域之間的信息傳遞。這要求在子區域的劃分過程中，為每個子區域的邊界外布置虛擬網格，以接收相鄰子區域傳遞來的信息。圖2中，整個計算區域被劃分成4個子區域，不同子區域需要被分配不同的線程分別進行計算，每個子區域在與其他子區域相鄰的邊界外布置虛擬網格以接收信息，如子區域1在右邊界和上邊界布置的虛擬網格將分別接收子區域2左邊界網格和子區域3下邊界網格上信息。

圖2　計算網格的劃分以及虛擬網格的布置Fig.2　Partition of computational mesh and arrangement　of ghost mesh

3.2具體并行算法

按照分區并行的基本思路，并行運算的具體步驟如下：

1)對于任何一個子區域，運用快速步進法重新初始化該區域網格Close節點的φ值；

2)進行子區域同步；

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

在并行計算過程中，子區域同步可按下面步驟執行：

①子區域虛擬網格節點接收相鄰區域邊界相應節點上的值；

②比較子區域邊界節點與相應虛擬節點上的φ值，若子區域邊界節點的φ值大于相應虛擬節點，表明該節點會受到相鄰區域的影響，重新標記該節點為Close節點，并將虛擬節點的φ值賦給該節點；

③找出所有受相鄰區域影響邊界節點中φ值最小者，標記為Alive節點，并保存其φ值為φmin；

④重新標記子區域內網格節點，若節點φ>φmin，則該節點標記為Far節點。

在圖3中，相鄰區域右邊界節點的φ值先賦給主區域虛擬網格節點；然后比較主區域邊界節點和虛擬節點的φ值，標記受影響邊界節點為Close節點，并更新其φ值；再找到邊界上φ值最小的受影響節點A，并以節點A的φ值為閾值，重新標記主區域內所有網格節點。

在同步過程中，如果所有子區域受相鄰區域影響邊界節點集均為空，則表明整個區域內網格節點φ值符合迎風特征，計算即可結束。

需要注意的是，由于快速步進算法不變更直接與交界面相鄰的Alive節點狀態，因此，在同步過程中，即使這些Alive節點的φ>φmin，也應該禁止對這些節點重新標記。如圖4所示，相鄰區域節點A'為Alive節點，主區域右邊界上節點A是受相鄰區域的影響節點，同步過程中，節點A的φ值為標記主區域節點的閾值φmin，而Alive節點B的φ>φmin，但由于節點B的φ值是快速步進算法的邊界條件，所以禁止計算流程改變其節點類型。

(a) 虛擬網格節點賦值

(b) 標記受影響的邊界節點

(c) 找到φ值最小的受影響節點A

(d) 重新標記主區域內網格節點圖3　子區域同步過程Fig.3　Synchronization of sub-domain

圖4　節點B禁止重新標記Fig.4　Reflag of B node is forbidden

4算例及分析

為了驗證上文所提出的并行策略的有效性，本節分別選取圓球、五葉管和圓環管3個算例對不同執行條件下快速步進重新初始化進行了數值對比實驗。所有數值算例均采用ThinkCentre M8300t主機進行計算，其處理器為4核8線程的Intel i7-2600，內存大小為16 G。

4.1圓球

本算例設定計算區域大小為1×1×1，分辨率為100×100×100，整個計算區域將劃分成左右兩部分。圓球位于左側計算區域，如圖5所示，其球心位置為[0.3,0.5,0.5]，球體半徑為0.2。

圓球表面的φ值等于零，計算區域內其他位置的φ值可以通過快速步進法獲得。為驗證子區域同步算法的有效性，計算區域分成[0,0.5] ×[0,1] ×[0,1]和[0.5,1] ×[0,1] ×[0,1]兩部分，圓球完全落入左側區域，右側區域內的φ值必須根據左側區域傳遞來的信息來確定。圖6給出分區域并行計算

給出的不同φ值等值面，可以看出，子區域間同步算法能夠有效地將左半區φ值信息傳遞給右半區。

圖5　圓球表面(φ的零等值面)Fig.5　Surface of sphere (isosurface at φ value of zero)

圖6　不同φ值等值面(圓球)Fig.6　Isosurfaces at different φ value (sphere)

表1給出了不同網格分辨率下并行化快速步進法獲得的φ等值面包圍的體積及φ值計算誤差。從表1中可以看出，并行化快速步進法計算得到的φ= 0.1等值面包圍的體積與精確解相差很小，在200×200×200網格分辨率下，體積誤差小于0.1%，且誤差階數為1階。

由于并行化快速步進法獲得的L1誤差量級為10-3，分區域計算時可設定邊界容許誤差為1.0×10-4。將圓球放置在計算區域中心，可比較不同線程數下實現φ值重新初始化所需的CPU計算時間，具體計算結果如表2所示。

表1不同網格分辨率下圓球計算誤差

Table 1Computational error of sphere under differentgrid resolution

分辨率體積損失/%L1誤差階數真值1.13×10-1———50×50×501.09×10-23.718.46×10-3—100×100×1001.11×10-21.944.09×10-31.05200×200×2001.13×10-20.092.11×10-30.96

從表2可以看出，并行化快速步進法能夠有效縮短計算時間，在網格分辨率較大的條件下，加速比能達到5左右，但子區域間同步開銷使得效率會隨著線程數增加而減小，當線程數達到一定數量時，計算時間甚至會略有增加。

表2不同線程數計算耗費CPU時間(圓球)

Table 2CPU time under different number ofthreads (sphere) s

分辨率1線程2線程4線程8線程50×50×500.870.430.260.28100×100×1009.114.512.741.79200×200×20086.0944.0826.1016.56

4.2五葉管

本算例的五葉管由y-z平面上五葉圖形在x方向延長所形成，五葉圖形表達式如下：

(5)

在本算例中，R= 0.2,r= 0.1，計算區域大小為2×1×1，網格分辨率200×100×100。圖7為φ= 0和 0.1等值面。

(a) 虛擬網格節點賦值

(b) 標記受影響的邊界節點圖7　不同φ值等值面(五葉管)Fig.7　Isosurfaces at different φvalue (pentafoil cube)

分辨率體積損失/%L1誤差階數真值7.89×10-1———50×25×257.66×10-12.857.87×10-3—100×50×507.74×10-11.884.64×10-30.76200×100×1007.87×10-10.222.60×10-30.83

表3給出了φ= 0.1等值面包圍的體積以及φ值計算誤差。從中可以看出，一方面，隨著網格分辨率的增加，體積損失逐漸減小，當網格分辨率為200×100×100時，體積損失僅為0.22%；另一方面，盡管五葉管表面某些區域曲率較大，但計算結果誤差階數仍大于0.75。

在x方向按等長度將整個計算區域劃分成多個子區域，分配給不同線程進行計算。表4給出不同網格分辨率下，多線程并行計算所耗費的CPU計算時間。從中可以看出，子區域間同步開銷抵消了雙線程帶來的運算加速，但隨著線程數增加，計算時間進一步縮減，在網格分辨率為200×100×100的條件下，加速比能達到2.5。

表4不同線程數計算耗費CPU時間(五葉管)

Table 4CPU time under different number of threads(pentafoil cube) s

分辨率1線程2線程4線程8線程50×50×500.870.430.260.28100×100×1009.114.512.741.79200×200×20086.0944.0826.1016.56

4.3圓環管

本算例將重新初始化圓環管外φ值分布,計算區域大小為1×1×1，分辨率為100×100×100，為平均分配每個線程計算量，圓環管放置在計算區域中心，其表達式為

(6)

其中，R=0.3,r= 0.05。φ零等值面如圖8所示。圖9(a)給出了φ= 0.1的等值面，圖9(b)～(i)則給出了分區域計算各子區域內φ= 0.1等值面。

圖8　圓環管表面(φ =0等值面)Fig.8　Surface of circular cube (isosurface at φ =0)

分辨率體積損失/%L1誤差階數真值1.33×10-1———50×50×501.28×10-14.116.82×10-3—100×100×1001.30×10-12.563.51×10-30.96200×200×2001.31×10-11.581.81×10-30.92

表5則給出了不同網格分辨率下并行化快速步進法獲得的φ= 0.1等值面包圍的體積及φ值計算誤差。通過對比可以看出，φ= 0.1等值面包圍體積與真值相差很小，體積損失與網格分辨率大小幾乎成線性關系，且并行化快速步進法的計算精度幾乎達到1階。

表6不同線程數計算耗費CPU時間(圓環管)

Table 6CPU time under different number of threads (circular cube)s

分辨率1線程2線程4線程8線程50×50×500.700.730.460.40100×100×1008.398.135.154.01200×200×20078.02149.0471.4553.35

依次在x、z、y方向上平分計算區域，給出不同線程數下實現φ值重新初始化所需的CPU計算時間(表6)。可以看出，隨著網格分辨率變大，子區域間同步開銷也會大幅增加，甚至大于快速步進法本身的開銷，但隨著線程數的增加，加速比能夠達到2。此外，計算區域的劃分方式也會對加速比造成影響(表7)。

表7不同區域劃分策略下的并行計算時間比較

Table 7Comparison of parallelized computation time under different partition strategiess

分辨率線程數x/x,yy/y,zz/z,x200×200×2002149.0482.0180.96493.3163.0171.45

5結論

通過分區并行，提出了一種level set函數重新初始化快速步進并行方法，通過圓球、五葉管和圓環管算例得到如下結論：

1)并行快速步進法能夠有效縮短計算時間，且保持計算精度基本不變，在8線程的條件下，最佳加速比能夠達到5左右；

2)并行快速步進法的重新初始化效率與區域劃分有關，較好的區域劃分能夠使效率提高更為顯著；

3)并行化快速步進法的計算精度會受到表面邊界曲率條件的影響，邊界越光滑，獲得的計算結果精度越接近1階。

參考文獻：

[1]LU Xinhua, ZHANG Xiaofeng, LU Junqing, et al. Numerical simulation of breaking wave generated sediment suspension and transport process based on CLSVOF Algorithm[J]. China ocean engineering, 2014, 28(5): 701-712.

[2]MARKUS D, ARNOLD M, WüCHNER R, et al. A virtual free surface (VFS) model for efficient wave-current CFD simulation of fully submerged structures[J]. Coastal engineering, 2014, 89: 85-98.

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

[4]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(1): 146-159.

[5]WANG Zhaoyuan, YANG Jianming, STERN F. An improved particle correction procedure for the particle level set method[J]. Journal of computational physics, 2009, 288(16): 5819-5837.

[6]RUSSO G, SMEREKA P. A remark on computing distance functions[J]. Journal of computational physics, 2000, 163(1): 51-67.

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

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

[9]ENRIGHT D, LOSASSO F, FEDKIW R. A fast and accurate Semi-Lagrangian particle level set method[J]. Computers & structure, 2005, 83(6/7): 479-490.

[10]SUSSMAN M. A parallelized, adaptive algorithm for multiphase flows in general geometries[J]. Computers & structures, 2005, 83(6/7): 435-444.

[11]WANG Kai, CHANG A, KALE L V, et al. Parallelization of a level set method for simulating dendritic growth[J]. Journal of parallel and distributed computing, 2006, 66(11): 1379-1386.

[12]HAJIHASHEMI M R, El-SHENAWEE M. High performance computing for the level-set reconstruction algorithm[J]. Journal of parallel and distributed computing, 2010, 70(6): 671-679.

本文引用格式：

黃筱云, 董國海, 趙利平，等. 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.

A parallelized fast marching method for reinitialization of level set function

HUANG Xiaoyun1,2,3，DONG Guohai1，ZHAO Liping2，CHENG Yongzhou2

(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 410004, China; 3. State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China)

Abstract：In order to increase computational efficiency of reinitializing level set function, a parallelization strategy of the fast marching method was proposed based on domain decomposition parallelization idea, and the fast parallelized reinitialization of level set function was achieved. Based on domain parallelization idea, a parallelization strategy of the fast marching method was proposed and the fast parallelized reinitialization of level set function was achieved so as to further increase computational efficiency of reinitializing level set function by the fast marching method. The accuracy and computational efficiency of the new parallel algorithm for level set function reinitialization were discussed through some examples of sphere, pentafoil cube and circular cube. It is shown that, compared with the serial fast marching method, the parallel algorithm maintains the accuracy of the serial algorithm of 1st order and remarkably decreases computational time of reinitialization in which the best speedup of the method can approach 5 under the thread number of 8.

Keywords：level set function; reinitialization; fast marching method; parallelization；domain decomposition；parallel algorithm; speedup

收稿日期：2015-02-02.

基金項目：國家自然科學基金青年基金資助項目(51109018, 41176072)；中國博士后科學基金資助項目(2014M561230)；水文水資源與水利工程科學國家重點實驗室開放研究基金資助項目(2013491411).

作者簡介：黃筱云(1980-), 男, 講師, 博士. 通信作者：黃筱云, E-mail:huangxiaoyun@csust.edu.cn.

DOI:10.11990/jheu.201502005

中圖分類號：TV131.2

文獻標志碼：A

文章編號：1006-7043(2016)05-0666-07

網絡出版時間：2016-04-11.