基于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核心間實現線程級并行?；赟PH方法的流體模擬包括粒子物理量計算及渲染等部分，所有部分都由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.