﻿ 快速粒子level set方法的并行化
«上一篇
 文章快速检索 高级检索

 哈尔滨工程大学学报  2018, Vol. 39 Issue (9): 1478-1484  DOI: 10.11990/jheu.201703082 0

### 引用本文

HUANG Xiaoyun, XIA Bo, CHENG Yongzhou, et al. Parallelization of the fast particle level set method[J]. Journal of Harbin Engineering University, 2018, 39(9), 1478-1484. DOI: 10.11990/jheu.201703082.

### 文章历史

1. 长沙理工大学 水利工程学院, 湖南 长沙 410114;
2. 大连理工大学 海岸和近海工程国家重点实验室, 辽宁 大连 116024

Parallelization of the fast particle level set method
HUANG Xiaoyun1,2, XIA Bo1, CHENG Yongzhou1, ZHAO Liping1
1. School of Hydraulic Engineering, Changsha University of Science and Technology, Changsha 410114, China;
2. State Key Laboratory of Coastal and offshore Engineering, Dalian University of Technology, Dalian 116024, China
Abstract: To shorten the computational time of the fast particle level set (FPLS) method for tracking interfaces, the domain and particle decomposition methods were respectively adopted to reduce the overhead of Euler (advection equation of level set function and reinitialization) and Lagrangian (time integration of particle and error correction) computations. The new algorithm was achieved by the OpenMP technique and demonstrated by typical examples, such as Zalesak's disk, single vortex, and 3D transformation. The final results showed that the speedup ratio at four threads can exceed 2 and that the speedup ratio at eight threads can approach 4. These findings reveal that the parallelized FPLS method has good applicability and scalability.
Keywords: fast particle level set method    parallelization    domain decomposition    particle decomposition    OpenMP    Zalesak's disk    single vortex

1 快速粒子level set方法的基本原理

Level set方法是将运动的交界面看作一个高阶符号距离函数的零等值面，其控制方程如下

 ${\varphi _t} + \mathit{\boldsymbol{u}} \cdot \nabla \varphi = 0$ (1)

 $\left| {\nabla \varphi } \right| = 1$ (2)

1.1 半拉格朗日方法

 $\begin{array}{l} \varphi _{i,j}^{n + 1} = \alpha \beta \varphi _{r + 1,s + 1}^n + \left( {1 - \alpha } \right)\beta \varphi _{r,s + 1}^n + \\ \;\;\;\;\;\;\;\;\;\;\alpha \left( {1 - \beta } \right)\varphi _{r + 1,s}^n + \left( {1 - \alpha } \right)\left( {1 - \beta } \right)\varphi _{r,s}^n \end{array}$ (3)

 $r = i - \left[ {{u_{i,j}}\Delta t/\Delta h} \right],\alpha = \left[ {\left( {j - r} \right)\Delta h - {u_{i,j}}\Delta t} \right]/\Delta h$
 $s = j - \left[ {{v_{i,j}}\Delta t/\Delta h} \right],\beta = \left[ {\left( {j - s} \right)\Delta h - {v_{i,j}}\Delta t} \right]/\Delta h$

1.2 快速步进法

 $\begin{array}{l} \left| {\nabla {\varphi _{i,j}}} \right| = \left( {\max {{\left( {D_{ij}^{ - x}\varphi , - \varphi ,0} \right)}^2} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\left. {\max {{\left( {D_{ij}^{ - y}\varphi , - \varphi ,0} \right)}^2}} \right)^{1/2}} = 1 \end{array}$ (4)

 $D_{ij}^{ - x}\varphi = \left( {{\varphi _{i,j}} - {\varphi _{i - 1,j}}} \right)/\Delta h$
 $D_{ij}^{ + x}\varphi = \left( {{\varphi _{i + 1,j}} - {\varphi _{i,j}}} \right)/\Delta h$
 $D_{ij}^{ - y}\varphi = \left( {{\varphi _{i,j}} - {\varphi _{i,j - 1}}} \right)/\Delta h$
 $D_{ij}^{ + y}\varphi = \left( {{\varphi _{i,j + 1}} - {\varphi _{i,j}}} \right)/\Delta h$

