黄筱云, 董国海, 赵利平, 程永舟

(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.