2. 长沙理工大学 水利工程学院, 湖南 长沙 410004;
3. 河海大学 水文水资源与水利工程科学国家重点实验室, 江苏 南京 210098
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
基于欧拉方法模拟波浪传播过程,要求精确捕捉或追踪自由表面位置。与其他界面捕捉或追踪方法不同,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函数重新初始化结果分析该并行化快速步进法的计算精度,探讨多线程并行计算效率。
1 Level set方法概述Level set方法采用的控制方程为
$\frac{{\partial \varphi }}{{\partial t}} + u \cdot \nabla \varphi = 0$ | (1) |
为保证式(1)计算的level set函数具有符号距离函数的特征,即
$\left| {\nabla \varphi } \right| = 1$ | (2) |
需要对式(1)的计算结果进行重新初始化。迭代重新初始化通过计算下面方程稳定解来完成:
$\frac{{\partial \varphi }}{{\partial \tau }} + S\left( \varphi \right)\left[{\left| {\nabla \varphi } \right| - 1} \right] = 0$ | (3) |
式中:S(φ)为光滑函数,τ为虚拟时间;而快速步进法则是运用一种迎风格式对方程(2)的边值问题直接求解。
2 快速步进法的基本原理快速步进方法通常采用一阶迎风格式来对式(2)进行求解和推进,其数值离散格式具体如下所示:
$\left| {\nabla \varphi } \right| = {\left[{\max {{\left( {D_{ijk}^{ - x}\varphi ,- D_{ijk}^{ + x}\varphi ,0} \right)}^2} + \max {{\left( {D_{ijk}^{ - y}\varphi ,- D_{ijk}^{ + y},0} \right)}^2} + \max {{\left( {D_{ijk}^{ - z},D_{ijk}^{ + z}\varphi ,0} \right)}^0}} \right]^{1/2}} = 1$ | (4) |
其中,$D_{ijk}^{ - x}\varphi = \left( {{\varphi _{i,j.k}} - {\varphi _{i - 1,j,k}}} \right)/\Delta h;D_{ijk}^{ + x}\varphi = \left( {{\varphi _{i + 1,j,k}} - {\varphi _{i,j,k}}} \right)/\Delta h;D_{ijk}^{ - y}\varphi = \left( {{\varphi _{i,j,k}} - {\varphi _{i,j - 1,k}}} \right)/\Delta h$;Dijk+yφ=φi+1,j,k-φi,j,k/Δh;Dijk-zφ=φi,j,k-φi,j,k-1/Δh;Dijk+zφ=φi,j,k+1-φi,j,k/Δh。
式中△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节点上的φ值。
3 快速步进法并行策略作为一种主流的并行策略,分区并行是将整个计算区域划分成多个子区域,并分配给不同线程计算;与单区计算相比,分区并行计算需要在边界上交换数据,以使计算结果尽量和单区计算的结果一致,要实现快速步进法的分区并行就要考虑步进过程中相邻区域的影响。并行快速步进算法首先需要将整个计算区域分成多个子区域,然后利用不同的线程来分别执行不同子区域的计算任务,再在子区域边界上交换数据,保证子区域边界处节点选择的迎风性。
3.1 子区域划分要实现节点选择的迎风性,需要在计算过程中保持子区域之间的信息传递。这要求在子区域的划分过程中,为每个子区域的边界外布置虚拟网格,以接收相邻子区域传递来的信息。图 2中,整个计算区域被划分成4个子区域,不同子区域需要被分配不同的线程分别进行计算,每个子区域在与其他子区域相邻的边界外布置虚拟网格以接收信息,如子区域1在右边界和上边界布置的虚拟网格将分别接收子区域2左边界网格和子区域3下边界网格上信息。
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的φ值是快速步进算法的边界条件,所以禁止计算流程改变其节点类型。
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给出分区域并行计算给出的不同φ值等值面,可以看出,子区域间同步算法能够有效地将左半区φ值信息传递给右半区。
表 1给出了不同网格分辨率下并行化快速步进法获得的φ等值面包围的体积及φ值计算误差。从表 1中可以看出,并行化快速步进法计算得到的φ= 0.1等值面包围的体积与精确解相差很小,在200×200×200网格分辨率下,体积误差小于0.1%,且误差阶数为1阶。
分辨率 | 体积 | 损失/% | L1误差 | 阶数 |
真值 | 1.13×10-1 | — | — | — |
50×50×50 | 1.09×10-2 | 3.71 | 8.46×10-3 | — |
100×100×100 | 1.11×10-2 | 1.94 | 4.09×10-3 | 1.05 |
200×200×200 | 1.13×10-2 | 0.09 | 2.11×10-3 | 0.96 |
由于并行化快速步进法获得的L1误差量级为10-3,分区域计算时可设定边界容许误差为1.0×10-4。将圆球放置在计算区域中心,可比较不同线程数下实现φ值重新初始化所需的CPU计算时间,具体计算结果如表 2所示。
s | ||||
分辨率 | 1线程 | 2线程 | 4线程 | 8线程 |
50×50×50 | 0.87 | 0.43 | 0.26 | 0.28 |
100×100×100 | 9.11 | 4.51 | 2.74 | 1.79 |
200×200×200 | 86.09 | 44.08 | 26.10 | 16.56 |
从表 2可以看出,并行化快速步进法能够有效缩短计算时间,在网格分辨率较大的条件下,加速比能达到5左右,但子区域间同步开销使得效率会随着线程数增加而减小,当线程数达到一定数量时,计算时间甚至会略有增加。
4.2 五叶管本算例的五叶管由y-z平面上五叶图形在x方向延长所形成,五叶图形表达式如下:
$\left\{ \begin{gathered} y = R\sin \theta + rsin\left( {5\theta } \right)\sin \theta + 0.5 \hfill \\ z = R\cos \theta + r\sin \left( {5\theta } \right)\cos \theta + 0.5 \hfill \\ \end{gathered} \right.$ | (5) |
在本算例中,R= 0.2,r= 0.1,计算区域大小为2×1×1,网格分辨率200×100×100。图 7为φ= 0和 0.1等值面。
表 3给出了φ= 0.1等值面包围的体积以及φ值计算误差。从中可以看出,一方面,随着网格分辨率的增加,体积损失逐渐减小,当网格分辨率为200×100×100时,体积损失仅为0.22%;另一方面,尽管五叶管表面某些区域曲率较大,但计算结果误差阶数仍大于0.75。
分辨率 | 体积 | 损失/% | L1误差 | 阶数 |
真值 | 7.89×10-1 | — | — | — |
50×25×25 | 7.66×10-1 | 2.85 | 7.87×10-3 | — |
100×50×50 | 7.74×10-1 | 1.88 | 4.64×10-3 | 0.76 |
200×100×100 | 7.87×10-1 | 0.22 | 2.60×10-3 | 0.83 |
在x方向按等长度将整个计算区域划分成多个子区域,分配给不同线程进行计算。表 4给出不同网格分辨率下,多线程并行计算所耗费的CPU计算时间。从中可以看出,子区域间同步开销抵消了双线程带来的运算加速,但随着线程数增加,计算时间进一步缩减,在网格分辨率为200×100×100的条件下,加速比能达到2.5。
分辨率 | 1线程 | 2线程 | 4线程 | 8线程 |
50×50×50 | 0.87 | 0.43 | 0.26 | 0.28 |
100×100×100 | 9.11 | 4.51 | 2.74 | 1.79 |
200×200×200 | 86.09 | 44.08 | 26.10 | 16.56 |
本算例将重新初始化圆环管外φ值分布,计算区域大小为1×1×1,分辨率为100×100×100,为平均分配每个线程计算量,圆环管放置在计算区域中心,其表达式为
$\left\{ \begin{gathered} x = \left( {R + r\cos t} \right)\sin \theta + 0.5 \hfill \\ y = \left( {R + r\cos t} \right)\cos \theta + 0.5 \hfill \\ z = r\sin t + 0.5 \hfill \\ \end{gathered} \right.$ | (6) |
其中,R=0.3,r= 0.05。φ零等值面如图 8所示。图 9(a)给出了φ= 0.1的等值面,图 9(b)~(i)则给出了分区域计算各子区域内φ= 0.1等值面。
表 5则给出了不同网格分辨率下并行化快速步进法获得的φ= 0.1等值面包围的体积及φ值计算误差。通过对比可以看出,φ= 0.1等值面包围体积与真值相差很小,体积损失与网格分辨率大小几乎成线性关系,且并行化快速步进法的计算精度几乎达到1阶。
分辨率 | 体积 | 损失/% | L1误差 | 阶数 |
真值 | 1.33×10-1 | — | — | — |
50×50×50 | 1.28×10-1 | 4.11 | 6.82×10-3 | — |
100×100×100 | 1.30×10-1 | 2.56 | 3.51×10-3 | 0.96 |
200×200×200 | 1.31×10-1 | 1.58 | 1.81×10-3 | 0.92 |
依次在x、z、y方向上平分计算区域,给出不同线程数下实现φ值重新初始化所需的CPU计算时间(表 6)。可以看出,随着网格分辨率变大,子区域间同步开销也会大幅增加,甚至大于快速步进法本身的开销,但随着线程数的增加,加速比能够达到2。此外,计算区域的划分方式也会对加速比造成影响(表 7)。
分辨率 | 1线程 | 2线程 | 4线程 | 8线程 |
50×50×50 | 0.70 | 0.73 | 0.46 | 0.40 |
100×100×100 | 8.39 | 8.13 | 5.15 | 4.01 |
200×200×200 | 78.02 | 149.04 | 71.45 | 53.35 |
通过分区并行,提出了一种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, WVCHNER 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. |