Level set函数快速步进重构的隐式分区并行

黄筱云，夏波，程永舟，江诗群

(1.长沙理工大学 水利工程学院，湖南 长沙 410114； 2.长沙理工大学 水沙科学与水灾害防治湖南省重点实验室，湖南 长沙 116024)

Level set方法[1]、VOF方法[2]、波前追踪法(front-tracking method)[3]和相场方法(phase-field method)[4]是追踪自由表面的4种常用方法。Level set方法的最大优势是能够容易处理自由表面的拓扑变化，但在相同网格分辨率下其计算精度不如显式的波前追踪法，其质量守恒性弱于VOF法[5]。为弥补level set方法的上述不足，主要的措施包括采用高精度的离散格式(如空间上5阶WENO格式、时间上3阶TVD Runge-Kutta格式)[6-8]、重构level set方程的解[9-11]等。重构将强制level set方程的解符合符号距离函数条件。为实现重新初始化，一种方式是迭代求解伪时间瞬态重构方程[9-12]，另一种是直接求解Eikonal方程。现阶段，Eikonal方程的快速求解方法有快速步进法(fast marching method)[13-15]、快速扫描法(fast sweeping method)[16-18]2类。对于几何特征简单的交界面，快速扫描法的计算速度要快于快速步进法[19]。不过，快速扫描法需要对计算全域进行遍历，而快速步进法只需遍历自由表面附近的区域。为进一步提高快速步进重构的效率，Herrmann[20]提出了基于区域分解的快速步进并行算法。Tugurlan[21]和Yang等[22]又分别改进了信息传递方式以提高分布式内存架构下并行计算的效率。黄筱云等[23-24]则提出了基于共享内存架构的快速步进重构算法，解决了相互重叠内边界节点重构值传递问题。在此基础上，夏波等[25]提出了一种自适应分区策略，确保各分区计算荷载平衡，进一步提高并行重构效率。

在level set方法中，重新初始化会定期实施，甚至是每一步都需进行level set函数的重构。若采用文献[25]的自适应分区技术，区域分解将会带来额外的计算开销。因此，本文将充分利用共享内存架构的优势，提出一种隐式分区并行策略，避免区域分解带来额外开销。

1 显式分区并行重构

1.1 快速步进分区并行重构

Level set函数的重构(又称为重新初始化)被认为是求解如下问题：

(1)

式中：Ω表示整个计算区域；Γ为交界面。快速步进重构采用1阶迎风格式离散方程(1)，即节点(i,j)上φi,j的空间导数由该节点与φ<φi,j的相邻节点所决定[26]。离散后的方程(1)则变成一个关于φi,j的二次方程。若给定那些φ<φi,j相邻节点的值，即可求出φi,j的值。因此，选择从交界面开始，向外依次求解每个节点上的二次方程，便能求出整个域内所有节点的值。具体操作时，先采用文献[27]中的方法重构交界面周围网格节点，将其作为起始边界；再试算起始边界外侧所有网格节点的值，并将这些节点根据其值放置到最小堆中，其中，处于根节点的网格节点被认为满足方程(1)，或者说重构成功；然后，将根节点上的节点取出，试算其相邻未重构的节点，并将这些节点放置到最小堆中，或调整它们在最小堆中的位置；最后，将根节点的子孙节点中的网格节点依次向上递补。

基于区域分解的并行重构预先将整个计算区域划分成p个子区域，并分配给对应的线程，每个线程独立完成各自子区域内的重构。这时，需要求解的问题形式为：

(2)

在文献[23-25]中，计算区域将沿网格线进行分割，故子区域边界节点将相互重叠。每个子区域独自实施重构无法确保这些重叠的内边界节点的重构值一致，故需要采取同步措施。在显式分区并行算法中，当重构到边界节点时，该节点重构值将与相邻子区域重叠节点比较，若其值小于重叠节点，则将该值传递给相邻子区域。每个子区域在接收到传递值后，除了更新重叠节点的值外，还将其作为基准，回滚域内重构过的网格节点。当所有子区域内所有节点都重构完毕后，整个区域的重构过程才算结束。具体过程见文献[24]。

1.2 自适应分区策略

Level set方法的重新初始化只需在交界面周围一定范围内执行，传统均分计算区域的策略可能导致子区域间重构计算开销不一致，如图1(a)所示。自适应分区策略则通过平分交界面长度(或面积)来划分计算区域，保证子区域间开销基本一致，如图1(b)所示。

图1 2种显式分区策略的比较Fig.1 Comparison between two explicit domain decompositions

自适应分区是一种递归剖分过程。即先将计算区域分成2个子区域，再分别将子区域分成2个更小的子区域，然后继续将这些更小的子区域分成2块，直至划分子区域的数量与线程数量相同。此外，文献[25]中的自适应分区策略还根据沿不同坐标轴方向划分所产生需要重构内边界节点的数量来确定子区域一分为二的操作。这是因为需要重构的内边界节点数量越少，子区域间信息传递次数以及域内节点回滚的次数也越少。

