CN102902514A - 半隐式类粒子法的大规模并行处理方法 - Google Patents

Abstract

本发明提供一种半隐式类粒子法的大规模并行处理方法：1)采用高效且易并行的对称兰乔斯算法（Symmetric Lanczos Algorithm,SLA）求解半隐式粒子法中压力泊松方程；2)在服务器上可采用共享内存的并行模型（如OpenMP、OpenCL）和/或消息传递的并行模型（MPI）以及它们的混合模型并行加速算法，本发明可以大规模并行加速需要求解压力泊松方程的半隐式类粒子法，并行加速效果显著。

Description

半隐式类粒子法的大规模并行处理方法

技术领域

本发明属于流体流动计算机仿真领域，涉及对不可压缩流体的自由表面流动的模拟，具体涉及一种半隐式类粒子法的大规模并行处理方法。

背景技术

自由表面流动是一种常见的流动，在溃坝流动、河流的明渠流动，以及蒸发沸腾过程的传热传质过程中都涉及到自由表面流动，采用计算机准确的模拟自由表面流动对工程应用有重大意义；粒子类方法是在拉格朗日框架下追踪粒子的运动，因此易于仿真自由表面流动。

随着计算机技术的发展，采用计算机模拟流体流动已经成为研究流体运动的一种重要的方法。网格类方法是最早采用的一类计算机模拟流体运动的方法，在处理自由表面流动的时候需要引入专门的处理界面运动的方法，如VOF方法和Level Set方法。对于剧烈波动的自由表面流动问题，大变形自由表面的处理是很复杂且很耗时的工作，并且有时会导致求解的错误；另外，网格类方法中对流项的数值耗散也是一个很严重的问题。

粒子类方法是近年来快速发展的一种算法，它们采用的是在拉格朗日框架下追踪粒子的运动，可以追踪自由表面的位置，因此易于仿真自由表面流动。光滑粒子动力学（SPH）和移动粒子半隐式算法（MPS）是两种常见的粒子类方法。前者采用显式算法求解压力，主要用于求解固体力学问题以及可压缩的流体力学问题；而后者采用半隐式算法通过求解压力泊松方程获得压力，主要用于不可压缩性流动的模拟。

通常来说，粒子类方法的计算量是大于网格类方法的计算量，同时工程应用领域也要求大规模粒子数的模拟，因此发展高效大规模并行的粒子类方法具有重大的现实意义。显式的粒子类方法（如SPH）已经得到了较好的并行。本发明主要是针对半隐式类粒子法（如MPS）的并行化。

发明内容

粒子类方法的算法可以分为显式算法和半隐式算法。在显式类算法（如SPH算法）中，粒子的压力是通过状态方程由该粒子所代表的流体微团的密度计算得到，每一个粒子压力的计算是独立于其他粒子的，因此显式类粒子方法是很容易并行的；而在半隐式类粒子算法（如MPS）中，压力是通过全局求解压力泊松方程而获得，因此每一个粒子压力的计算都是依赖于其他粒子的信息的，并且全局求解压力泊松方程是半隐式类粒子法中最为耗时的部分，因此半隐式类粒子法并行的关键在于找到一个高效且易于并行的用于求解压力泊松方程的求解器。

本发明的目的在于提供一种半隐式类粒子法的大规模并行处理方法，该方法先给出一个高效且易于并行的求解器，然后通过并行实现大规模粒子数的模拟。

为达到上述目的，本发明采用了以下技术方案：

该处理方法在并行服务器上采用并行模型并行加速半隐式类粒子法，所述半隐式类粒子法是指黏性项和外力源项显式计算，而压力项则通过求解压力泊松方程隐式计算的粒子类算法，半隐式类粒子法中的压力泊松方程采用对称兰乔斯算法求解。

所述半隐式类粒子法包括移动粒子半隐式算法（Moving Particle Semi-implicitmethod，MPS）和不可压缩的光滑粒子动力学方法（Incompressible Smoothed ParticleHydrodynamics，ISPH）。

