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

胡鵬飛，袁志勇，廖祥云，鄭 奇，陳二虎

(武漢大學(xué)計(jì)算機(jī)學(xué)院，湖北 武漢 430072)

1 引言

自然場(chǎng)景的真實(shí)感流體場(chǎng)景模擬在數(shù)字娛樂(lè)產(chǎn)業(yè)、軍事、醫(yī)學(xué)等領(lǐng)域中(如電影游戲特效、計(jì)算機(jī)動(dòng)畫(huà)和虛擬手術(shù))有著廣泛的應(yīng)用，帶人機(jī)交互功能的流體模擬還需是在線實(shí)時(shí)模擬。真實(shí)感流體模擬通常采用物理模擬方法，可分為基于網(wǎng)格的歐拉方法和基于粒子的拉格朗日方法兩大類[1]。歐拉方法的坐標(biāo)系是固定的，流體的速度、壓力和密度等屬性保存在這些固定點(diǎn)中，一般用于離線應(yīng)用。而拉格朗日方法將流體視為流動(dòng)的單元，用大量離散的粒子逼近模擬流體的運(yùn)動(dòng)，每個(gè)粒子都保有自身的相關(guān)物理屬性，粒子間的相互作用通過(guò)光滑核函數(shù)來(lái)實(shí)現(xiàn)。近年來(lái)，基于粒子的拉格朗日方法得到了眾多學(xué)者的青睞，取得了較好的研究進(jìn)展。

Foster N和Metaxas D[2]介紹了一種運(yùn)用納維-斯托克斯N-S(Navier-Stokes)方程模擬流體的粒子方法，該方法是用來(lái)模擬基于物理的流體運(yùn)動(dòng)的一個(gè)常用的模型，但該方程計(jì)算量大。光滑粒子流體動(dòng)力學(xué)SPH(Smoothed Particle Hydrodynamics)方法是基于粒子的流體模擬方法，運(yùn)用SPH方法能大大簡(jiǎn)化N-S方程求解。盡管如此，隨著流體模擬的真實(shí)性和實(shí)時(shí)性要求越來(lái)越高，用于模擬流體流動(dòng)的粒子數(shù)量也越來(lái)越多，從而導(dǎo)致計(jì)算量大大增大，僅在CPU上進(jìn)行SPH方法求解難以達(dá)到流體模擬的真實(shí)性和實(shí)時(shí)性要求。

隨著圖形硬件的發(fā)展，圖形處理單元GPU(Graphic Processing Unit)不僅只負(fù)責(zé)圖形渲染工作，還能采用GPU完成一些通用計(jì)算工作。在處理能力和存儲(chǔ)器帶寬上，GPU相對(duì)CPU有明顯優(yōu)勢(shì)，且基于粒子的流體模擬問(wèn)題具有高度的數(shù)據(jù)平行性和高強(qiáng)度計(jì)算需求，這類實(shí)時(shí)流體模擬問(wèn)題非常適合采用GPU加速策略解決。

近幾年來(lái)，CPU的發(fā)展也從單核發(fā)展到多核，多線程編程可在多個(gè)CPU核心間實(shí)現(xiàn)線程級(jí)并行。基于SPH方法的流體模擬包括粒子物理量計(jì)算及渲染等部分，所有部分都由GPU并行處理完成還存在較多的困難，組合CPU-GPU的混合加速流體模擬比傳統(tǒng)的CPU或GPU計(jì)算更快且能簡(jiǎn)化算法設(shè)計(jì)。

基于上述考慮，本文提出一種基于CPU-GPU混合加速的SPH流體仿真方法。該混合加速方法包括流體模擬計(jì)算和渲染部分：(1)在流體模擬計(jì)算部分，采用GPU并行加速技術(shù)，并運(yùn)用劃分均勻空間網(wǎng)格方法加速構(gòu)建鄰居粒子表，同時(shí)對(duì)SPH各過(guò)程間的數(shù)據(jù)傳輸關(guān)聯(lián)進(jìn)行優(yōu)化，以減少數(shù)據(jù)傳輸帶來(lái)的加速損失；(2)在流體渲染部分，采用基于多核CPU的OpenMP加速技術(shù)對(duì)粒子采用點(diǎn)、球體、透明液體等進(jìn)行渲染處理。

2 相關(guān)工作

