基于CPU-GPU混合加速的SPH流体仿真方法

胡鹏飞，袁志勇，廖祥云，郑 奇，陈二虎

(武汉大学计算机学院，湖北 武汉 430072)

1 引言

自然场景的真实感流体场景模拟在数字娱乐产业、军事、医学等领域中(如电影游戏特效、计算机动画和虚拟手术)有着广泛的应用，带人机交互功能的流体模拟还需是在线实时模拟。真实感流体模拟通常采用物理模拟方法，可分为基于网格的欧拉方法和基于粒子的拉格朗日方法两大类[1]。欧拉方法的坐标系是固定的，流体的速度、压力和密度等属性保存在这些固定点中，一般用于离线应用。而拉格朗日方法将流体视为流动的单元，用大量离散的粒子逼近模拟流体的运动，每个粒子都保有自身的相关物理属性，粒子间的相互作用通过光滑核函数来实现。近年来，基于粒子的拉格朗日方法得到了众多学者的青睐，取得了较好的研究进展。

Foster N和Metaxas D[2]介绍了一种运用纳维-斯托克斯N-S(Navier-Stokes)方程模拟流体的粒子方法，该方法是用来模拟基于物理的流体运动的一个常用的模型，但该方程计算量大。光滑粒子流体动力学SPH(Smoothed Particle Hydrodynamics)方法是基于粒子的流体模拟方法，运用SPH方法能大大简化N-S方程求解。尽管如此，随着流体模拟的真实性和实时性要求越来越高，用于模拟流体流动的粒子数量也越来越多，从而导致计算量大大增大，仅在CPU上进行SPH方法求解难以达到流体模拟的真实性和实时性要求。

随着图形硬件的发展，图形处理单元GPU(Graphic Processing Unit)不仅只负责图形渲染工作，还能采用GPU完成一些通用计算工作。在处理能力和存储器带宽上，GPU相对CPU有明显优势，且基于粒子的流体模拟问题具有高度的数据平行性和高强度计算需求，这类实时流体模拟问题非常适合采用GPU加速策略解决。

近几年来，CPU的发展也从单核发展到多核，多线程编程可在多个CPU核心间实现线程级并行。基于SPH方法的流体模拟包括粒子物理量计算及渲染等部分，所有部分都由GPU并行处理完成还存在较多的困难，组合CPU-GPU的混合加速流体模拟比传统的CPU或GPU计算更快且能简化算法设计。

基于上述考虑，本文提出一种基于CPU-GPU混合加速的SPH流体仿真方法。该混合加速方法包括流体模拟计算和渲染部分：(1)在流体模拟计算部分，采用GPU并行加速技术，并运用划分均匀空间网格方法加速构建邻居粒子表，同时对SPH各过程间的数据传输关联进行优化，以减少数据传输带来的加速损失；(2)在流体渲染部分，采用基于多核CPU的OpenMP加速技术对粒子采用点、球体、透明液体等进行渲染处理。

2 相关工作

Gingold R A等[3]提出的SPH方法是一种常用于模拟流体流动的计算方法，最初主要用于解决天体物理学中的流体质团无边界情况下的三维空间任意流动的计算问题。Stam J等[4]最早引入SPH方法模拟气体和火焰效果。Müller M等[5]利用SPH方法模拟流体的运动，首次将SPH方法应用于水面的绘制，提出了用一个光滑核函数来计算粒子间的密度、粘度、速度等物理量。之后的相关研究大多是基于Müller M的基础工作，并且相关改进主要集中于计算的简化以及稳定性的提高方面[6]。近几年的工作大多围绕不同模态的多相耦合展开，Schechter H等[7]提出在流体、固体周围设置虚粒子用于解决流-固耦合问题。Akinci N等[8]提出了一种可同时用于单双向耦合的流-固耦合边界处理方法。

SPH方法需要对用于模拟流体流动的离散粒子状态进行计算，需要大量的计算资源，如果仅依靠CPU处理计算工作，对于大数量级的粒子模拟，现有的CPU还不能完成逼真、实时的流体运动模拟。

