SPH-GPU并行计算在风沙流中的应用

梁岚博，金阿芳，闻腾腾

新疆大学机械工程学院，乌鲁木齐 830047

风沙运动为主要标志的荒漠化及其带来的风沙环境是全球最突出的环境问题之一。风沙流活动对农业、牧业、沙漠公路交通、沙漠石油开采等造成了极大的危害，沙害造成的经济损失高达数百亿元人民币。而风沙两相流是一门既古老又年轻的科学，是认识风沙运动规律的重要基础，同时也是开展沙漠化治理和促进沙漠化区域可持续发展的科学依据，因此风沙两相流的研究具有重要的理论和实践意义[1]。

风沙两相流的研究方法主要有风洞模拟、野外观测和数值模拟。风洞试验观测是研究风沙运动最常用的手段，通常采用非接触测量仪器如粒子图像测速仪（particle image velocimeter，PIV）、相位多普勒粒子分析仪（phase Doppler particle analyzer，PDPA）、高速摄影仪等，可以对沙粒的微细运动和风洞输沙的时空变化进行分析。野外观测更多采用风速仪和集沙仪以及各种类型的输沙传感器等对风沙运动及输沙量进行测量[2]。随着计算机硬件和软件以及计算方法的发展，数值模拟近年来在风沙两相流研究中的作用越来越重要。数值模拟方法便于对风沙运动进行定量模拟和预测，易于控制参数，相对于风洞实验和野外观测有成本较低和便于实现的优点。

目前，流体数值模拟研究包括基于网格的欧拉法和基于粒子的拉格朗日法。在基于欧拉网格的数值方法中如FDM，首先需要在问题域上生成网格，要在不规则或者很复杂的几何形状上构造规律的网格是一件非常困难的事情，在流体模拟的整个过程中还会出现很多如大变形、变形边界、自由表面和物质交界面的问题，给基于网格的数值模拟方法带来了很大的挑战。拉格朗日无网格法主要有SPH、MPS和DEM等方法，以上粒子之间不需要进行网格划分（不是预先定义好的网格系统），是一系列任意分布的离散点，因此可以克服流体模拟中网格畸变和网格移动等引起的各种问题。目前，SPH方法被广泛应用于各种领域如溃坝、轮胎划水、飞机开沟、血流血管和船的撞击等[3-4]，且该方法更加稳定和真实，但SPH方法运用在风沙两相流中还较为缺乏。文献[5]首次将SPH 方法运用在风沙流中，建立了基于SPH 方法的风沙流控制方程，提出了模拟过程中一些关键参数的选择，并使用虚粒子法阻止粒子穿透边界。文献[6]利用文献[5]的算法，对沙粒的跃移运动进行了模拟，并在模拟过程中引入了起风和停风模块，分析了在风的加载和卸载时，沙粒的运动轨迹。文献[7]利用文献[5]的算法，建立了具有坡面形态的计算模型，得到了坡面不同位置沙粒的跃移高度和水平距离，验证了沙粒跃移的平均距离是单个沙波纹长度的1.16倍。文献[8]针对欧拉-欧拉流体计算模型在求解气沙两项耦合流动问题中存在的缺陷，提出了基于SPH-FVM 耦合的方法模拟风沙运动，用SPH 法离散沙粒相进行求解，用有限体积法（finite volume method，FVM）求解气相粒子，气沙两相通过拖曳力以及体积分数等参数进行耦合，模拟了自由来流风作用下沙粒的运动过程以及沙丘缓慢向前运动的过程。由于以上模拟过程中需要将计算区域离散成数量庞大的单个粒子，其计算时间长和模拟效率低一直是该方法的瓶颈。

近年来，随着科学技术的迅速发展，科学研究领域对计算机性能的要求也在不断增长，更是对计算能力提出了极高的要求[9]。传统计算机的单核处理器CPU，已经远远不能满足计算要求，即使是目前广泛使用的多核处理器CPU，在庞大的计算量压力下，也捉襟见肘[10]。为了满足巨大计算需求，众核架构图形处理器（GPU）的出现成为了必然。目前，由于GPU 超强的计算能力且适合计算密集型和易于并行的程序，以广泛应用到SPH 水体仿真算法中[11-12]，但应用到风沙运动中还比较匮乏。

因此本文引入了GPU 众核架构，并且通过对SPH串行算法的改进，提出了一种SPH-GPU并行计算方法，用以在计算机上快速的模拟风沙运动。

1 风沙流的SPH描述和离散

1.1 SPH基本思想