Gingold R A等[3]提出的SPH方法是一種常用于模擬流體流動(dòng)的計(jì)算方法，最初主要用于解決天體物理學(xué)中的流體質(zhì)團(tuán)無(wú)邊界情況下的三維空間任意流動(dòng)的計(jì)算問(wèn)題。Stam J等[4]最早引入SPH方法模擬氣體和火焰效果。Müller M等[5]利用SPH方法模擬流體的運(yùn)動(dòng)，首次將SPH方法應(yīng)用于水面的繪制，提出了用一個(gè)光滑核函數(shù)來(lái)計(jì)算粒子間的密度、粘度、速度等物理量。之后的相關(guān)研究大多是基于Müller M的基礎(chǔ)工作，并且相關(guān)改進(jìn)主要集中于計(jì)算的簡(jiǎn)化以及穩(wěn)定性的提高方面[6]。近幾年的工作大多圍繞不同模態(tài)的多相耦合展開(kāi)，Schechter H等[7]提出在流體、固體周圍設(shè)置虛粒子用于解決流-固耦合問(wèn)題。Akinci N等[8]提出了一種可同時(shí)用于單雙向耦合的流-固耦合邊界處理方法。

SPH方法需要對(duì)用于模擬流體流動(dòng)的離散粒子狀態(tài)進(jìn)行計(jì)算，需要大量的計(jì)算資源，如果僅依靠CPU處理計(jì)算工作，對(duì)于大數(shù)量級(jí)的粒子模擬，現(xiàn)有的CPU還不能完成逼真、實(shí)時(shí)的流體運(yùn)動(dòng)模擬。

隨著圖形硬件的發(fā)展，GPU通用計(jì)算GPGPU(General Purposed GPU)的概念被提出，GPGPU使得許多算法能在GPU上并行實(shí)現(xiàn)。Bolz J等[9]在GPU上實(shí)現(xiàn)了稀疏矩陣的求解。Moreland K等[10]實(shí)現(xiàn)了基于GPU的快速傅里葉變換。同時(shí)，研究人員也提出了一些采用GPU加速實(shí)現(xiàn)流體流動(dòng)模擬的算法。Amada T等[11]用GPU實(shí)現(xiàn)了對(duì)SPH流體的模擬，但該算法中，鄰居粒子表的構(gòu)建是在CPU端進(jìn)行的。Harada T等[12]提出了在GPU上構(gòu)造均勻空間網(wǎng)格來(lái)覆蓋整個(gè)離散粒子系統(tǒng)的計(jì)算區(qū)域，從而很大程度上加快了鄰居粒子表的構(gòu)建。

另外，多核處理器已成為未來(lái)處理器技術(shù)的發(fā)展方向，研究人員開(kāi)始利用多核CPU來(lái)進(jìn)行流體模擬。Ihmsen M等[13]提出一種基于多核CPU的并行計(jì)算方法對(duì)SPH粒子狀態(tài)屬性進(jìn)行計(jì)算，大大提高了粒子的處理速度。多核多線程開(kāi)發(fā)工具OpenMP[14]是一種共享存儲(chǔ)并行計(jì)算機(jī)系統(tǒng)上的應(yīng)用程序接口，它能通過(guò)與標(biāo)準(zhǔn)Fortran、C和C++結(jié)合進(jìn)行工作，能有效支持同步共享變量以及負(fù)載的合理分配等任務(wù)且使用便捷。

3 SPH流體建模

SPH方法將流場(chǎng)離散成一系列粒子，把連續(xù)的物理量用多個(gè)離散粒子的集合的插值來(lái)表示，離散化形式為：

(1)

其中，A(r)代表密度、壓力、粘度力等粒子的物理屬性，r是Ω空間的任意一點(diǎn)，W是以h為作用半徑的對(duì)稱的光滑核函數(shù)，j是離散粒子系統(tǒng)中迭代所有的粒子質(zhì)點(diǎn)，mj是粒子j的質(zhì)量，Aj是粒子的某種屬性，如粘稠力、密度等，ρj是粒子j的密度。

當(dāng)使用密度值作為屬性Aj時(shí)，得到各流體粒子的密度值為：

(2)

為使位于不同壓強(qiáng)區(qū)的兩個(gè)粒子相互壓力相等，采用雙方粒子壓強(qiáng)的算術(shù)平均值代替單個(gè)粒子壓力，得到：