所述并行服务器为对称多处理器、GPU服务器、由多CPU核的服务器组成的多节点的服务器、由多GPU核的服务器组成的多节点的服务器或由多CPU核和多GPU核服务器组成的混合多节点服务器。

所述并行模型为共享内存的并行模型、消息传递的并行模型或二者的混合模型。

所述并行加速具体包括以下几种优选情况：

a）如果并行服务器是对称多处理器，采用OpenMP模型（http://www.openmp.org）、MPI模型（张武生，薛巍，李建江，郑伟民《MPI并行程序设计实例教程》清华大学出版社，北京，2009）或采用OpenMP和MPI的混合模型来并行加速；OpenMP模型是基于共享内存的并行模型，由开源网站维护，网站上给出了具体的编程方法；MPI模型是基于消息传递的并行模型；

b）如果并行服务器是GPU服务器，通过OpenCL模型（http://www.khronos.org/opencl/）或CUDA模型（http://developer.nvidia.com/category/zone/cuda-zone）并行加速；OpenCL模型适用于GPU的开源并行模型，由开源网站维护；CUDA是由Nvidia公司开发使用于GPU的并行模型，在公司官网上给出了具体的编程方法；

c）如果并行服务器是多节点的服务器，每个节点是多CPU核的服务器，在每一个节点上，采用OpenMP模型并行加速，而在节点之间采用MPI模型来并行加速；

d）如果并行服务器是多节点的服务器，每个节点是多核的GPU服务器，在每一个节点上，采用OpenCL模型或CUDA模型并行加速，而在节点之间采用MPI模型来并行加速；如果节点中既有多核的CPU节点，又有多核GPU节点，对CPU节点采用OpenMP模型并行加速，对GPU节点采用OpenCL模型或CUDA模型并行加速，在节点之间仍采用MPI模型并行加速。

所述处理方法的具体步骤如下：

1）设置计算区域内粒子的初始分布，将计算区域划分为多个子区域，将每一个子区域分配给一个进程负责计算；

2）将进程上靠近划分边界的粒子信息发送给邻居进程，并且邻居进程上要完成相应的接受操作；

3）在每一个进程上显式计算外力源项和黏性项，然后计算临时速度场和临时位置场，从而计算压力泊松方程的源项；步骤3）的操作是针对子区域上的每一个粒子在进程内并行加速；

4）采用并行对称兰乔斯算法全局求解压力泊松方程，获得压力值；

5）将获得的压力值发送到邻居进程上；

6）传送压力值的通信完成之后，并行计算压力梯度项，用于修正临时速度场，从而获得最终的速度场；

7）计算过程中要动态调节区域划分边界，以保证每个进程上的负载都近似相同；

8）判断本进程上由于粒子移动产生的粒子越界的情况，将已经越过本进程边界的粒子信息发送给相应的进程；

9）判断是否已达到模拟设定时间，如果已达到模拟设定时间则退出；如果没有达到模拟设定时间，返回步骤2），继续下一个时间步长的计算。

粒子法的压力泊松方程离散后通常得到一个主对角占优的大型对称稀疏矩阵，不完全乔莱斯基分解算法（ICCG）是求解大型对称稀疏矩阵的理想迭代方法，但是难以并行。为了实现粒子法的大规模并行运算，需要引入一种高效且易于并行的求解线性方程组的方法。通过多次尝试，本发明最终采用对称兰乔斯算法（SLA）求解压力泊松方程，这种算法易于并行，并且其求解效率优于ICCG算法，并行加速效果显著。

附图说明

图1为背景网格的划分以及中心粒子的作用范围。

图2为子区域划分。

图3为不同子区域（进程）间相互通信。

图4为溃坝问题的模型。

图5为在求解溃坝问题中压力泊松方程时SLA算法和ICCG算法单次求解执行时间的比较。

图6为SLA算法和ICCG算法中易于并行的操作所占的百分比。

图7为OpenMP模型在不同CPU核数下的并行效率。