SPH 方法是在1977 年由Gingold 和Monaghan 等人[13-14]提出，最初用于解决天体物理学问题，后来经过其他学者的不断改进和完善，推广到流体力学各个方面。该无网格方法是将连续的流体离散成一系列任意分布的SPH 粒子，每个粒子携带质量、速度、密度等物理信息，通过求解SPH控制方程来描述整个系统的演变过程。

1.2 SPH基本方程

SPH方法的构造按两个关键步骤进行：第一步为积分表示法，又称场函数核近似法；第二步为粒子近似法[15]。场函数核近似法如下式：

式中，Ω是支持域，f(x)为三维坐标向量x的函数，δ(x-x′)是狄拉克函数，性质如下：

用光滑函数W(x-x′)取代f(x)积分表示式中的δ(x-x′)则：

式中，角括号＜＞在SPH的习惯用法中是核近似算子的标记；W(x-x′)即光滑函数；h是光滑长度，用来定义光滑函数影响区域的范围。

式（3）可转化为支持域内所有粒子叠加求和的离散化形式，称为粒子近似法。在粒子i处的粒子近似形式可写为：

其中，N是粒子i的支持域内的粒子总数；xi、xj分别为i和j粒子的坐标向量；ρj、mj为j粒子的密度和质量。

1.3 N-S方程的SPH离散

1.3.1 连续密度法

对连续性方程进行速度散度SPH近似：

1.3.2 动量方程的粒子近似法

对动量方程的梯度项直接应用SPH 粒子近似法进行变换可得：

1.4 邻近粒子的搜索

在模拟风沙运动的过程中，气沙两相粒子之间的位置是不断变化的，所以需要在每个时间步内都要重新计算每个粒子作用域中所包含的粒子，可以说邻近粒子搜索法的选择是非常重要的一步。