随着图形硬件的发展，GPU通用计算GPGPU(General Purposed GPU)的概念被提出，GPGPU使得许多算法能在GPU上并行实现。Bolz J等[9]在GPU上实现了稀疏矩阵的求解。Moreland K等[10]实现了基于GPU的快速傅里叶变换。同时，研究人员也提出了一些采用GPU加速实现流体流动模拟的算法。Amada T等[11]用GPU实现了对SPH流体的模拟，但该算法中，邻居粒子表的构建是在CPU端进行的。Harada T等[12]提出了在GPU上构造均匀空间网格来覆盖整个离散粒子系统的计算区域，从而很大程度上加快了邻居粒子表的构建。

另外，多核处理器已成为未来处理器技术的发展方向，研究人员开始利用多核CPU来进行流体模拟。Ihmsen M等[13]提出一种基于多核CPU的并行计算方法对SPH粒子状态属性进行计算，大大提高了粒子的处理速度。多核多线程开发工具OpenMP[14]是一种共享存储并行计算机系统上的应用程序接口，它能通过与标准Fortran、C和C++结合进行工作，能有效支持同步共享变量以及负载的合理分配等任务且使用便捷。

3 SPH流体建模

SPH方法将流场离散成一系列粒子，把连续的物理量用多个离散粒子的集合的插值来表示，离散化形式为：

(1)

其中，A(r)代表密度、压力、粘度力等粒子的物理属性，r是Ω空间的任意一点，W是以h为作用半径的对称的光滑核函数，j是离散粒子系统中迭代所有的粒子质点，mj是粒子j的质量，Aj是粒子的某种属性，如粘稠力、密度等，ρj是粒子j的密度。

当使用密度值作为属性Aj时，得到各流体粒子的密度值为：

(2)

为使位于不同压强区的两个粒子相互压力相等，采用双方粒子压强的算术平均值代替单个粒子压力，得到：

(3)

对于单个粒子产生的压力p，可以用理想气体状态方程pi=K(ρi-ρ0)进行求解，ρ0表示流体的初始密度值，K是与流体温度有关的常数。

流体的粘性力通过公式(4)求出：

(4)

其中,μ是流体的粘度系数，vi是粒子i的速度。

单个粒子在初始状态下外力仅有重力，碰撞后，还会受到来自其他粒子以及容器的碰撞力和摩擦力：

(5)

通过上述力的分别计算，得到单个流体粒子的合力，再运用牛顿第二定律求出粒子运动加速度，通过加速度变化可以求出粒子速度变化，从而求出粒子位置改变。

4 CPU-GPU混合加速

基于CPU-GPU的混合计算系统是在现有CPU并行系统基础上增加了GPU作为辅助并行加速部件，从而充分利用系统内的可用资源，达到一种综合加速的效果。本文所述的基于CPU-GPU混合加速的SPH流体仿真框架如图1所示。该混合模拟仿真方法主要包括GPU流体计算和CPU流体渲染两部分。

Figure 1 Framework for SPH fluid simulation based on CPU-GPU hybrid acceleration图1 基于CPU-GPU混合加速的SPH流体模拟框架

4.1 GPU流体计算

SPH流体模拟方法采用离散的粒子来模拟连续的流体运动，为了得到具有丰富细节的流体，需有大数量级粒子参与模拟计算。而采用纯CPU计算方法进行流体迭代，不能达到实时的模拟效果。由于SPH模拟方法本身具有高度的并行性，它非常适合用GPU进行加速计算。借助GPU的高度并行性和可编程性，本文便能设计出在GPU上执行的SPH流体模拟并行算法，从而实现大数量级粒子的流体系统的逼真、实时仿真。

4.1.1 构建邻居粒子表