2 隐式分区并行重构

上面所述的并行算法需提前将整个计算区域沿网格线显式地划分成多个子区域。为保存重叠节点在不同线程下的计算值，实现相互重叠节点间的信息比较和传递，需要为这些节点分配额外的存储空间。当频繁地实施自适应分区并行重构时，存储空间的分配和回收操作会造成不小的计算开销。

本文提出的隐式分区并行重构不直接划分计算区域，而是划分交界面，将问题(1)将变成：

(3)

这意味着将一个以交界面为起始边界的重构问题变成以交界面一部分为起始边界的重构问题的集合。每个线程将独立求解一个交界面部分的重构问题。图1所描述的椭圆形交界面边界重构问题若按照隐式分区并行策略处理，将变成图2所示的情况。即将交界面分成长度一致的4个部分，每个线程完成以交界面的一部分为起始边界的重构。

图2 图1重构的隐式分区并行策略Fig.2 Implicit domain decomposition method for the problem shown in Fig.1

在具体的操作时，无需为每个线程复制整个计算区域，而是允许每个线程均能获取所有节点。在线程更新网格节点值之前，应预先比较写入值与节点既有值的大小，若写入值更小，则更新该节点，当多个线程更新同一节点出现冲突时，可通过OpenMP中的原子操作解决；若节点值更小，写入值将不被采纳，而且线程将不能通过该节点向外继续求解其他未重构的节点。此时，线程自动将该节点作为其子区域的边界节点。

与显式分区并行重构比较，隐式分区并行重构没有子区域间节点值传递环节和回滚环节。隐式分区并行重构方法的特点是允许线程重构其他线程重构过的节点，其中，节点被重新重构称为重演。

在快速步进重构法中，节点根据重构状态分为3类：Far表示未重构、Close表示重构中、Alive表示已重构。而在新的隐式分区并行重构算法中，节点状态将被设定为线程私有。这样，被一个线程重构过的节点可以被另一个线程继续重构。当一个线程重构结束时，所有节点会被标记成Alive或保留Far标记。节点保留Far标记是因为该节点在其他线程下的重构值更小。当重构结束时，重构的节点势必属于某个线程，即以某个线程的计算结果作为其重构值，而属于某个线程重构节点的集合构成该线程控制的区域，这样，整个重构区域分解成多个子区域。

为了保证线程间重构开销平衡，减小节点重复重构的次数，隐式分区并行重构同样采用自适应策略划分交界面。定义交界面为Γ，线程总数为NP，具体步骤如下：

1) 令o=1,n=Np,Γo=Γ。

2) 确定交界面Γo的坐标中值(xc,yc)。

3) 确定交界面Γ0分别与x=xc及y=yc相交的交点个数mx和my。

4) 当mxxc部分移到Γ0+n/2，其余部分保留在Γ0中；否则，将交界面中y>yc的部分移给Γ0+n/2。

5) 若n>2，则o←o,n←n/2,Γo←Γo和o←o+n/2,n←n/2,Γo←Γo+n/2，

6) 返回步骤2)执行递归操作。

3 算例及讨论

3.1 圆球

圆球位于尺寸为100×100×100的计算区域内，其球心位置为(40,60,60)，半径为25。重构区域的范围设定为φ∈[-12,12]。图3为8线程下重构等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10与平面x=40以及平面z=60的交线。图4给出了y=50平面上不同线程控制的分区，可以看出，分区大小基本相等，且其边界垂直于交界面边界，这是因为信息沿交界面法向方向传播。表1给出φ=10等值面所包围的体积以及重构的计算误差，其中，体积计算和误差估计采用文献[25]中的方法，整体计算精度为1阶。图5比较了隐式分区和显式分区下不同线程数量带来的加速比。8线程下显式分区并行的加速比达到5，而隐式分区并行的加速比接近4。虽然相同线程数量下隐式分区并行重构加速比略小于显式分区计算，但前者实际计算时间要小于后者。

图5 隐式与显式分区下圆球并行重构加速比与耗时比较Fig.5 Comparison of speedup & runtime between the implicit and explicit domain decomposition methods in parallel reconstruction of the sphere

表1 圆球等值面φ=10并行重构的计算损失和误差Table 1 Computational loss and error of the parallel reconstructed sphere isosurface of φ=10

图4 平面y=50上不同线程控制的区域(圆球)Fig.4 Regions governed by different threads at the plane of y=50(sphere)

图3 8线程下不同等值面的重构结果与平面x=40和平面z=60的交线 (圆球)Fig.3 Intersections between sphere′s reconstructed isosurfaces and two planes of x=40 & z=60 under 8 threads (sphere)

实际上，隐式分区并行重构中的节点重演和显式分区并行算法中的节点回滚都是为了保证线程间计算的同步。在显式分区并行算法中，定义了回滚参数：

(4)