(3)

對(duì)于單個(gè)粒子產(chǎn)生的壓力p，可以用理想氣體狀態(tài)方程pi=K(ρi-ρ0)進(jìn)行求解，ρ0表示流體的初始密度值，K是與流體溫度有關(guān)的常數(shù)。

流體的粘性力通過(guò)公式(4)求出：

(4)

其中,μ是流體的粘度系數(shù)，vi是粒子i的速度。

單個(gè)粒子在初始狀態(tài)下外力僅有重力，碰撞后，還會(huì)受到來(lái)自其他粒子以及容器的碰撞力和摩擦力：

(5)

通過(guò)上述力的分別計(jì)算，得到單個(gè)流體粒子的合力，再運(yùn)用牛頓第二定律求出粒子運(yùn)動(dòng)加速度，通過(guò)加速度變化可以求出粒子速度變化，從而求出粒子位置改變。

4 CPU-GPU混合加速

基于CPU-GPU的混合計(jì)算系統(tǒng)是在現(xiàn)有CPU并行系統(tǒng)基礎(chǔ)上增加了GPU作為輔助并行加速部件，從而充分利用系統(tǒng)內(nèi)的可用資源，達(dá)到一種綜合加速的效果。本文所述的基于CPU-GPU混合加速的SPH流體仿真框架如圖1所示。該混合模擬仿真方法主要包括GPU流體計(jì)算和CPU流體渲染兩部分。

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

4.1 GPU流體計(jì)算

SPH流體模擬方法采用離散的粒子來(lái)模擬連續(xù)的流體運(yùn)動(dòng)，為了得到具有豐富細(xì)節(jié)的流體，需有大數(shù)量級(jí)粒子參與模擬計(jì)算。而采用純CPU計(jì)算方法進(jìn)行流體迭代，不能達(dá)到實(shí)時(shí)的模擬效果。由于SPH模擬方法本身具有高度的并行性，它非常適合用GPU進(jìn)行加速計(jì)算。借助GPU的高度并行性和可編程性，本文便能設(shè)計(jì)出在GPU上執(zhí)行的SPH流體模擬并行算法，從而實(shí)現(xiàn)大數(shù)量級(jí)粒子的流體系統(tǒng)的逼真、實(shí)時(shí)仿真。

4.1.1 構(gòu)建鄰居粒子表

SPH流體的模擬是通過(guò)離散的粒子實(shí)現(xiàn)的，離散粒子間的相互作用通過(guò)一個(gè)對(duì)稱的、具有一定作用半徑的光滑核函數(shù)計(jì)算得到。粒子間的相互作用具有一定作用半徑，對(duì)粒子的屬性進(jìn)行更新時(shí)，需遍歷系統(tǒng)中的其他所有粒子，這會(huì)導(dǎo)致非常低的查找效率。為此，Amada T等[11]提出構(gòu)造鄰居粒子表的方法，采用該方法計(jì)算粒子物理量時(shí)，可只查找鄰居粒子來(lái)進(jìn)行計(jì)算。其優(yōu)點(diǎn)在于每個(gè)迭代步，只需構(gòu)建一次鄰居粒子表，在粒子物理量計(jì)算過(guò)程中，都是查找表中數(shù)據(jù)進(jìn)行計(jì)算。但是，采用這種方法構(gòu)建鄰居粒子表時(shí)，采用的方法仍是遍歷所有粒子，導(dǎo)致建表時(shí)間過(guò)長(zhǎng)，影響模擬的實(shí)時(shí)性。后來(lái)一些學(xué)者在原有鄰居粒子算法基礎(chǔ)上，進(jìn)行了CPU端的優(yōu)化工作，如基于粒子對(duì)的鄰居粒子表構(gòu)建，這種改進(jìn)能使表的構(gòu)建降低一半工作量。Harada T等[12]提出了構(gòu)造空間均勻網(wǎng)格覆蓋整個(gè)粒子位置空間的算法，設(shè)置均勻網(wǎng)格的半徑與光滑核函數(shù)的作用半徑相同，這樣，在進(jìn)行粒子鄰居搜索時(shí)，只需要對(duì)當(dāng)前粒子所在網(wǎng)格及周圍26個(gè)網(wǎng)格進(jìn)行搜索即可。這種改進(jìn)極大提高了粒子的鄰居表構(gòu)建速度。

