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

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

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

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

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

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

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

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

關(guān)鍵詞：level set函數(shù)；重新初始化；快速步進法；并行；分區(qū)；并行算法；加速比

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/23.1390.U.20160411.1031.038.html

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

1Level set方法概述

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

(1)

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

(2)

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

(3)

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

2快速步進法的基本原理

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

(4)

式中△h為單元網(wǎng)格尺寸?？梢钥闯?，該格式對網(wǎng)格節(jié)點的選擇類似于根據(jù)流的信息域選擇節(jié)點，這種節(jié)點選擇方式可以保證信息始終向一個方向傳遞，即從較小值的節(jié)點向較大值的節(jié)點傳遞。當(dāng)獲得交界面周圍網(wǎng)格節(jié)點的φ值后，就可以向外構(gòu)造出φ> 0區(qū)域內(nèi)網(wǎng)格節(jié)點上的φ值，若改變區(qū)域內(nèi)網(wǎng)格節(jié)點上φ值的符號，就可構(gòu)造出φ< 0區(qū)域內(nèi)網(wǎng)格節(jié)點上的φ值。快速步進法的具體實現(xiàn)步驟可以歸納如下：

1)標(biāo)記交界面周圍的φ值已知網(wǎng)格節(jié)點為Alive節(jié)點，與Alive節(jié)點相鄰其他節(jié)點為Close節(jié)點，剩余節(jié)點為Far節(jié)點；

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

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

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

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

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

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

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

3快速步進法并行策略

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

3.1子區(qū)域劃分

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

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

3.2具體并行算法

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

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

2)進行子區(qū)域同步；

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

在并行計算過程中，子區(qū)域同步可按下面步驟執(zhí)行：

①子區(qū)域虛擬網(wǎng)格節(jié)點接收相鄰區(qū)域邊界相應(yīng)節(jié)點上的值；

②比較子區(qū)域邊界節(jié)點與相應(yīng)虛擬節(jié)點上的φ值，若子區(qū)域邊界節(jié)點的φ值大于相應(yīng)虛擬節(jié)點，表明該節(jié)點會受到相鄰區(qū)域的影響，重新標(biāo)記該節(jié)點為Close節(jié)點，并將虛擬節(jié)點的φ值賦給該節(jié)點；

③找出所有受相鄰區(qū)域影響邊界節(jié)點中φ值最小者，標(biāo)記為Alive節(jié)點，并保存其φ值為φmin；

④重新標(biāo)記子區(qū)域內(nèi)網(wǎng)格節(jié)點，若節(jié)點φ>φmin，則該節(jié)點標(biāo)記為Far節(jié)點。

在圖3中，相鄰區(qū)域右邊界節(jié)點的φ值先賦給主區(qū)域虛擬網(wǎng)格節(jié)點；然后比較主區(qū)域邊界節(jié)點和虛擬節(jié)點的φ值，標(biāo)記受影響邊界節(jié)點為Close節(jié)點，并更新其φ值；再找到邊界上φ值最小的受影響節(jié)點A，并以節(jié)點A的φ值為閾值，重新標(biāo)記主區(qū)域內(nèi)所有網(wǎng)格節(jié)點。

在同步過程中，如果所有子區(qū)域受相鄰區(qū)域影響邊界節(jié)點集均為空，則表明整個區(qū)域內(nèi)網(wǎng)格節(jié)點φ值符合迎風(fēng)特征，計算即可結(jié)束。

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

(a) 虛擬網(wǎng)格節(jié)點賦值

(b) 標(biāo)記受影響的邊界節(jié)點

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

(d) 重新標(biāo)記主區(qū)域內(nèi)網(wǎng)格節(jié)點圖3　子區(qū)域同步過程Fig.3　Synchronization of sub-domain

圖4　節(jié)點B禁止重新標(biāo)記Fig.4　Reflag of B node is forbidden

4算例及分析

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

4.1圓球

本算例設(shè)定計算區(qū)域大小為1×1×1，分辨率為100×100×100，整個計算區(qū)域?qū)澐殖勺笥覂刹糠?。圓球位于左側(cè)計算區(qū)域，如圖5所示，其球心位置為[0.3,0.5,0.5]，球體半徑為0.2。

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

給出的不同φ值等值面，可以看出，子區(qū)域間同步算法能夠有效地將左半?yún)^(qū)φ值信息傳遞給右半?yún)^(qū)。

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

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

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

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

表1不同網(wǎng)格分辨率下圓球計算誤差

Table 1Computational error of sphere under differentgrid resolution

分辨率體積損失/%L1誤差階數(shù)真值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可以看出，并行化快速步進法能夠有效縮短計算時間，在網(wǎng)格分辨率較大的條件下，加速比能達到5左右，但子區(qū)域間同步開銷使得效率會隨著線程數(shù)增加而減小，當(dāng)線程數(shù)達到一定數(shù)量時，計算時間甚至?xí)杂性黾印?/p>

表2不同線程數(shù)計算耗費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，計算區(qū)域大小為2×1×1，網(wǎng)格分辨率200×100×100。圖7為φ= 0和 0.1等值面。

(a) 虛擬網(wǎng)格節(jié)點賦值

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

分辨率體積損失/%L1誤差階數(shù)真值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等值面包圍的體積以及φ值計算誤差。從中可以看出，一方面，隨著網(wǎng)格分辨率的增加，體積損失逐漸減小，當(dāng)網(wǎng)格分辨率為200×100×100時，體積損失僅為0.22%；另一方面，盡管五葉管表面某些區(qū)域曲率較大，但計算結(jié)果誤差階數(shù)仍大于0.75。

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

表4不同線程數(shù)計算耗費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圓環(huán)管

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

(6)

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

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

分辨率體積損失/%L1誤差階數(shù)真值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則給出了不同網(wǎng)格分辨率下并行化快速步進法獲得的φ= 0.1等值面包圍的體積及φ值計算誤差。通過對比可以看出，φ= 0.1等值面包圍體積與真值相差很小，體積損失與網(wǎng)格分辨率大小幾乎成線性關(guān)系，且并行化快速步進法的計算精度幾乎達到1階。

表6不同線程數(shù)計算耗費CPU時間(圓環(huán)管)

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方向上平分計算區(qū)域，給出不同線程數(shù)下實現(xiàn)φ值重新初始化所需的CPU計算時間(表6)?？梢钥闯?，隨著網(wǎng)格分辨率變大，子區(qū)域間同步開銷也會大幅增加，甚至大于快速步進法本身的開銷，但隨著線程數(shù)的增加，加速比能夠達到2。此外，計算區(qū)域的劃分方式也會對加速比造成影響(表7)。

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

Table 7Comparison of parallelized computation time under different partition strategiess

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

5結(jié)論

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

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

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

3)并行化快速步進法的計算精度會受到表面邊界曲率條件的影響，邊界越光滑，獲得的計算結(jié)果精度越接近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函數(shù)重新初始化的并行快速步進法[J]. 哈爾濱工程大學(xué)學(xué)報， 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.

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

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

DOI:10.11990/jheu.201502005

中圖分類號：TV131.2

文獻標(biāo)志碼：A

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

網(wǎng)絡(luò)出版時間：2016-04-11.