图8为MPI模型在不同节点数下的并行效率。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

SLA算法求解压力泊松方程的步骤如下：

SLA算法实际上是一种Krylov子空间方法，设m是Krylov子空间的维度，对称兰乔斯算法求解对称矩阵方程组Ax=b的具体算法如下：

1）给定误差极限ε>0,初始解向量x⁰∈Rⁿ;

2）计算r⁰=b-Ax⁰,β=||r⁰||，v₁=r⁰/β;

3）for k＝1to n do

β₁=0

for i=1to m do

α_i=(Av_i,v_i)

v_t=Av_i-α_iv_i-β_iv_i-1

if(i<m)β_i+1=||v_t||,v_i+1=v_t/β_i+1

End for

4）用追赶法求解如下的m维的三对角方程组:

解得y,计算 $z^m=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{v}_j{y}_j;$

5）如果||r⁰-Az^m||<ε,取x=x⁰+z^m,跳出循环；否则x⁰=x⁰+z^m，转1）

End for

本发明涉及的是流体运动的模拟，所以计算区域就是流体运动所能达到的地方。比如说，水在水箱中流动，整个水箱就是计算区域，模拟是针对计算区域中流体进行的。

本发明采用OpenMP/MPI的混合模型完成了半隐式粒子类算法的并行。其中OpenMP模型是多线程的加速方案，通常应用于共享内存的多核服务器；而MPI是进程级别的加速方案，通常应用于通过网络连接的集群服务器的多个节点上。MPS是一种典型的半隐式粒子类方法，并行化MPS具体的实施步骤如下：

1）设置计算区域内粒子的初始分布，并完成计算参数（如r_e）的初始化；如图2所示，将计算区域划分为多个子区域，将每一个子区域分配给一个进程负责计算，通常每一个进程都是位于一个多核集群服务器节点上；

2）将进程（子区域）上靠近划分边界的一列背景网格内的粒子信息（如图3所示）发送给邻居进程（子区域），并且邻居进程上要完成相应的接受操作；发送和接受的操作按如下步骤完成：

a）根据粒子的位置信息选出靠近划分边界处一列背景网格内的粒子（如图3所示）并通过MPI_Type_indexed()函数将这些粒子设置为发送粒子类型，然后将发送粒子类型的粒子信息采用非堵塞发送函数MPI_Isend()发送往相应的邻居进程；

b）在本进程中利用MPI_Type_indexed()函数在空闲内存上创建接受粒子类型，采用MPI中的堵塞接受函数MPI_Recv()接受相应邻居进程发送的信息，确保每一次发送操作在邻居进程上都对应一次接受操作；

3）在每一个进程上显式计算外力源项和黏性项，然后计算临时速度场和临时位置场，从而计算压力泊松方程的源项；这些操作是通过for循环针对子区域上的每一个粒子进行的，采用OpenMP模型中的编译制导语句通过派生多线程完成对for循环的并行，如下所示：

#pragma omp parallel for schedule(static)shared(…)private(…)

其中，schedule(static)指明的是采用静态方式来分配并行任务，这可以自动完成多线程间的负载均衡；shared(…)指明共享变数列表，private(…)指明的是每个线程的私有变数列表；

4）采用并行对称兰乔斯算法（SLA）全局求解压力泊松方程，获得压力值；这里的“并行SLA算法”指的是采用OpenMP/MPI的混合模型来加速的算法，具体来说，在每一个进程上的for循环操作都采用如步骤3所示的OpenMP方案来加速，而SLA算法中的用于计算向量内积的归约操作用OpenMP中的reduction()子句完成；在进程之间每一次计算粒子间相互作用之前，都采用步骤2所示的MPI接受发送方案来完成通信，内积操作采用MPI_Allreduce()函数完成归约；

5）采用与步骤2相同的MPI接受发送方案将获得的压力发送到邻居进程上；

6）传送压力的通信完成之后，采用类似于步骤3的OpenMP加速for循环的方案并行计算压力梯度项，用于修正临时速度场，从而获得最终的速度场；