傳統(tǒng)的粒子對(duì)鄰居粒子表構(gòu)建方法是基于CPU的，其建表過(guò)程如算法1所示：

算法1基于CPU的鄰居粒子表構(gòu)建

遍歷每個(gè)粒子:

(1)將當(dāng)前粒子i增加到鄰居粒子表的i行第一個(gè)位置，并確定當(dāng)前粒子i所在的網(wǎng)格g;

(2)遍歷網(wǎng)格g及周圍網(wǎng)格中的粒子j;

① 如果j不在i的影響域內(nèi)，則遍歷下一粒子;

② 如果j在i的影響域內(nèi)，則i、j互為鄰居關(guān)系，將j添加到i的鄰居粒子表中。

(3)將當(dāng)前粒子i添加到網(wǎng)格g中，確定粒子i與網(wǎng)格g的歸屬關(guān)系。

結(jié)束遍歷。

該方法將確定粒子所在網(wǎng)格和搜索粒子鄰居兩項(xiàng)操作同步交錯(cuò)執(zhí)行。如果直接移植到GPU運(yùn)行，這種交錯(cuò)執(zhí)行方式會(huì)涉及到線程同步及共享資源加鎖問(wèn)題，這些操作會(huì)極大地降低GPU加速效果，且該方法構(gòu)造過(guò)程采用了鄰居對(duì)方法，鄰居粒子表最終存放的是只是單向的鄰居關(guān)系。這種單向的鄰居關(guān)系在后續(xù)粒子物理屬性求解階段會(huì)涉及資源加鎖問(wèn)題。

為此，本文設(shè)計(jì)出一種新穎的鄰居粒子表構(gòu)建算法，它將確定粒子所在網(wǎng)格與建立鄰居粒子表兩項(xiàng)操作分開(kāi)串行執(zhí)行，其中每項(xiàng)操作是在GPU中并行執(zhí)行。這種處理策略易于兩項(xiàng)操作的并行化并且所建立的鄰居粒子表中存放的雙向鄰居信息，能對(duì)后續(xù)粒子物理屬性進(jìn)行并行求解，從而能夠避免線程爭(zhēng)奪資源?；贕PU的建立鄰居粒子表算法如算法2所示。

Figure 2 Schematic diagram for SPH fluid simulation to optimize data transmission based on GPU acceleration圖2 基于GPU加速的SPH流體模擬數(shù)據(jù)傳輸優(yōu)化示意圖

算法2基于GPU的鄰居粒子表構(gòu)建算法

操作1 確定粒子的網(wǎng)格歸屬:

1：for 粒子表中粒子pido

2： 確定粒子pi所在的網(wǎng)格g;

3： 向網(wǎng)格g中添加粒子pi;

4： 網(wǎng)格g容量增1;

5：end for

操作2 建立鄰居粒子表:

6：for粒子表中粒子pido

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

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

9： 確定粒子pi所在的網(wǎng)格g;

10： for 網(wǎng)格g周圍網(wǎng)格g_neighbourdo

11： for 網(wǎng)格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 計(jì)算粒子屬性

SPH流體模擬方法中粒子屬性計(jì)算過(guò)程具有高度并行性，適合采用GPU并行計(jì)算。但是，該方法中屬性間可能有依賴關(guān)系，如粒子的加速度屬性是需要根據(jù)計(jì)算得到的密度和壓強(qiáng)屬性求得，這種依賴關(guān)系表明粒子的某些依賴屬性間是不能并行計(jì)算的。對(duì)于沒(méi)有依賴關(guān)系的粒子屬性，其計(jì)算可以放到一個(gè)核函數(shù)中進(jìn)行，以減少不必要的GPU調(diào)用開(kāi)銷；對(duì)于有依賴關(guān)系的粒子屬性，其計(jì)算需放到不同的有先后關(guān)系的核函數(shù)中并行計(jì)算。

4.1.3 數(shù)據(jù)傳輸優(yōu)化

