基于CUDA的粒子系统的实现

基于CUDA的粒子系统的实现

用途:

这篇文章做为代码实现的先导手册,以全局的方式概览一下粒子系统的实现大纲.

科普:

对粒子进行模拟有两种基本方法:

1. Eulerian(grid-based) method : 相对于固定坐标系的点进行属性计算
2. Lagrangian(particle) method : 相对于移动坐标系的点进行属性计算

这是两个彻底不一样的坐标系参考下的方法,并且这两种方法各自发展出了本身的方法派系.

在刚体动力学中的应用：

1. 欧拉方法多被应用于固定坐标系下的动力学模拟,好比吊车,机械臂等.
2. 拉格朗日方法多被应用于移动坐标系下的动力学模拟,好比航天器等航天器.

由于在各自不一样的领域两种坐标系都能以最简方程式表达动力方程,因此两种坐标系下的方法并不打算合并统一.
若读者不太了解的话请思考,xyz坐标系和极坐标系的关系.

Particle-based method(拉格朗日方法)对比欧拉方法在粒子模拟领域
优势:

1. 只有必要时才进行计算
2. 更小的储存空间和带宽,由于只需计算粒子位置的模型属性而不用计算整个空间内任意点的属性.
3. 不须要被限制在一个有限的盒子内(由于欧拉方法计算空间中任意点的状态,因此必须限制范围)
4. 能简单实现质量守恒(每一个粒子的质量固定)

缺点:

1. 每一个粒子须要储存本身的计算结果(当粒子很是多时须要很大的内存开销)
   幸运的是,这对于显卡这种高速且大量并行运算的架构来讲这种模拟方法正合适.

这篇文章会介绍如何实现一个粒子模拟系统.

实现：

咱们把实现过程分为三步:

1. 积分位移
2. 创建网格结构数据
3. 处理碰撞

1.积分位移

计算积分位移是一个简单的过程, 咱们采集粒子的属性(位置,速度),而后根据每一个粒子的位置和速度计算下一个时间步长中粒子移动到的位置.
咱们使用欧拉方法(固定坐标系)根据粒子受到的力计算其速度 F=mvt .t为时间步长. 碰撞检测随后在这一环节进行计算.

2.粒子间交互

