基于GPU的面向SPH流體模擬的鄰居查找算法*

趙相坤 李鳳霞 戰(zhàn)守義

(北京理工大學北京市智能信息技術實驗室，北京100081)

基于物理的模擬已廣泛用于計算機圖形學中，越來越多的自然現(xiàn)象如雨、泥漿、煙圈、海浪等都可以被逼真地模擬顯示出來.這些模擬問題都可以歸結為復雜的流體模擬.一個典型的流體模擬可以用Navie-Strokes方程表示:

式中，u為流體速度，ρ為密度，p為壓強，ν為動力粘滯系數(shù)，F(xiàn)為作用于流體的所有外力之和.

目前主要有兩種跟蹤流體移動的方法:(1)歐拉方法，該方法一般基于網(wǎng)格，在固定的網(wǎng)格點監(jiān)測流體的參數(shù)(密度、速度、壓強等)隨時間的變化;(2)拉格朗日方法，該方法用一系列的粒子表示流體，流體在指定位置的參數(shù)通過對一定范圍內(nèi)粒子屬性值的加權求和得到.

平滑粒子動力學(SPH)［1］是一種拉格朗日形式的流體模擬方法.在SPH方法中，流體用一系列的攜帶有各種流體屬性的粒子表示，任意位置的流體屬性通過對附近的粒子屬性值進行加權求和得到.當使用的粒子數(shù)達到幾十萬甚至上百萬時，查找指定位置的粒子的所有鄰居就成為流體模擬的瓶頸.

為了查找粒子在指定位置處的鄰居，一般采用空間劃分的方法，將空間劃分成均勻的網(wǎng)格，如圖1所示.不同的粒子根據(jù)其位置屬性分別屬于不同的網(wǎng)格單元，位于同一網(wǎng)格內(nèi)的粒子互為鄰居，當然任一粒子的鄰居還包括與粒子所在網(wǎng)格相鄰的其余26個網(wǎng)格內(nèi)的粒子.

圖1 三維空間的網(wǎng)格剖分Fig.1 Grid partition in 3D space

為了實現(xiàn)粒子鄰居的快速查找，已經(jīng)有很多使用圖形處理器(GPU)加速鄰居查找的算法.Amada等［2］將GPU用于SPH的加速計算，但由于所使用的鄰居查找算法是在CPU上實現(xiàn)的，每一次模擬都需要把鄰居表從CPU傳輸?shù)紾PU，使得該方法無法充分使用GPU，CPU和GPU間的數(shù)據(jù)傳輸幾乎抵消了后續(xù)的GPU計算方面的優(yōu)勢.文獻［3］中將一種稱為模板路由的數(shù)據(jù)結構用于最近鄰居的查找，該方法需要預先建立一個復雜的紋理軌跡.由于任意空間網(wǎng)格用多個紋元表示，因此用于較大計算區(qū)域時需要一個大的紋理.Kipfer等［4］使用均勻的網(wǎng)格和排序方法在GPU上實現(xiàn)了粒子間的碰撞檢測.Purcell等［5］使用基于排序的網(wǎng)格方法和模板緩沖來處理位于同一個單元內(nèi)的多個光子.Harada等［6］提出了一種基于SPH的流體模擬系統(tǒng)，該系統(tǒng)使用桶紋理來表示3D的網(wǎng)格結構并實現(xiàn)了一種高效的鄰居查找算法，但每一個網(wǎng)格內(nèi)最多只能有4個粒子.Bayraktar等［7］提出了一種在GPU上實現(xiàn)的基于網(wǎng)格的鄰居查找算法，復雜的紋理產(chǎn)生過程使得該算法雖然可在GPU上實現(xiàn)，但是速度低于同類基于排序網(wǎng)格的方法.

為此，文中提出了一種新的在GPU上實現(xiàn)的基于網(wǎng)格的鄰居查找算法.首先使用粒子初始位置信息建立粒子的網(wǎng)格紋理圖，然后對粒子的網(wǎng)格紋理按照粒子的網(wǎng)格id在GPU上排序，最后利用排序的網(wǎng)格紋理和原始的網(wǎng)格紋理在GPU上建立粒子的鄰接紋理圖，實現(xiàn)粒子鄰居的快速查找.

1 算法原理及步驟

為了實現(xiàn)粒子鄰居的快速查找，需要建立一張如圖2所示的鄰接紋理圖:水平方向上保存了特定粒子的所有鄰居，其中R存儲當前粒子的索引號，G存儲當前粒子的鄰接粒子索引號，B未用，故用-1填充;在垂直方向上粒子按索引號依次排列.

圖2 鄰接紋理圖Fig.2 Neighbor texturemap