1.3 误差校正

 ${\varphi _p}\left( {{\mathit{\boldsymbol{x}}_{i,j}}} \right) = {s_p}\left( {{\mathit{\boldsymbol{r}}_p} - \left| {{\mathit{\boldsymbol{x}}_{i,j}} - {\mathit{\boldsymbol{x}}_p}} \right|} \right)$ (5)

 ${\varphi ^ + } = \mathop {\max }\limits_{\forall p \in {E^ + }} \left( {{\varphi _p},{\varphi ^ + }} \right)$ (6)

 ${\varphi ^ - } = \mathop {\min }\limits_{\forall p \in {E^ + }} \left( {{\varphi _p},{\varphi ^ - }} \right)$ (7)

 $\varphi = \left\{ \begin{array}{l} {\varphi ^ + },\;\;\;\;\left| {{\varphi ^ + }} \right| \le \left| {{\varphi ^ - }} \right|\\ {\varphi ^ - },\;\;\;\;\left| {{\varphi ^ + }} \right| > \left| {{\varphi ^ - }} \right| \end{array} \right.$ (8)

2 快速粒子level set方法的并行 2.1 计算区域分解与粒子分解

2.2 对流方程与重新初始化

1) 确定紧邻交界面网格节点(下面称之为交界面节点)的$φ$值，将其他节点初始化；

2) 最小堆的初始化，即试算与交界面节点相邻节点，并放置到最小堆中；

3) 若最小堆非空，取出$φ$值最小的节点，对其相邻的节点进行更新，并将相邻节点放置在最小堆中；若取出的节点为边界节点，则将该节点的信息传递给相邻区域；否则，前往6)；

4) 根据接收的信息，排查需要重演的节点，并将这些节点重新放回至最小堆中；

5) 返回步骤3)；

6) 等待其他线程的最小堆为空；若接受到新的节点信息，返回4)。

2.3 粒子计算