SPH流体的模拟是通过离散的粒子实现的，离散粒子间的相互作用通过一个对称的、具有一定作用半径的光滑核函数计算得到。粒子间的相互作用具有一定作用半径，对粒子的属性进行更新时，需遍历系统中的其他所有粒子，这会导致非常低的查找效率。为此，Amada T等[11]提出构造邻居粒子表的方法，采用该方法计算粒子物理量时，可只查找邻居粒子来进行计算。其优点在于每个迭代步，只需构建一次邻居粒子表，在粒子物理量计算过程中，都是查找表中数据进行计算。但是，采用这种方法构建邻居粒子表时，采用的方法仍是遍历所有粒子，导致建表时间过长，影响模拟的实时性。后来一些学者在原有邻居粒子算法基础上，进行了CPU端的优化工作，如基于粒子对的邻居粒子表构建，这种改进能使表的构建降低一半工作量。Harada T等[12]提出了构造空间均匀网格覆盖整个粒子位置空间的算法，设置均匀网格的半径与光滑核函数的作用半径相同，这样，在进行粒子邻居搜索时，只需要对当前粒子所在网格及周围26个网格进行搜索即可。这种改进极大提高了粒子的邻居表构建速度。

传统的粒子对邻居粒子表构建方法是基于CPU的，其建表过程如算法1所示：

算法1基于CPU的邻居粒子表构建

遍历每个粒子:

(1)将当前粒子i增加到邻居粒子表的i行第一个位置，并确定当前粒子i所在的网格g;

(2)遍历网格g及周围网格中的粒子j;

① 如果j不在i的影响域内，则遍历下一粒子;

② 如果j在i的影响域内，则i、j互为邻居关系，将j添加到i的邻居粒子表中。

(3)将当前粒子i添加到网格g中，确定粒子i与网格g的归属关系。

结束遍历。

该方法将确定粒子所在网格和搜索粒子邻居两项操作同步交错执行。如果直接移植到GPU运行，这种交错执行方式会涉及到线程同步及共享资源加锁问题，这些操作会极大地降低GPU加速效果，且该方法构造过程采用了邻居对方法，邻居粒子表最终存放的是只是单向的邻居关系。这种单向的邻居关系在后续粒子物理属性求解阶段会涉及资源加锁问题。

为此，本文设计出一种新颖的邻居粒子表构建算法，它将确定粒子所在网格与建立邻居粒子表两项操作分开串行执行，其中每项操作是在GPU中并行执行。这种处理策略易于两项操作的并行化并且所建立的邻居粒子表中存放的双向邻居信息，能对后续粒子物理属性进行并行求解，从而能够避免线程争夺资源。基于GPU的建立邻居粒子表算法如算法2所示。

Figure 2 Schematic diagram for SPH fluid simulation to optimize data transmission based on GPU acceleration图2 基于GPU加速的SPH流体模拟数据传输优化示意图

算法2基于GPU的邻居粒子表构建算法

操作1 确定粒子的网格归属:

1：for 粒子表中粒子pido

2： 确定粒子pi所在的网格g;

3： 向网格g中添加粒子pi;

4： 网格g容量增1;

5：end for

操作2 建立邻居粒子表:

6：for粒子表中粒子pido

7：neighbourList[pi][0].index=pi;

8：neighbourList[pi][0].dist= 0;

9： 确定粒子pi所在的网格g;

10： for 网格g周围网格g_neighbourdo

11： for 网格g_neighbour中粒子pjdo

12：dist=position[pi]-position[pj];

13： ifdist

14：nindex=neighbourList_Size[pi];

15：neighbourList[pi][nindex] =pj;

16：neighbourList[pi][nindex] =dist;

17：neighbourList_Size[pi] +=1;

18： end for

19： end for

20：end for

4.1.2 计算粒子属性

SPH流体模拟方法中粒子属性计算过程具有高度并行性，适合采用GPU并行计算。但是，该方法中属性间可能有依赖关系，如粒子的加速度属性是需要根据计算得到的密度和压强属性求得，这种依赖关系表明粒子的某些依赖属性间是不能并行计算的。对于没有依赖关系的粒子属性，其计算可以放到一个核函数中进行，以减少不必要的GPU调用开销；对于有依赖关系的粒子属性，其计算需放到不同的有先后关系的核函数中并行计算。

4.1.3 数据传输优化

拥有数十至数百个核心的众核GPU与几个核心的CPU相比，GPU能更快地进行SPH流体模拟计算且具有高效的显存读写能力。但是，一般基于GPU的并行计算程序中，CPU与GPU间的数据传输是并行能力提高的瓶颈。因此，本文在设计基于GPU加速的SPH流体模拟方法时，对CPU与GPU间的数据传输进行了优化处理。