圖3給出了基于GPU的鄰居查找算法流程圖，在開始時需要建立一張存儲粒子位置信息的紋理圖，用作網(wǎng)格紋理的輸入，以后所有的操作全部在GPU上實現(xiàn).按照鄰居表的建立順序，下面給出詳細的實現(xiàn)方法.

圖3 基于GPU的鄰居查找算法流程圖Fig.3 Flowchart of GPU-based neighbor search algorithm

1.1 建立粒子網(wǎng)格紋理圖

初始時的粒子位置紋理如圖4(a)所示，其中(xpos，ypos，zpos)為粒子的坐標，partIdx為粒子的索引號，A為紋理的透明位.為了能迅速地找到粒子的鄰居，通常將包圍流體的容器空間劃分為三維的網(wǎng)格，任意粒子都屬于某一個網(wǎng)格.

圖4 粒子位置紋理和網(wǎng)格紋理Fig.4 Particle position texture and grid texture

為了得到粒子所屬網(wǎng)格以及和所屬網(wǎng)格直接相鄰的其它網(wǎng)格，首先由粒子的位置信息計算出粒子的所屬網(wǎng)格，即以粒子的坐標除以單元格的大小得到每個粒子的網(wǎng)格位置信息(ix，iy，iz).按照式(3)進行相應的轉換:

式中:b、h和d分別為x，y，z3個方向上網(wǎng)格單元的寬度、高度和深度.按照

計算出粒子的一維網(wǎng)格id，其中Nx和Ny分別為水平和垂直方向上的網(wǎng)格總數(shù).

將粒子所屬容器的幾何信息(位置和空間大小)以及粒子位置紋理圖作為輸入?yún)?shù)傳遞給像素著色器，就可以很容易地建立一張如圖4(b)所示的粒子網(wǎng)格紋理圖.位置紋理和網(wǎng)格紋理都是按照粒子的自然順序排列的，即按照partIdx有序.

1.2 備份網(wǎng)格紋理圖

初始的網(wǎng)格紋理圖按照粒子的partIdx有序排列.為了快速查找位于同一網(wǎng)格內(nèi)的所有粒子，可以使用基于GPU的折半查找算法［8］.使用基于GPU的查找算法需要網(wǎng)格紋理中的粒子按照gridIdx有序，而在建立粒子鄰接紋理圖時，需要借助原始的網(wǎng)格紋理圖得到粒子的partIdx，為此在對網(wǎng)格紋理排序前需要事先備份網(wǎng)格紋理圖.使用像素著色器可以將原來的網(wǎng)格紋理作為輸入寫入到備份紋理中.

1.3 排序網(wǎng)格紋理圖

為了快速查找具有相同網(wǎng)格id的粒子，需要對網(wǎng)格紋理按照gridIdx進行排序.目前存在很多在GPU上實現(xiàn)的排序算法，如奇偶歸并排序［8］、Bitonic歸并排序［8］、改進的Bitonic歸并排序［9］、GPUSort［10］、桶排序［6］、快速排序［11］等.文中采用奇偶歸并排序方法，單步的奇偶歸并排序GPU著色程序可以參考文獻［8］.經(jīng)過排序后，具有相同網(wǎng)格id的粒子相鄰排列.

1.4 建立粒子鄰接紋理圖

圖5給出了生成粒子鄰接紋理圖的過程，其中傳輸?shù)臄?shù)據(jù)標示在箭頭頭部，數(shù)字表示執(zhí)行的先后順序.

圖5 鄰接紋理圖的生成過程Fig.5 Creation procedure of neighbor texturemap

圖2中鄰接紋理的R位的粒子索引從紋理坐標y計算得到.如果每張紋理圖存儲4096個粒子，那么紋理的R位保存的粒子索引號為

式中，P為頁序號，介于0～Np/4096之間，Np為粒子總數(shù).

在得到粒子索引號partIdx后，在備份的原始網(wǎng)格紋理圖上得到該粒子的網(wǎng)格(文中稱為直屬網(wǎng)格)序號gridIdx，這樣根據(jù)gridIdx，在排序后的網(wǎng)格紋理圖中，使用基于GPU的折半查找算法［8］就可以得到與當前粒子網(wǎng)格id相同的其它粒子首次出現(xiàn)的位置.

在CPU上預先計算出相鄰的27個網(wǎng)格id間的間隔，然后傳輸?shù)紾PU，疊加到當前粒子的grid Idx就可以得到任意粒子所屬的其余26個網(wǎng)格.

雖然得到一個粒子的27個網(wǎng)格(直屬網(wǎng)格及其26個相鄰的鄰居網(wǎng)格)內(nèi)粒子首次出現(xiàn)的位置gridIdxStart，如何得到后續(xù)的相鄰粒子以及得到后寫入到紋理的什么位置是GPU鄰居查找方法的關鍵.這是因為和CPU不同的是，在GPU上，像素著色器不允許將指定的值寫到紋理的任意指定的位置上.為此，文中提出了一種解決像素著色器不能將任意值寫入到紋理的指定位置的方法.