7）每隔一定的迭代次数（20-100）调节区域划分边界，以保证每个进程上的负载都近似相同，具体的调节步骤如下：

a）统计本进程上要计算的粒子的数目；

b）如果右侧进程（子区域）上的粒子数目明显多于本进程计算的粒子数目，则将本进程的右边界向右移，反正左移；

c）将本进程的右边界采用类似步骤2的通信操作发送给本进程右侧的进程，作为其左边界；

8）判断本进程上由于粒子移动产生的粒子越界的情况，将已经越过本进程边界的粒子信息采用类似步骤2的操作发送给相应的进程；

9）判断是否已达到模拟设定时间，如果已达到模拟设定时间则退出；如果没有达到模拟设定时间，返回步骤2，继续下一个时间步长的迭代。

半隐式类粒子法中除去计算压力的其他操作多是针对单个粒子独立进行的，这部分操作本来就是易于并行的；求解压力的对称兰乔斯算法主要由矩阵向量乘法、向量的线性运算和向量的内积运算组成，这些操作也是易于并行的。因此在引入对称兰乔斯算法（SLA）求解压力泊松方程之后，半隐式类方法就很容易实现并行化。本发明首先通过粒子划分的方法采用OpenMP模型在多核的服务器节点上完成了并行；为了进一步提高计算的规模，本发明又通过区域划分的方法采用MPI模型在多节点的服务器上完成了并行。从而使得半隐式类粒子法计算千万粒子数目的工况成为可能。

为了展示并行的效果，本发明采用并行后的程序模拟了如图4所示的溃坝问题。首先，图5说明了SLA算法在模拟溃坝问题中的计算效率比ICCG方法还要高；同时图6则说明了SLA算法中易于并行的部分所占的比例远远高于ICCG算法；因此，SLA算法是一种高效且易于并行的算法，很适合半隐式粒子法的并行；其次，图7给出了采用OpenMP模型的程序在采用不同个数的CPU核时的并行效率；最后，图8给出了采用MPI并行后程序的在不同计算节点下并行效率。

参见图1，为了减少计算量，每一个粒子只与周围有限个粒子发生相互作用。图中r_e是一个粒子的有效作用半径，该粒子只与在r_e为半径的圆内的粒子发生相互作用。为了减小粒子搜索邻居的时间，首先在整个计算区域内建立以r_e为间距的背景网格，每一个粒子的邻居只需要在它周围的九个背景网格（图1中阴影所示）内搜索（三维的情况是在周围的27个背景网格内寻找）。

参见图2，在MPI并行模型中，首先将计算区域按照图1中的背景网格线划分成紧密相邻的多个子区域；然后将子区域分配给不同的进程，由每一个进程负责每一个子区域内粒子的计算。

参见图3，当一个进程负责一个子区域的计算之后，子区域边界处的粒子必然会与相邻的子区域内的粒子发生相互作用，因此子区域（进程）间的通信操作不可避免。如图1所示，一个粒子只与它紧密相邻的背景网格内的粒子相互作用，因此只需要将靠近边界处的一列背景网格中的粒子的信息发给与之相邻的子区域即可。

参见图4，溃坝问题的模拟用于展示并行的效果。在水箱的左端，有一用挡板固定的液柱，初始时刻，抽取挡板，液柱在重力作用下倒塌，沿水箱底面向右流动。液柱的高度H为0.5m，液柱的宽度L为0.25m。粒子的配置采用如下方式进行标示：dam N×M。其中dam表示算例属于溃坝问题的模拟，N表示在宽度（图4中L方向）方向采用N个粒子进行离散，M表示在高度（图4中H方向）方向采用M个粒子进行离散。

参见图5，ICCG算法是在串行程序中的一种高效求解压力泊松方程的求解器。本算例采用的粒子配置为dam100x200，可以看出SLA算法甚至比ICCG算法有更高的求解效率。