3 算例及结果讨论 3.1 Zalesak圆盘

 $\left\{ \begin{array}{l} u = \left( {{\rm{ \mathsf{ π} /314}}} \right)\left( {{\rm{50}} - y} \right)\\ v = \left( {{\rm{ \mathsf{ π} /314}}} \right)\left( {x - 50} \right) \end{array} \right.$ (9)
 Download: 图 2 Zalesak圆盘数值结果与真解的比较 Fig. 2 Comparisons of the Zalesak′s disk between the numerical results and the exact solution

 Download: 图 3 不同线程数目下FPLS主要计算环节耗时占比 Fig. 3 Percentage of execution time of the main steps in FPLS method under the different total number of threads
 Download: 图 4 各环节的加速比以及总加速比 Fig. 4 Speedups of the main steps and the total execution
3.2 单涡

 $\left\{ \begin{array}{l} u = - {\sin ^2}\left( {{\rm{ \mathsf{ π} }}x} \right)\sin \left( {2{\rm{ \mathsf{ π} }}y} \right)\cos \left( {{\rm{ \mathsf{ π} }}t/T} \right)\\ v = \sin \left( {{\rm{2 \mathsf{ π} }}x} \right){\sin ^2}\left( {{\rm{ \mathsf{ π} }}y} \right)\cos \left( {{\rm{ \mathsf{ π} }}t/T} \right) \end{array} \right.$ (10)

 Download: 图 5 并行FPLS方法的单涡测试 Fig. 5 Validation of the parallel FPLS method by the single vortex flow

3.3 Zalesak圆球

Zalesak圆球追踪是为了检验FPLS方法在三维条件下的耗散性。本算例中，半径为15的Zalesak球Zalesak圆球布置在尺寸大小为100×100×100的计算区域内，缺口的深度和宽度分别为12.5和5。在式(9)给定的流场作用下，圆球在z＝50的平面上从初始位置(50, 75, 50)开始，绕区域中心作圆周运动。旋转一周同样需628个单位时间。三维算例所采用的最大线程数为8。图 6给出了8线程下FPLS方法的计算结果，最终的体积损失小于1%。不同线程数下各主要计算环节加速比及总加速比见图 4(c)，当线程数增加时，加速比也随之增加，8线程所获得总加速比能达到4.0。

 Download: 图 6 Zalesak圆球在不同时刻的位置 Fig. 6 Snapshots of time evolution of a Zalesak′s sphere at different time
3.4 三维变形

 $\left\{ \begin{array}{l} u = 2{\sin ^2}\left( {{\rm{ \mathsf{ π} }}x} \right)\sin \left( {{\rm{2 \mathsf{ π} }}y} \right)\sin \left( {{\rm{2 \mathsf{ π} }}z} \right)\cos \left( {{\rm{ \mathsf{ π} }}t/T} \right)\\ v = - \sin \left( {{\rm{2 \mathsf{ π} }}x} \right){\sin ^2}\left( {{\rm{ \mathsf{ π} }}y} \right)\sin \left( {{\rm{2 \mathsf{ π} }}z} \right)\cos \left( {{\rm{ \mathsf{ π} }}t/T} \right)\\ w = - \sin \left( {{\rm{2 \mathsf{ π} }}x} \right)\sin \left( {{\rm{ \mathsf{ π} }}y} \right){\sin ^2}\left( {{\rm{2 \mathsf{ π} }}z} \right)\cos \left( {{\rm{ \mathsf{ π} }}t/T} \right) \end{array} \right.$ (10)

 Download: 图 7 圆球在涡旋流场中的变形 Fig. 7 Snapshots of the time evolution of a sphere in the vortex flow
4 结论

1) 并行的粒子Level set方法保留了原有串行算法的精度和守恒性，而计算时间大大缩减，在二维算例中，4线程能够达到的加速比超过2，而在三维算例中，8线程的加速比也接近4。

2) 粒子Level set方法中，粒子时程积分和误差校正的计算开销比重占总计算开销的7成以上，将粒子计算开销在线程间平均，能够显著提高计算效率。

3) 快速粒子Level set方法的对流计算和重新初始化只在界面附近区域进行，区域分解无法保证这两步计算在线程间保持荷载平衡。但由于其计算开销比重较低，对整体并行计算效率影响较小。

 [1] XIE Zhihua. Numerical modelling of wind effects on breaking solitary waves[J]. European journal of mechanics-B/fluids, 2014, 43: 135-147. DOI:10.1016/j.euromechflu.2013.08.001 (0) [2] LUPIERI G, CONTENTO G. Numerical simulations of 2-D steady and unsteady breaking waves[J]. Ocean engineering, 2015, 106: 298-316. DOI:10.1016/j.oceaneng.2015.07.014 (0) [3] GAETA M G, LAMBERTI A. The role of air modeling on the numerical investigation of coastal dynamics and wave-structure interactions[J]. Computers & fluids, 2015, 111: 114-126. (0) [4] HUANG Xiaoyun, LI Shaowu. A two-dimensional numerical wave flume based on SA-MPLS method[J]. Acta oceanologica sinica, 2012, 31(3): 18-30. DOI:10.1007/s13131-012-0203-2 (0) [5] LI Shaowu, ZHUANG Qian, HUANG Xiaoyun, et al. 3D simulation of flow with free surface based on adaptive octree mesh system[J]. Transactions of Tianjin university, 2015, 21(1): 32-40. DOI:10.1007/s12209-015-2314-2 (0) [6] MCSHERRY R J, CHUA K V, STOESSER T. Large eddy simulation of free-surface flows[J]. Journal of hydrodynamics, 2017, 29(1): 1-12. DOI:10.1016/S1001-6058(16)60712-6 (0) [7] ENRIGHT D, FEDKIW R, FERZIGER J, et al. A hybrid particle level set method for improved interface capturing[J]. Journal of computational physics, 2002, 183(1): 83-116. DOI:10.1006/jcph.2002.7166 (0) [8] 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. DOI:10.1006/jcph.1994.1155 (0) [9] SUSSMAN M, PUCKETT E G. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[J]. Journal of computational physics, 2000, 162(2): 301-337. DOI:10.1006/jcph.2000.6537 (0) [10] VAN DER PIJL S P, SEGAL A, VUIK C, et al. A mass-conserving level-set method for modelling of multi-phase flows[J]. International journal for numerical methods in fluids, 2005, 47(4): 339-361. DOI:10.1002/(ISSN)1097-0363 (0) [11] LV Xin, ZOU Qingping, ZHAO Yong, et al. A novel coupled level set and volume of fluid method for sharp interface capturing on 3D tetrahedral grids[J]. Journal of computational physics, 2010, 229(7): 2573-2604. DOI:10.1016/j.jcp.2009.12.005 (0) [12] WANG Y, SIMAKHINA S, SUSSMAN M. A hybrid level set-volume constraint method for incompressible two-phase flow[J]. Journal of computational physics, 2012, 231(19): 6438-6471. DOI:10.1016/j.jcp.2012.06.014 (0) [13] ENRIGHT D, LOSASSO F, FEDKIW R. A fast and accurate semi-Lagrangian particle level set method[J]. Computers & structures, 2005, 83(6/7): 479-490. (0) [14] HIEBER S E, KOUMOUTSAKOS P. A Lagrangian particle level set method[J]. Journal of computational physics, 2005, 210(1): 342-367. DOI:10.1016/j.jcp.2005.04.013 (0) [15] WANG Zhaoyuan, YANG Jianming, STERN F. An improved particle correction procedure for the particle level set method[J]. Journal of computational physics, 2009, 228(16): 5819-5837. DOI:10.1016/j.jcp.2009.04.045 (0) [16] IANNIELLO S, MASCIO A D. A self-adaptive oriented particles level-set method for tracking interfaces[J]. Journal of computational physics, 2010, 229(4): 1353-1380. DOI:10.1016/j.jcp.2009.10.034 (0) [17] VARTDAL M, BØCKMAN A. An oriented particle level set method based on surface coordinates[J]. Journal of computational physics, 2013, 251: 237-250. DOI:10.1016/j.jcp.2013.05.044 (0) [18] 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. DOI:10.1016/j.jcp.2015.06.039 (0) [19] 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. DOI:10.1016/j.jpdc.2006.02.005 (0) [20] MIRZADEH M, GUITTET A, BURSTEDDE C, et al. Parallel level-set methods on adaptive tree-based grids[J]. Journal of computational physics, 2016, 322: 345-364. DOI:10.1016/j.jcp.2016.06.017 (0) [21] CUMINATO J, FILHO A C, BOAVENTURA M, et al. Simulation of free surface flows in a distributed memory environment[J]. Journal of computational and applied mathematics, 1999, 103(1): 77-92. DOI:10.1016/S0377-0427(98)00242-8 (0) [22] QIANG Ji, FURMAN M A, RYNE R D. A parallel particle-in-cell model for beam-beam interaction in high energy ring colliders[J]. Journal of computational physics, 2004, 198(1): 278-294. DOI:10.1016/j.jcp.2004.01.008 (0) [23] QIANG Ji, LI Xiaoye. Particle-field decomposition and domain decomposition in parallel particle-in-cell beam dynamics simulation[J]. Computer physics communications, 2010, 181(12): 2024-2034. DOI:10.1016/j.cpc.2010.08.021 (0) [24] 黄筱云, 董国海, 赵利平, 等. 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. (0) [25] 黄筱云, 董国海, 常佳夫, 等. Level set函数快速步进重构并行算法的改进[J]. 哈尔滨工程大学学报, 2017, 38(6): 836-842. HUANG Xiaoyun, DONG Guohai, CHANG Jiafu, et al. Improvement of parallel fast marching method for reconstruction of level set function[J]. Journal of Harbin Engineering University, 2017, 38(6): 836-842. (0)