首先，假定任意粒子擁有相同數(shù)目的鄰居并且任意網(wǎng)格都擁有相同數(shù)目的粒子，這樣任意粒子都會在鄰接紋理中出現(xiàn)27次，使用式(6)將紋理的橫坐標x除以每一個網(wǎng)格內(nèi)的最大粒子數(shù)Nmax，p，就得到該粒子所屬的片段step Idx(介于0～26之間)，利用式(7)可以得到粒子在片段內(nèi)的偏移量offset，將這個偏移量疊加到利用折半查找算法得到的首次出現(xiàn)的位置，就可以得到位于同一網(wǎng)格內(nèi)的后續(xù)粒子的網(wǎng)格索引號.

式中，int(·)為取整.

將網(wǎng)格索引號轉化為紋理坐標，然后利用紋理坐標從已排序的網(wǎng)格紋理中的R位得到當前粒子(索引號為partIdx)的可能鄰居粒子partIdxNeighbor.這些可能的鄰居粒子在最終寫入到紋理之前還要做如下的判別:

(1)如果partIdxNeighbor和partIdx相等，則說明找到的粒子和當前粒子是同一個粒子而不是鄰居，否則就是潛在鄰居;

(2)如果是潛在鄰居，則求取兩者之間的距離，如果該距離小于平滑半徑則是鄰居，否則不是鄰居.

按照上面的思想，最終建立了一張如圖2所示的粒子鄰接紋理圖.

2 文中算法與其它算法的比較

2.1 與基于CPU的鄰居查找算法的比較

在CPU上可以采用二叉空間分割(BSP)樹、多維樹、三維網(wǎng)格等來實現(xiàn)鄰居的查找.文中采用三維網(wǎng)格來實現(xiàn)基于CPU的鄰居查找，該算法和文中提出的基于GPU的鄰居查找算法的時間比較如圖6所示.從圖6可知:(1)基于GPU的鄰居查找算法的性能明顯優(yōu)于基于CPU的鄰居查找算法;(2)隨著粒子規(guī)模的不斷增加，基于CPU的鄰居查找算法的查找時間迅速上升，而且其增幅要高于基于GPU的鄰居查找算法，說明基于CPU的鄰居查找算法更適用于小規(guī)模的鄰居查找，在粒子越多的情況下，基于GPU的鄰居查找算法性能的優(yōu)越性越明顯.

圖6 2種算法建立鄰接紋理圖所用時間的比較Fig.6 Comparison of time cost for creating neighbor texture map by two algorithms

除了查找時間上的差異外，將文中提出的基于GPU的鄰居查找算法應用于基于SPH的流體模擬，可以保證流體模擬完全運行在GPU上，從而解決了在CPU上建立鄰接表，在CPU和GPU之間頻繁傳輸數(shù)據(jù)的低效問題［2］.

2.2 與其它基于GPU的算法的比較

與目前兩個主流的基于GPU的鄰居查找算法——Bayraktar算法［7］和Harada算法［6］相比，文中算法主要做了如下改進:

(1)Bayraktar算法在生成網(wǎng)格紋理圖時，需要統(tǒng)計網(wǎng)格紋理圖中具有相同grid Idx的粒子總數(shù)和具有相同grid Idx的粒子首次出現(xiàn)的位置，并保存到紋理中，文中算法丟棄了這一步;

(2)同一網(wǎng)格內(nèi)的粒子數(shù)通過改變鄰接紋理圖的寬度靈活設置，彌補了Harada算法只允許同一網(wǎng)格最多有4個粒子的局限，提高了計算精度.

3 應用實例

采用文中提出的基于GPU的鄰居查找算法來實現(xiàn)完全基于GPU的SPH流體模擬，其流體模擬框架如圖7所示.圖8給出了文中算法使用6萬個粒子進行模擬時粒子在容器內(nèi)的波動效果，表1給出了在NVIDIA 9600顯卡上文中算法和Bayraktar算法于不同粒子數(shù)下的幀速率對比.從表1可知，在NVIDIA 9600顯卡上，文中算法可以獲得11.6 f/s的幀速率，而Bayraktar算法在同樣的硬件條件下只能達到2.8 f/s的幀速率，說明文中算法因在GPU鄰接紋理圖生成過程方面的改進而大大提高了模擬的幀速率.從表1還可以知道，文中算法的模擬速度大約是Bayraktar算法的3～4倍.

圖7 基于GPU的SPH流體模擬框架Fig.7 Framework of GPU-based SPH fluid simulation

圖8 使用文中算法模擬6萬個粒子在容器內(nèi)的波動效果Fig.8 Waving simulation results in a poolwith 6×104 particles using the proposed algorithm

