基於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)