参见图6，本算例采用的粒子配置为dam100x200，SLA算法中易于并行的部分所占的百分比要远远高于ICCG算法中的百分比，因此SLA算法也是更容易并行的。

参见图7，OpenMP的并行效率被定义为采用OpenMP多核并行执行时相对于单核串行执行的加速比与所用核数的比值。可以看到随着模拟规模的增大采用OpenMP模型并行的算法的并行效率也在不断提高；当10个CPU核用于并行时，并行效率最高可达75%。

参见图8，MPI的并行效率被定义为采用MPI多个节点加速时相对于单节点执行时的加速比与所用节点数的比值。可以看到随着模拟规模的增大，采用MPI模型并行算法的并行效率也在不断提高；当10个节点用于并行时，并行效率最高可达55%。

Claims (6)

1.一种半隐式类粒子法的大规模并行处理方法，其特征在于：该处理方法在并行服务器上采用并行模型并行加速半隐式类粒子法，半隐式类粒子法中的压力泊松方程采用对称兰乔斯算法求解。

2.根据权利要求1所述一种半隐式类粒子法的大规模并行处理方法，其特征在于：所述半隐式类粒子法包括移动粒子半隐式算法和不可压缩的光滑粒子动力学方法。

3.根据权利要求1所述一种半隐式类粒子法的大规模并行处理方法，其特征在于：所述并行服务器为对称多处理器、GPU服务器、由多CPU核的服务器组成的多节点的服务器、由多GPU核的服务器组成的多节点的服务器或由多CPU核和多GPU核服务器组成的混合多节点服务器。

4.根据权利要求1所述一种半隐式类粒子法的大规模并行处理方法，其特征在于：所述并行模型为共享内存的并行模型、消息传递的并行模型或二者的混合模型。

5.根据权利要求1、2、3或4所述一种半隐式类粒子法的大规模并行处理方法，其特征在于：所述并行加速具体包括以下几种优选情况：

a）如果并行服务器是对称多处理器，采用OpenMP模型、MPI模型或采用OpenMP和MPI的混合模型来并行加速；

b）如果并行服务器是GPU服务器，通过OpenCL模型或CUDA模型并行加速；

c）如果并行服务器是多节点的服务器，每个节点是多CPU核的服务器，在每一个节点上，采用OpenMP模型并行加速，而在节点之间采用MPI模型来并行加速；

d）如果并行服务器是多节点的服务器，每个节点是多核的GPU服务器，在每一个节点上，采用OpenCL模型或CUDA模型并行加速，而在节点之间采用MPI模型来并行加速；如果节点中既有多核的CPU节点，又有多核GPU节点，对CPU节点采用OpenMP模型并行加速，对GPU节点采用OpenCL模型或CUDA模型并行加速，在节点之间仍采用MPI模型并行加速。

6.根据权利要求1所述一种半隐式类粒子法的大规模并行处理方法，其特征在于：所述处理方法的具体步骤如下：

1）设置计算区域内粒子的初始分布，将计算区域划分为多个子区域，将每一个子区域分配给一个进程负责计算；

2）将进程上靠近划分边界的粒子信息发送给邻居进程，并且邻居进程上要完成相应的接受操作；

3）在每一个进程上显式计算外力源项和黏性项，然后计算临时速度场和临时位置场，从而计算压力泊松方程的源项；步骤3）的操作是针对子区域上的每一个粒子在进程内并行加速；

4）采用并行对称兰乔斯算法全局求解压力泊松方程，获得压力值；

5）将获得的压力值发送到邻居进程上；

6）传送压力值的通信完成之后，并行计算压力梯度项，用于修正临时速度场，从而获得最终的速度场；

7）计算过程中要动态调节区域划分边界，以保证每个进程上的负载都近似相同；

8）判断本进程上由于粒子移动产生的粒子越界的情况，将已经越过本进程边界的粒子信息发送给相应的进程；

9）判断是否已达到模拟设定时间，如果已达到模拟设定时间则退出；如果没有达到模拟设定时间，返回步骤2），继续下一个时间步长的计算。

Images