基于GPU加速的SPH流体模拟数据传输优化如图2所示，图2中小矩形框表示SPH方法的粒子物理属性，椭圆形框表示SPH方法的主要计算工作。外层黑线框内表示数据流在GPU内部的流动过程，它不涉及CPU与GPU间的数据交互；在粒子的各种物理量中，最终需要读出的只有处于外层黑色框外的粒子的位置信息，但其他的速度、加速度、密度、压强等物理量信息不必读出而一直保存于GPU中，直接进行下一步计算，从而大大减少了CPU与GPU的数据传输次数。

4.2 CPU加速渲染

为了逼真地展现模拟流体的运动，渲染与可视化显得尤为重要，若采用球体对粒子系统渲染，在设定的光照条件下，便能清晰地展现出粒子与粒子间、粒子与容器边界间的碰撞过程。在现有的基于GPU加速的SPH计算方法中，通常是对流体的计算过程进行GPU并行加速，但渲染过程采用CPU串行执行。本文所述SPH流体仿真方法中，利用多核CPU的硬件条件并结合OpenMP实现了CPU加速渲染。

OpenMP采用了Fork-Join并行执行模式，如图3所示。Fork-Join并行执行模式的基本思想是，运行时库维护一个线程池，当主线程遇到并行结构时，从线程池创建一个并行的从线程组；当从线程组完成并行区域的工作后，各从线程间进行同步，从属线程返回到池中，主线程继续执行并行构造后的操作。

Figure 3 Fork-Join parallel execution model of OpenMP图3 OpenMP Fork-Join并行执行模型

5 实验结果

实验对64 k和128 k粒子数的球体渲染操作在CPU串行和OpenMP多核CPU并行环境下进行单帧平均耗时测试，结果对比情况如图4所示(实验平台的CPU具体配置是：双核CPU采用Intel® CoreTM2 Duo CPU E7500 2.93 GHz，四核CPU采用Intel® Xeon® CPU E3-1230 V2 3.30 GHz)。从图4可以看出，基于CPU的OpenMP加速渲染能使粒子球体的渲染效率成倍提高。

Figure 4 Comparison of rendering time using OpenMP图4 OpenMP加速前后单帧球体渲染时间对比

在双核CPU和带GPU显卡的微机实验平台(CPU：Intel® CoreTM2 Duo E7500 双核2.93 GHz；GPU显卡：ATI Radeon HD 5850,1 GB)上，通过对不同粒子数(128 k、64 k、27 k和8 k)SPH流体模拟方法中各部分操作进行数据传输优化得出如图5所示的数据对比关系。图5中四对两条相邻柱状图分别表示不同粒子数在单帧数据优化前后的CPU与GPU间平均数据传输时间消耗。从图5可以看出，通过对SPH流体模拟方法中各部分操作进行数据传输优化能很大程度上减小CPU-GPU间的数据交换消耗，从而加速SPH流体仿真的计算。

Figure 5 Comparison of transmission time between CPU and GPU after data optimization图5 数据优化前后CPU与GPU间单帧数据传输时间对比

粒子邻居搜索消耗是SPH流体仿真方法的主要消耗之一，通过对粒子网格归属及邻居粒子表构建过程进行改进并经GPU并行处理(实验平台同图5)，GPU并行前后粒子邻居搜索消耗如表1所示。表1是一个迭代步构建邻居粒子表在纯CPU端、GPU并行(未经数据传输优化)及GPU并行(经数据传输优化)的不同时间消耗对比，实验分别对8 k、27 k、128 k不同粒子数的单帧平均邻居搜索时间进行测试。

Table 1 Comparison of neighbor search time after parallel表1 并行前后邻居搜索时间对比

从表1可知，经数据传输优化的GPU并行邻居搜索效率比CPU串行邻居搜索效率提高了近千倍。

实验对流体从圆柱形容器高处下落到容器中的情况进行了SPH流体粒子模拟，模拟效果如图6所示(实验平台同图5)。图6a～图6C分别是粒子数为8 k、27 k、128 k时的CPU-GPU混合加速流体模拟效果。

Figure 6 SPH fluid simulation based on CPU-GPU hybrid acceleration with different particles number图6 CPU-GPU混合加速后不同粒子数流体模拟效果