擁有數(shù)十至數(shù)百個(gè)核心的眾核GPU與幾個(gè)核心的CPU相比，GPU能更快地進(jìn)行SPH流體模擬計(jì)算且具有高效的顯存讀寫能力。但是，一般基于GPU的并行計(jì)算程序中，CPU與GPU間的數(shù)據(jù)傳輸是并行能力提高的瓶頸。因此，本文在設(shè)計(jì)基于GPU加速的SPH流體模擬方法時(shí)，對(duì)CPU與GPU間的數(shù)據(jù)傳輸進(jìn)行了優(yōu)化處理。

基于GPU加速的SPH流體模擬數(shù)據(jù)傳輸優(yōu)化如圖2所示，圖2中小矩形框表示SPH方法的粒子物理屬性，橢圓形框表示SPH方法的主要計(jì)算工作。外層黑線框內(nèi)表示數(shù)據(jù)流在GPU內(nèi)部的流動(dòng)過(guò)程，它不涉及CPU與GPU間的數(shù)據(jù)交互；在粒子的各種物理量中，最終需要讀出的只有處于外層黑色框外的粒子的位置信息，但其他的速度、加速度、密度、壓強(qiáng)等物理量信息不必讀出而一直保存于GPU中，直接進(jìn)行下一步計(jì)算，從而大大減少了CPU與GPU的數(shù)據(jù)傳輸次數(shù)。

4.2 CPU加速渲染

為了逼真地展現(xiàn)模擬流體的運(yùn)動(dòng)，渲染與可視化顯得尤為重要，若采用球體對(duì)粒子系統(tǒng)渲染，在設(shè)定的光照條件下，便能清晰地展現(xiàn)出粒子與粒子間、粒子與容器邊界間的碰撞過(guò)程。在現(xiàn)有的基于GPU加速的SPH計(jì)算方法中，通常是對(duì)流體的計(jì)算過(guò)程進(jìn)行GPU并行加速，但渲染過(guò)程采用CPU串行執(zhí)行。本文所述SPH流體仿真方法中，利用多核CPU的硬件條件并結(jié)合OpenMP實(shí)現(xiàn)了CPU加速渲染。

OpenMP采用了Fork-Join并行執(zhí)行模式，如圖3所示。Fork-Join并行執(zhí)行模式的基本思想是，運(yùn)行時(shí)庫(kù)維護(hù)一個(gè)線程池，當(dāng)主線程遇到并行結(jié)構(gòu)時(shí)，從線程池創(chuàng)建一個(gè)并行的從線程組；當(dāng)從線程組完成并行區(qū)域的工作后，各從線程間進(jìn)行同步，從屬線程返回到池中，主線程繼續(xù)執(zhí)行并行構(gòu)造后的操作。

Figure 3 Fork-Join parallel execution model of OpenMP圖3 OpenMP Fork-Join并行執(zhí)行模型