表1 不同粒子數(shù)下兩種算法的幀速率Table 1 Frame-rates of two algorithmswith differentnumbers of particles

在同樣的硬件條件下，使用文中算法與Harada算法［6］進行了對比實驗.如果每個網(wǎng)格內(nèi)的最大粒子數(shù)只能是4，那么在粒子特別多時，位于一個網(wǎng)格內(nèi)的粒子往往超過4個，就會出現(xiàn)因沒有將粒子納入鄰接紋理圖而導致粒子逃逸容器的情況;增加網(wǎng)格剖分的粒度，可以不受最多只有4個粒子的限制，但會降低模擬的速度.因此文中提出的控制粒子鄰居的方法更加靈活.

圖9是一組將一個球丟入容器內(nèi)時使用6萬個粒子模擬的情況，流體的自由表面采用文獻［12］中算法構建.由于文獻［12］中算法運行于CPU上，在GPU上產(chǎn)生的粒子位置信息不得不從紋理中讀回CPU中，因此模擬的速度大大降低，在NVIDIA 9600顯卡和Intel Core2 Duo CPU E4400硬件條件下，可以得到2.3 f/s的幀速率，這和Harada算法［6］在同樣條件下的幀速率相當.

圖9 使用6萬個粒子模擬的球丟入容器內(nèi)的效果Fig.9 Simulation results of throwing a ball into a tank with 6×104 particles

4 結語

文中提出了一種新的基于GPU的適合SPH流體模擬的鄰居查找算法，該算法保證了SPH方法可以完全運行于GPU上，大大改進了傳統(tǒng)方法在CPU實現(xiàn)鄰居查找時的CPU和GPU間數(shù)據(jù)傳輸?shù)牡托栴}.實驗表明，和其它基于GPU的鄰居查找算法相比，文中算法不但查找速度提高了3～4倍，而且鄰居粒子數(shù)最大值靈活可控，保證了算法的精度.

文中算法將鄰居粒子索引號寫入紋理位置，克服了GPU無法將指定數(shù)據(jù)寫入紋理任意位置的缺點，大大改進了GPU算法的性能.同樣，這種基于GPU的鄰居查找算法也適合于其它基于粒子的應用.但文中算法使用了大量的紋理緩存，在一定程度上降低了算法的效率.

［1］Müller M，Charypar D，Gross M.Particle-based fluid simulation for interactive applications［C］∥Proceedings of ACM SIGGRAPH/Eurographics Symposium on Computer Animation.San Diego:ACM，2003:154-159.

［2］Amada T，Imura M，Yasumuro Y，et al.Particle-based fluid simulation on GPU［C］∥Proceedings of ACM Workshop on General-Purpose Computing on Graphics Processors.Los Angeles:ACM，2004:1-2.

［3］Müller M，Solenthaler B，Keiser R，et al.Particle based fluid-fluid interaction［C］∥Proceedings of ACM SIGGRAPH/Eurographics Symposium on Computer Animation.Los Angeles:ACM，2005:237-244.

［4］Kipfer P，Segal M，Westermann R.Uberflow:a GPU based particle engine［C］∥Proceedings of ACM SIGGRAPH/Eurographics Conference on Graphics Hardware.New York:ACM，2004:115-122.

［5］Purcell T J，Donner C，Cammarano M，et al.Photon mapping on programmable graphics hardware［C］∥Proceedings of ACM SIGGRAPH/Eurographics Conference on Graphics Hardware.San Diego:ACM，2003:41-50.

［6］Harada T，Koshizuka S，Kawaguchi Y.Smoothed particle hydrodynamics on GPUs［C］∥Proceedings of Computer Graphics International.Berlin/Heidelberg:Springer-Verlag，2007:63-70.

［7］Bayraktar S，Güdükbay U，?zgüc B.GPU-based neighborsearch algorithm for particle simulations［J］.Journal of Graphics，GPU，and Game Tools，2007，14(1):31-42.

［8］Fernando Randima.GPU gems:programming techniques，tips and tricks for real-time graphics［M］.Boston:Addison-Wesley，2004:107-121.

［9］Pharr Matt，F(xiàn)ernando Randima.GPU gems2:programming techniques for high-performance graphics and general-purpose computation［M］.Boston:Addison-Wesley，2005:733-746.

［10］Govindaraju N K，Manocha D，Raghuvanshi N，et al.GPUSort:high performance sorting using graphics processors［R］.North Carolina:Department of Computer Science，University of North Carolina at Chapel Hill，2006.

［11］Cederman D，Tsigas P.A practical quick sort algorithm for graphics processors［R］.G?teborg:Department of Computer Science and Engineering，Chalmers University of Technology，2008.

［12］Lorensen W E，Cline H E.Marching cubes:a high resolution 3D surface construction algorithm［C］∥Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques.New York:ACM，1987:163-169.