若是粒子间没有交互(好比实现烟花,烟雾,火焰等游戏特效),那整个系统就太简单了.
固然咱们讨论的粒子间是有交互的,基本交互就是碰撞(固然还有粘滞,引力吸引等),咱们使用空间分割方法来进行碰撞检测(固然这也是惟一的高效碰撞检测方法,之后咱们单独讨论不一样碰撞算法的区别和实现)
但首先咱们认可一个基本事实:
*. 当距离越远做用力越小,距离越近做用力越大 (你有可能认为这不符合宏观现实,但请你考虑往分子中插入一个质子须要多大力?）
这样咱们就能够只计算给定粒子周围半径内其余粒子给它的协力.空间分割方法能高效的获取粒子的最近临(你能够暂且看成kd-tree)

3.均分空间网格

咱们使用均分网格做为空间分割方法(固然也可使用其余的空间分割方法)
均分网格吧空间分割成一个个 cell. 为了简单起见咱们认为 cell的大小只能容下一个粒子.那么每一个粒子最多被空间中8个 cell所包围.咱们暂且认为粒子间没有穿透.那么每一个 cell最多最多只能接触到4个粒子.

当一个粒子刚好在一个 cell的中心的时候,咱们须要检测空间中 3x3x3-1=26个格子内是否包含其余粒子.因此咱们把粒子中心点位置的 cell的索引号做为粒子的索引号,这样根据每一个粒子的 cell编号就能够轻易的找到附近是否有其余粒子,由于 cell的编号是有排序的.

还有一种实现方式是把粒子半径中碰到的全部 cell编号都记录下来(更多细节被储存),这样碰撞处理所需计算量就减小了,但创建这个索引模型所需的开销就变大了,而根据实践结果而言,增长模型创建的开销比减小碰撞处理的时间来得更大,这是一种得不偿失的作法.

由于每一个时间步长(你能够暂且看成每帧,若是你对时间步长没有概念的话),网格数据都要进行更新,你可能想到在原有的数据结构中去更新,而非从新创建,但一般来讲更新模型并不比从新创建模型须要的时间更少.

4.创建空间网格模型

4.1 基于原子操做的方法

在支持原子操做的GPU上(nvdia显卡计算能力超过1.1),咱们能够用相对简单的算法实现这个创建过程.原子操做容许多个threads同时更新 global内存块中的同一个数据,而没必要担忧其发生读写冲突和竞争.

咱们使用两个数组,存放在 global memory中：

1. gridCounters - 这个数组储存着每一个 cell中含有多少个粒子.且每帧都重置为0,而后从新计算.
2. gridCells - 这个数据中储存着每一个 cell中包含粒子的索引号,且具备最大包含粒子个数限制.

updateGrid 这个 kernel函数负责计算并更新每个网格结构. GPU上每一个 thread跟踪计算一个粒子.
计算流程包括:

1. 粒子所在 cell
2. 使用 atomicAdd 函数增长 gridCounters中对应 cell的 counter值.
3. 把本身(粒子)的索引号加入 gridCells中对应 cell之中(使用scattered global write(散列全局内存写入)技术)

咱们限制每一个 cell设置最大包含粒子数,这样它就不会致使内存溢出问题,这是为了防范粒子穿透等状况发生时程序崩溃,也为了防止程序性能降低.(当多个粒子处在同一个 cell时,写入相关 cell的操做会进行序列化,这会致使程序性能的下降.),固然也能够分两步走:

1. 只计算 counter值
2. 扫描所有 grid 并写入粒子索引.

区别在于第二步每一个 thread计算对象是 grid,因此不会存在排序写入的问题.
下图显示了一个例子:

4.2 基于排序的方法

另外一种方法使用排序而不是原子操做来创建网格结构模型.

1. 使用 kernel calcHash 计算每一个粒子所在 cell的索引哈希值.而后根据该哈希值进行排序.这样就得到了以 cell索引为顺序的列表,随后寻找并合并相同 cell id的项,并散列写入 cell start列表中,由于已经排序,因此不会发生写入冲突或者排队写入的问题.
   这种方法解决了原子写入的排队等待问题的瓶颈,但增长了排序算法的开销.(实践结果显示使用Satish, N., Harris, M方法
   (fast radix sort 快速树排序在 CUDPP库中实现并提供)排序方法的算法开销比使用原子操做致使的排队写入开销来的小,因此性能更加)
   注意:该方法在不支持CUDA散列内存写入技术的硬件上只能使用二叉树来实现.

5. 粒子碰撞

一旦咱们创建起网格空间模型,咱们就可使用它来加速粒子之间的交互计算,咱们使用DEM 方法计算碰撞. 这种碰撞方法由两个粒子间的弹簧阻尼函数计算,弹簧力使两个粒子分开,阻尼力使粒子粘滞.
每一个粒子的计算首先在以前创建的空间模型内寻找附近 27个 cell中是否有其余粒子,若是存在其余粒子就进行碰撞的计算.若是没有其余粒子,就直接由粒子的速度计算下一时间的位置并刷新便可.
咱们能够看到若是咱们不创建空间网格模型,若是咱们有1万个粒子,则须要进行1W x 1W = 1亿次计算,而一旦创建空间模型只需进行 1W x 27 = 27万次计算便可,差了将近 400倍. N^2 和 27*N 的算法复杂度.

性能

在大多数的硬件上, 4.2 排序法有更好的性能.由于排完序可直接进行碰撞检测,在列表中的顺序便是空间上的位置结构.
而使用原子操做的方法对粒子的分布更加敏感,以前咱们已经讨论过排队等待写入的问题.
在实践中可使用 texture memory 对 global array进行映射,来获取粒子的速度和位置数据,由于 texture缓存对二维读取是优化过的, 一个 warp(32条并行 threads)内的 threads刚好要读取 texture缓存中的数据地址很是相近的话, texture缓存会为读取操做提供最佳性能和带宽,这是 texture 缓存硬件设计上的优点所决定的. (使用 texture缓存为读取数据增长了45%的性能)
这更加说明了排序方法的优点,由于它就是按照位置来排序的.

总结

CUDA的散列内存写入技术为建立动态数据结构提供了强力的支持.排序使得内存连续读取特性获得了充分的发挥,这一系列的组合使得GPU能够在一个实时的环境下模拟大型粒子系统. 咱们还可使用如下这些算法来进一步优化程序性能：

1. Qiming Hou, Kun Zhou, Baining Guo, BSGP: Bulk-Synchronous GPU Programming, ACM TOG (SIGGRAPH 2008)
2. Ericson, C., Real-Time Collision Detection, Morgan Kaufmann 2005
3. Satish, N., Harris, M., Garland, M., Designing Efficient Sorting Algorithms for Manycore GPUs, 2009.
4. Joshua A. Anderson, Chris D. Lorenz, and Alex Travesset General purpose molecular dynamics simulations fully implemented on graphics processing units, Journal of Computational Physics 227 (2008)

---

脚注:

scattered-global-write 散列全局内存写入技术是CUDA引入的技术,一般来讲只有获取或者写入连续物理内存地址时,写入和读取操做才能被合并成为一个，由于处理器一个指令只能获取一个内存地址(好比SIMD),若是内存地址连续那么就能够一次性拷贝连续数据至内存(内存对齐后,intel的CPU芯片也支持最多4个 float4类型的 SIMD).但当地址不连续时由于只有一个地址,因此也只能一个个指令顺序进行读取,显然这种方式成为了其余部分的性能瓶颈. 散列全局内存写入就是为了解决那些不连续内存写入的问题,它把多个不连续的内存写入合并为一个命令,而后并发写入,解决了这个瓶颈,固然也多亏了Nvidia显卡的MIMD架构(多指令多数据流).
DEM Harada, T.: Real-Time Rigid Body Simulation on GPUs. GPU Gems 3. Addison Wesley, 2007
Satish, N., Harris, M., Garland, M., Designing Efficient Sorting Algorithms for Manycore GPUs, 2009.
Ian Buck and Tim Purcell, A Toolkit for Computation on GPUs, GPU Gems, Addison-Wesley, 2004

参考文献:

1. Reeves, W. T. 1983. Particle Systems—a Technique for Modeling a Class of Fuzzy Objects. ACM Trans. Graph. 2, 2 (Apr. 1983), 91-108.
2. Monaghan J.: Smoothed particle hydrodynamics. Annu. Rev. Astron. Physics 30 (1992), 543. 12, 13
3. Müller M., Charypar D., Gross M.: Particle-based fluid simulation for interactive applications. Proceedings of 2003 ACM SIGGRAPH Symposium on Computer Animation (2003), 154–159.
4. Müller, M., Heidelberger, B., Hennix, M., and Ratcliff, J. 2007. Position based dynamics. J. Vis. Comun. Image Represent. 18, 2 (Apr. 2007), 109-118.
5. Harada, T.: Real-Time Rigid Body Simulation on GPUs. GPU Gems 3. Addison Wesley, 2007
6. Le Grand, S.: Broad-Phase Collision Detection with CUDA. GPU Gems 3, Addison Wesley, 2007
7. Nyland, L., Harris, M., Prins, J.: Fast N-Body Simulation with CUDA. GPU Gems 3. Addison Wesley, 2007
8. Z-order (curve), Wikipedia http://en.wikipedia.org/wiki/Z-order_(curve)
9. Ian Buck and Tim Purcell, A Toolkit for Computation on GPUs, GPU Gems, Addison-Wesley, 2004
10. Qiming Hou, Kun Zhou, Baining Guo, BSGP: Bulk-Synchronous GPU Programming, ACM TOG (SIGGRAPH 2008)
11. Ericson, C., Real-Time Collision Detection, Morgan Kaufmann 2005
12. Satish, N., Harris, M., Garland, M., Designing Efficient Sorting Algorithms for Manycore GPUs, 2009.
13. Joshua A. Anderson, Chris D. Lorenz, and Alex Travesset General purpose molecular dynamics simulations fully implemented on graphics processing units, Journal of Computational Physics 227 (2008)