5 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)對(duì)64 k和128 k粒子數(shù)的球體渲染操作在CPU串行和OpenMP多核CPU并行環(huán)境下進(jìn)行單幀平均耗時(shí)測(cè)試，結(jié)果對(duì)比情況如圖4所示(實(shí)驗(yàn)平臺(tái)的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加速前后單幀球體渲染時(shí)間對(duì)比

在雙核CPU和帶GPU顯卡的微機(jī)實(shí)驗(yàn)平臺(tái)(CPU：Intel? CoreTM2 Duo E7500 雙核2.93 GHz；GPU顯卡：ATI Radeon HD 5850,1 GB)上，通過(guò)對(duì)不同粒子數(shù)(128 k、64 k、27 k和8 k)SPH流體模擬方法中各部分操作進(jìn)行數(shù)據(jù)傳輸優(yōu)化得出如圖5所示的數(shù)據(jù)對(duì)比關(guān)系。圖5中四對(duì)兩條相鄰柱狀圖分別表示不同粒子數(shù)在單幀數(shù)據(jù)優(yōu)化前后的CPU與GPU間平均數(shù)據(jù)傳輸時(shí)間消耗。從圖5可以看出，通過(guò)對(duì)SPH流體模擬方法中各部分操作進(jìn)行數(shù)據(jù)傳輸優(yōu)化能很大程度上減小CPU-GPU間的數(shù)據(jù)交換消耗，從而加速SPH流體仿真的計(jì)算。

Figure 5 Comparison of transmission time between CPU and GPU after data optimization圖5 數(shù)據(jù)優(yōu)化前后CPU與GPU間單幀數(shù)據(jù)傳輸時(shí)間對(duì)比

粒子鄰居搜索消耗是SPH流體仿真方法的主要消耗之一，通過(guò)對(duì)粒子網(wǎng)格歸屬及鄰居粒子表構(gòu)建過(guò)程進(jìn)行改進(jìn)并經(jīng)GPU并行處理(實(shí)驗(yàn)平臺(tái)同圖5)，GPU并行前后粒子鄰居搜索消耗如表1所示。表1是一個(gè)迭代步構(gòu)建鄰居粒子表在純CPU端、GPU并行(未經(jīng)數(shù)據(jù)傳輸優(yōu)化)及GPU并行(經(jīng)數(shù)據(jù)傳輸優(yōu)化)的不同時(shí)間消耗對(duì)比，實(shí)驗(yàn)分別對(duì)8 k、27 k、128 k不同粒子數(shù)的單幀平均鄰居搜索時(shí)間進(jìn)行測(cè)試。

Table 1 Comparison of neighbor search time after parallel表1 并行前后鄰居搜索時(shí)間對(duì)比

從表1可知，經(jīng)數(shù)據(jù)傳輸優(yōu)化的GPU并行鄰居搜索效率比CPU串行鄰居搜索效率提高了近千倍。

實(shí)驗(yàn)對(duì)流體從圓柱形容器高處下落到容器中的情況進(jìn)行了SPH流體粒子模擬，模擬效果如圖6所示(實(shí)驗(yàn)平臺(tái)同圖5)。圖6a～圖6C分別是粒子數(shù)為8 k、27 k、128 k時(shí)的CPU-GPU混合加速流體模擬效果。

Figure 6 SPH fluid simulation based on CPU-GPU hybrid acceleration with different particles number圖6 CPU-GPU混合加速后不同粒子數(shù)流體模擬效果

表2給出了流體仿真中不同粒子數(shù)分別在CPU和CPU-GPU上進(jìn)行10 k個(gè)系統(tǒng)迭代的平均幀率對(duì)比情況(實(shí)驗(yàn)平臺(tái)同圖5)。從表2可以看出，即使粒子數(shù)達(dá)到64 k的大規(guī)模粒子場(chǎng)景，本文所述方法仍能夠達(dá)到30 FPS以上的實(shí)時(shí)模擬速度，與視頻采集的幀速率相當(dāng)。

Table 2 Comparison of frame rate with different particles number after parallel表2 不同粒子數(shù)并行前后幀率對(duì)比

為了更直觀、逼真地展現(xiàn)加速后SPH流體仿真的效果，本文對(duì)SPH粒子采用不同方式進(jìn)行渲染，圖7所示是粒子數(shù)為27 k時(shí)的加速后某一時(shí)刻渲染效果圖。圖7a～圖7c分別是采用點(diǎn)、小球和透明液體三種不同的渲染方式對(duì)SPH流體粒子進(jìn)行渲染。

Figure 7 Renderings of different render way with 27 000 particles圖7 粒子數(shù)為27 k時(shí)的三種不同渲染方式渲染效果圖

6 結(jié)束語(yǔ)

本文在介紹SPH流體模擬國(guó)內(nèi)外研究現(xiàn)狀的基礎(chǔ)上，提出一種基于CPU-GPU混合加速的SPH實(shí)時(shí)流體仿真方法，流體計(jì)算部分采用GPU并行加速,流體渲染部分采用基于CPU的OpenMP加速。通過(guò)前面的實(shí)驗(yàn)結(jié)果可知，運(yùn)用GPU的并行計(jì)算能力相比于純CPU實(shí)現(xiàn)而言，能實(shí)現(xiàn)更多粒子數(shù)的實(shí)時(shí)模擬；同時(shí)，利用基于多核CPU的OpenMP加速功能，能完成快速的渲染過(guò)程，從而實(shí)現(xiàn)逼真、實(shí)時(shí)的流體模擬仿真。

在未來(lái)的研究工作中，我們將研究CPU與GPU間的協(xié)同操作和負(fù)載均衡，以擴(kuò)大粒子的并行規(guī)模；研究高質(zhì)量的逼真實(shí)時(shí)流體渲染算法，并將本文所述方法應(yīng)用于各種實(shí)際的流體模擬工程領(lǐng)域中。

[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.