图1（a）全配对搜索法是最简单和直接的搜索方法。对于目标粒子i（实心红色圆点），要计算它到问题域中每一个粒子j(j=1,2,…,N,N是问题域中的粒子总数，图中用蓝色实心圆点表示）之间的距离，若距离小于粒子的支持域的半径kh，则粒子j为粒子i的邻近粒子（所有在红色圆圈之内的点）。在本文中因考虑对称光滑长度，则粒子i和粒子j为一组相互作用的粒子对，即粒子i也在粒子j的支持域内，那么，在搜索粒子j的邻近粒子时就不需要再考虑粒子i，因此可以减少一半的搜索计算量。该搜索方法对计算域内的所有粒子进行，虽然这个方法的应用较简单，但显然，这种方法的时间复杂度是O(N2)。

图1 邻近粒子的搜索Fig.1 Search for nearby particles

在运用链表搜索法时，在问题域上要铺设一临时网格，如图1（b）所示。将每个粒子都分布在网格单元内，并通过简单的存储规则将每个网格内的所有粒子连接起来。网格单元的大小与粒子支持域的大小相同，此网格只是用来划分区域，不参与计算。即，若粒子i的支持域大小为kh，则网格单元的大小也为kh。那么，对于给定的粒子i，其邻近粒子只能在同一网格或相邻网格里，由于本文的问题域是二维的，所以相邻网格的数量为8。若每个单元内的平均粒子数量很小，则链表搜索法的复杂度为O(N)。

上述算法在构建相邻粒子对时所使用到的输入数据仅与当前粒子（沙粒、气体以及虚粒子）的空间位置有关，与其他变量无关。每次计算输出的结果是粒子对的信息，并不改变当前时刻的粒子位置数据。在构建粒子对循环函数体中计算出的结果不影响输入数据，那么函数体中每次迭代的过程都不相互依赖。因此任一时刻某个粒子在寻找相邻粒子对时都是相互独立的，也就是说让GPU中的每一个线程携带粒子空间位置信息进行计算，每个线程执行时的区别仅是输入位置数据的不同，而线程的计算指令保持一致，符合GPU大规模并行计算的要求。

2 计算条件

文中主要是基于气沙两相耦合之间的受力机制来模拟沙粒的运动过程。模拟过程中涉及的各类参数及运行环境分别为表1和表2所示。由于在边界上的粒子积分时会被边界截断，会导致求解结果出错，所以在计算区域上部和下部采用固体边界，进出口采用周期性边界。在计算区域顶部和沙粒的底部分别设置一层虚粒子，从而阻止粒子非物理穿透边界如图2所示。初始床面上的风速廓线公式：

图2 计算区域初始状态示意图Fig.2 Schematic diagram of initial state of calculation area

表1 沙相、气相物性和计算参数Tabel 1 Sand，gas phase physical properties and calculation parameters

表2 计算环境Tabel 2 Computing environment

式中，μz为床面高程z处的风速；z0为沙床面粗糙度，z0=Ds/30。

3 SPH方法在CUDA计算架构下的实现

为了实现并行计算，需要对SPH 串行程序进行分析和修改，首先寻找程序中的热点函数（程序中耗时最多的函数），此外还要考虑数据的依赖性（前面的计算结果紧接着被后面的函数用到）分析哪部分程序适合串行和并行计算等等。最后用C语言编写主流程，采用CUDA C 语言调用GPU 设备实现并行计算，使改进后的程序，能使CPU与GPU协同合作完成任务，最大限度的发挥计算机的并行处理能力，从而提高SPH数值计算效率。

如图3 所示模拟过程分为4 个部分：读取需要模拟的数据、搜寻邻居粒子、计算粒子间作用力以及粒子位置的更新。CUDA 程序开始计算时，CPU 程序占主导地位，主机端创建输入和输出数组，为输入数据和结果提供存储空间，将已有的粒子信息文件（初始的位置、速度、密度等属性值）读取到CPU 内存里；随后CPU加载cudaMalloc函数分配GPU显存和cudaMemcpy函数把主机内存里的粒子信息拷贝至GPU 显存里，接下来CPU 加载核函数给GPU 做计算，GPU 计算时CPU阻塞不执行其他任务，核函数在设备端进行相邻粒子搜索，搜索之后计算粒子属性值，在计算粒子属性值时要考虑数据之间的依赖性，具有依赖性的数据不适合进行并行计算，否则会导致计算结果不准确甚至程序崩溃；最后进行粒子位置更新，更新完成后，跳至邻居粒子搜索进行循环执行，在完成所有计算任务之后，才将所得结果拷贝到内存中，这样就减少了内存和显存数据传输所花费的时间，提高了并行计算的效率。这个模拟过程中还要考虑合理利用设备资源，比如使SM 始终保持计算状态，合理的分配线程，选择合适的存储器等。

图3 GPU加速的SPH计算流程Fig.3 Accelerated SPH calculation process through GPU

4 CPU-GPU并行计算结果与分析

4.1 SPH串行算法热点分析

并行算法实现之前，要对整个串行程序进行热点分析，在进行并行加速时，要考虑串行程序中计算时间很长的程序段，也就是热点程序。如果首先对热点程序进行并行实现，那么对整个程序计算效率的提升是非常巨大的。上结所说模拟过程中耗时的部分主要有3个：搜寻邻居粒子、计算粒子间作用力以及粒子更新。

如图4是模拟25 000步得到的结果，从图中可以很明显地看出，整个程序耗时最大的部分在于搜寻邻居粒子上，这部分占据了整个程序一半以上的时间，是串行程序的主要热点所在。计算粒子间作用力，花费了接近40%的时间，也是热点程序的部分。相比较而言，粒子信息更新在总时间中所占据的比例仅仅2%。以上可得，在进行并行算法实现时，应该首先分析搜寻邻居粒子这段程序上，寻找更好的粒子搜寻方法，其次在优化计算粒子间相互作用力的程序，最后对于粒子更新可以选择保留串行计算，也可以进行并行计算。

图4 SPH串行热点程序Fig.4 SPH serial hotspot program

4.2 不同搜索方法的串行SPH计算时间

由上一节得出在搜寻邻近粒子上的耗时占用了整个程序的大部分时间，因此表3比较了全配对法和链表法所占的耗时，由于时间复杂度由O(N2)降低到O(N)，可见链表法是全配对搜索法计算速度的2.3倍。

表3 不同搜索法的SPH计算时间Tabel 3 SPH calculation time for different search methods

4.3 GPU并行计算模型验证

4.3.1 沙粒群运动的时空变化

本算例重点讨论流体起动，不考虑碰撞和冲击产生的瞬时冲击力的作用。

如图5 给出了在SPH-GPU 并行计算下时间步为1 500步图（a），5 100步图（b），8 400步图（c），9 900步图（d），19 200 步图（e），21 600 步图（f）的计算结果。从图中可知，沙粒的起动是在重力、气流拖曳力等各种力共同作用的结果，计算结果较好地描述了沙粒群的运动规律。如图（a）、（b）所示，在气流的作用下，平坦沙床上层沙粒产生加速度，当风速足够大时，沙粒从地面跳起开始运动。如图（c）、（d）中椭圆区域所示，随着时间的推移，当时间步进行到8 400步时，沙床此时开始发生松动或堆积，沙床中有些沙粒会发生离开沙床起跃的趋势，也有些沙粒其速度朝下，这些粒子会造成沙床更加紧实，沙床的松动和紧实是造成沙床产生沙波纹的直接原因之一。如图（e）、（f）所示，起跳沙粒在气流中不断的加速，又在重力作用下有了明显的下落。当沙床面上起跳的沙粒数与落回床面的沙粒数相差不大时，风沙运动进入了比较稳定的状态。

图5 沙粒群运动的时空变化Fig.5 Temporal and spatial changes of movement of sand particles

4.3.2 沙粒的运动轨迹

如图6为床面上典型沙粒的跃移运动轨迹，从图中可以看到，沙粒的运动轨迹具有上升段比较急剧，降落段比较平缓的特点，都呈现出一种非对称的抛物线形状，观察到靠近右周期边界的跃移轨迹的下落段并没有显示完整，这是由于进出口采用周期性边界，右边界出去的粒子从左边界进来，以上沙粒运动轨迹与文献[16]风洞实验中高速摄影记录的沙粒跃移运动相符。分析这是由于沙粒在跃移过程中始终受到气流拖曳力的作用，导致了在水平方向的分速度越来越大，当沙粒飞到一定高度，又在气流阻力和重力的作用下，迅速向下运动，因此导致降落段的水平位移较大。

图6 典型沙粒跃移轨迹Fig.6 Typical sand grain trajectories

如图7 表示沙粒运动中具有抛物线形状但同时又有些变异的轨迹，其跃移轨迹并不是非常光滑的抛物线形状，而是在转折处出现尖角，即沙粒会急剧降落，造成此现象的原因可能是局部气流纵向脉动较大，也可能是沙粒在空中与其他沙粒发生碰撞后获得了一个向下的冲量。

图7 变异的尖角轨迹Fig.7 Mutated sharp trajectories

4.4 GPU-CPU计算效率的对比

如图8 显示了粒子数分别为2 500、10 000、40 000来计算不同粒子数下GPU与CPU的效率，可以看出，使用不同的邻近粒子搜索法在GPU中的计算效率远高于在CPU 中的计算效率。GPU 并行链表搜索法相较于CPU 链表搜索法最高可以获得3.1 倍的加速比，而相对于CPU全配对搜索法最高可以获得20倍的加速比。当粒子数目为2 500 时，GPU 并行计算的加速效果并不明显，其原因是GPU的单线程计算能力远弱于CPU，随着粒子数目的增加GPU多线程并行计算能力的优势就发挥了出来，能够带来更大的加速比。

图8 CPU-GPU计算时间Fig.8 CPU-GPU calculation time

为了更好地理解CPU 与GPU 的加速机制，以10 000 粒子为例，统计了搜寻邻居粒子、计算粒子间作用力以及粒子更新在计算时占总时间的比例。如图9所示，GPU并行计算相较于CPU计算时，搜寻邻居粒子和粒子更新所占时间比例都减少了，然而计算粒子间作用的时间比例增加了。综合以上，可以推测出搜寻邻居粒子和粒子更新相对于计算粒子间作用力有更好的并行计算效率。

图9 SPH各部分占总时间的比例Fig.9 Proportion of each parts of SPH in total time

5 结论

本文首先对串行SPH程序进行分析和修改，最后采用SPH-GPU大规模并行加速技术对风沙流动过程进行数值模拟，结果表明：

（1）对SPH 串行算法代码进行了热点程序分析，得到在搜寻邻居粒子和计算粒子间作用力时耗时最长，因此优先优化这两部分程序段。由于搜寻邻居粒子是最耗时的部分，比较了全配对搜索法和链表搜索法的计算效率，得到链表法是全配对搜索法计算速度的2.3 倍。

（2）对SPH-GPU 并行计算模型进行验证。宏观上得到了沙粒群在各种力的作用下有起跳、上升和回落三个阶段的时空变化规律以及沙波纹形成过程，微观上得到了典型的沙粒跃移轨迹和变异的尖角轨迹。

（3）比较了CPU 与GPU 的计算效率，得到GPU 大规模并行计算的执行时间远远低于CPU 的执行时间。在粒子数较少时GPU 加速效果并不大，随着粒子数的增加，加速效果越加明显，最高可以获得20 倍的加速比。为了更好地理解CPU 与GPU 的加速机制，比较了搜寻邻居粒子、计算粒子间作用力以及粒子更新在计算时占总时间的比例，推测出搜寻邻居粒子和粒子更新相对于计算粒子间作用力有更好的并行计算效率。