式中：Nt节点回滚操作总次数;M计算区域内重构节点总数。

本文同样定义一个节点重演参数，也用Fr表示，即:

(5)

式中Mr节点重演操作总次数。

图6给出了采用隐式和显式分区在不同线程数量下的Fr值。可以看出，采用隐式分区算法，节点重演的比例明显低于显式分区算法。

图6 隐式与显式分区下圆球并行重构Fr值的对比Fig.6 Comparison of the Fr value between the implicit and explicit domain decomposition methods in parallel reconstruction of the sphere

3.2 Zalesak球

同样在100×100×100计算区域的(40,60,60)位置上，放置半径为25、缺口宽度为10，深度为25, 缺口在y平面上的投影与x轴夹角呈-45°的Zalesak球，如图7所示。本算例的整个重构范围仍取φ∈[-12,12]。图7给出了8线程下重构等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10与平面y=60以及平面z=60的交线。

图7 8线程下不同等值面的重构结果与平面y=60和平面z=60的交线(Zalesak圆球)Fig.7 Intersections between Zalesak sphere′s reconstructed isosurfaces and two planes of y=60 & z=60 under 8 threads (Zalesak sphere)

不同线程在y=50平面上的控制分区如图8所示，可以看出，与圆球算例同一平面上线程控制区域分布(图4)相比，本算例中不同线程控制区域分布存在差异。φ=5等值面所包围的体积以及重构值的计算误差见表2。

图8 平面y=50上不同线程控制的区域(Zalesak球)Fig.8 Regions governed by different threads at the plane of y=50(Zalesak sphere)

表2 Zalesak球等值面φ=10并行重构的计算损失和误差Table 2 Computational loss and error of parallel reconstruction of the Zalesak′s sphere isosurface of φ=5

Zalesak球隐式与显式分区并行重构的加速比和计算时间如图9所示，容易看出，2种并行算法的加速比相近，8线程下的加速比均接近4，但隐式分区并行重构算法的计算时间明显要少于显式分区，几乎只有前者耗时的1/2。从图10中同样可以看出，在Zalesak球等值面的并行重构中，隐式分区的节点重演次数也少于显式算法的节点回滚次数，前者的Fr值不到后者的1/3。

图10 隐式与显式分区下Zalesak球并行重构Fr值的对比Fig.10 Comparison of the Fr value between the implicit and explicit domain decomposition methods in parallel reconstruction of the Zalesak sphere

图9 隐式与显式分区下Zalesak球并行重构加速比与耗时比较Fig.9 Comparison of speedup & runtime between the implicit and explicit domain decomposition methods in parallel reconstruction of the Zalesak sphere

3.3 哑铃

放置在100×100×100计算区域内的哑铃由半径相同的2个球体和1个圆柱体构成。2个球体分别位于(40,60,60)和(70,30,40)处，半径为20，连接2个球心的圆柱的横截面半径为10。哑铃等值面同样采用8个线程重构，其重构范围依然限制在φ∈[-12,12]。图11给出8线程下哑铃重构等值面f=-10, -8, -6, -4, -2, 2, 4, 6, 8,10与平面y=45以及平面z=50的交线。图12显示平面y=50上不同线程控制的计算区域。

图11 8线程下不同等值面的重构结果与与平面y=45和平面z=50的交线(哑铃)Fig.11 Intersections between dumbbell′s reconstructed isosurfaces and two planes of y=45 & z=50 under 8 threads (dumbbell)

图12 平面y=50上不同线程控制的区域(哑铃)Fig.12 Regions governed by different threads at the plane of y=50(dumbbell)

表3给出了φ=5等值面所包围体积以及重构的计算误差。从图13可以看出，相同线程数量下隐式分区并行计算的加速比略大于显式分区，且隐式分区并行计算耗时几乎只有显式分区的1/2。如图14所示，在哑铃等值面的并行重构中，隐式分区的节点重演次数也少于显式算法的节点回滚次数。

表3 哑铃等值面φ=5并行重构的计算损失和误差Table 3 Computational loss and error of the parallel isosurface reconstruction of the dumbbell isosurface of φ=5

图13 隐式与显式分区下哑铃并行重构加速比与耗时比较Fig.13 Comparison of speedup & runtime between the implicit and explicit domain decomposition methods in parallel reconstruction of the dumbbell

图14 隐式与显式分区下哑铃并行重构Fr值的对比Fig.14 Comparison of the Fr value between the implicit and explicit domain decomposition methods in parallel reconstruction of the dumbbell

4 结论

1)基于共享内存架构的隐式分区并行算法能够有效节约计算时间。

2)相同线程数量下，隐式分区并行算法的加速比与显式分区并行算法基本一致，8线程下并行算法的加速比接近4。

3)在相同线程数量下，隐式分区并行算法的实际耗时要小于显式分区并行算法。

4)隐式分区并行重构的节点重演次数要少于显式分区并行重构的节点回滚次数。