表2给出了流体仿真中不同粒子数分别在CPU和CPU-GPU上进行10 k个系统迭代的平均帧率对比情况(实验平台同图5)。从表2可以看出，即使粒子数达到64 k的大规模粒子场景，本文所述方法仍能够达到30 FPS以上的实时模拟速度，与视频采集的帧速率相当。

Table 2 Comparison of frame rate with different particles number after parallel表2 不同粒子数并行前后帧率对比

为了更直观、逼真地展现加速后SPH流体仿真的效果，本文对SPH粒子采用不同方式进行渲染，图7所示是粒子数为27 k时的加速后某一时刻渲染效果图。图7a～图7c分别是采用点、小球和透明液体三种不同的渲染方式对SPH流体粒子进行渲染。

Figure 7 Renderings of different render way with 27 000 particles图7 粒子数为27 k时的三种不同渲染方式渲染效果图

6 结束语

本文在介绍SPH流体模拟国内外研究现状的基础上，提出一种基于CPU-GPU混合加速的SPH实时流体仿真方法，流体计算部分采用GPU并行加速,流体渲染部分采用基于CPU的OpenMP加速。通过前面的实验结果可知，运用GPU的并行计算能力相比于纯CPU实现而言，能实现更多粒子数的实时模拟；同时，利用基于多核CPU的OpenMP加速功能，能完成快速的渲染过程，从而实现逼真、实时的流体模拟仿真。

在未来的研究工作中，我们将研究CPU与GPU间的协同操作和负载均衡，以扩大粒子的并行规模；研究高质量的逼真实时流体渲染算法，并将本文所述方法应用于各种实际的流体模拟工程领域中。

[1] Liu M B,Liu G R.Smoothed particle hydrodynamics (SPH):An overview and recent developments[J]. Archives of Computational Methods in Engineering, 2010, 17(1):25-76.

[2] Foster N, Metaxas D. Controlling fluid animation[C]∥Proc of the 1997 International Conference on Computer Graphics, 1997:178-188.

[3] Gingold R A,Monaghan J J.Smoothed particle hydrodynamics-theory and application to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society, 1977,181:375-389.

[4] Stam J. Stable fluids[C]∥Proc of the 26th Annual Conference on Computer Graphics and Interactive Techniques, 1999:121-128.

[5] Müller M, Charypar D, Gross M. Particle-based fluid simulation for interactive applications[C]∥Proc of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2003:154-159.

[6] Gerszewski D,Bhattacharya H,Bargteil A W.A point-based method for animating elastoplastic solids[C]∥Proc of the 2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2009:133-138.

[7] Schechter H, Bridson R. Ghost SPH for animating water[J]. ACM Transactions on Graphics (TOG), 2012, 31(4):61:1-61:8.

[8] Akinci N, Ihmsen M, Akinci G, et al. Versatile rigid-fluid coupling for incompressible SPH[J]. ACM Transactions on Graphics (TOG), 2012, 31(4):62:1-62:8.

[9] Bolz J, Farmer I, Grinspun E, et al. Sparse matrix solvers on the GPU:Conjugate gradients and multigrid[J]. ACM Transactions on Graphics (TOG), 2003, 22(3):917-924.

[10] Moreland K, Angel E. The FFT on a GPU[C]∥Proc of the ACM SIGGRAPH/Eurographics Conference on Graphics Hardware, 2003:112-119.

[11] Amada T, Imura M, Yasumuro Y, et al. Particle-based fluid simulation on GPU[C]∥Proc of ACM Workshop on General-Purpose Computing on Graphics Processors and SIGGRAPH, 2004:1.

[12] Harada T, Koshizuka S, Kawaguchi Y. Smoothed particle hydrodynamics on GPUs[C]∥Proc of Computer Graphics International, 2007:63-70.

[13] Ihmsen M, Akinci N, Becker M, et al. A parallel SPH implementation on multi-core CPUs[J].Computer Graphics Forum, 2011, 30(1):99-112.

[14] Dagum L, Menon R. OpenMP:An industry standard API for shared-memory programming[J]. Computing in Science and Engineering, 1998, 5(1):46-55.