基于CUDA的并行沙粒模擬方法

吳蔚

摘 ?要：沙粒模擬是計(jì)算機(jī)圖形學(xué)中重要的研究?jī)?nèi)容。文章提出一種基于CUDA的并行化沙粒模擬方法，使用PCISPH流體模擬的算法實(shí)現(xiàn)沙粒流動(dòng)過(guò)程，并加入摩擦力模型保持沙粒穩(wěn)定堆積的效果。算法利用多線程進(jìn)行數(shù)值求解，充分挖掘了GPU并行計(jì)算的性能優(yōu)勢(shì)。實(shí)驗(yàn)表明，該算法可以實(shí)現(xiàn)高效的沙粒模擬，在動(dòng)畫、電影、工程等領(lǐng)域具有重要的應(yīng)用價(jià)值。

關(guān)鍵詞：沙粒模擬;SPH;并行計(jì)算;CUDA

中圖分類號(hào)：TP391.7 ? ? ? 文獻(xiàn)標(biāo)志碼：A ? ? ? ? ? ? ?文章編號(hào)：2095-2945（2020）25-0009-02

Abstract： Sand simulation is the key research content in computer graphics. A parallel sand simulation method based on CUDA is proposed， which use the PCISPH fluid simulation method to realize the free sand flow and adopt the friction model to keep in stable sandpile. The numerical solution is processed by multithread， which helps to fully excavate the performance advantage of GPU parallel computing. Experiments suggest that the method can realize the sand simulation efficiently， and bear important application value in the field of animation， movie， engineering， etc.

Keywords： sand simulation; SPH; parallel computing; CUDA

1 概述

沙粒等固體顆粒的模擬，是計(jì)算機(jī)圖形學(xué)中重要的研究?jī)?nèi)容，被廣泛應(yīng)用在動(dòng)畫、電影、工程等領(lǐng)域。在現(xiàn)實(shí)中可以觀察到，不穩(wěn)定的沙堆可以像流體一樣流動(dòng)。但與流體不同的是，沙粒內(nèi)部的能量可以迅速消散，當(dāng)這種能量達(dá)到一個(gè)較低值時(shí)，沙堆便會(huì)保持靜止的堆積狀態(tài)。之所以會(huì)表現(xiàn)出這樣復(fù)雜的特性，正是由于沙粒顆粒之間相互作用的結(jié)果。對(duì)于沙粒所表現(xiàn)出的流動(dòng)的效果，可以采用流體模擬的方法。而沙粒最終會(huì)趨于靜止，卻是由于沙粒之間的摩擦力抵消掉了重力的影響。因此沙粒模擬可以采用流體模擬的方法中加入粒子間摩擦力實(shí)現(xiàn)。

現(xiàn)今對(duì)于流體的模擬通常有兩類，即基于網(wǎng)格的歐拉方法和基于粒子的拉格朗日方法。歐拉方法將模擬區(qū)域劃分為網(wǎng)格，基于固定的網(wǎng)格點(diǎn)研究各時(shí)刻物理量在流場(chǎng)中的變化規(guī)律，其計(jì)算結(jié)果具有較高的數(shù)值精度;而拉格朗日方法則跟隨質(zhì)點(diǎn)運(yùn)動(dòng)，追蹤它們的物理量隨時(shí)間變化的規(guī)律，其更適用于需要表現(xiàn)表面細(xì)節(jié)或發(fā)生大變形的流體模擬應(yīng)用中。SPH（Smoothed Particle Hydrodynamics）方法是一種典型的基于拉格朗日視角的流體模擬方法，近年來(lái)成為最具前景的流體模擬方法之一。該方法將連續(xù)的流體離散成粒子，每個(gè)粒子攜帶一定的物理量，比如密度、壓強(qiáng)等，通過(guò)追蹤每個(gè)粒子的物理量變化從而表現(xiàn)整個(gè)流體的運(yùn)動(dòng)過(guò)程。在沙粒模擬中，單純使用SPH流體模擬方法并不能保持沙粒靜止堆積的狀態(tài)。因此沙粒粒子之間除了受到壓強(qiáng)力和速度差導(dǎo)致的黏性力等外力，還需要加入摩擦力的影響。Alduan等人[1]提出了一種基于SPH的沙粒模擬方法，使用結(jié)合單側(cè)不可壓縮性的PCISPH方法模擬沙粒的流動(dòng)，并加入摩擦力模型實(shí)現(xiàn)了沙粒的堆積效果。

通常來(lái)說(shuō)，基于粒子的數(shù)值計(jì)算方法比基于網(wǎng)格的方法更加高效，然而SPH方法中模擬介質(zhì)通常被離散成大量的粒子，每個(gè)粒子的物理量求解又依賴于對(duì)鄰域粒子的搜索以完成插值計(jì)算，采用傳統(tǒng)的CPU進(jìn)行迭代式地求解很難達(dá)到實(shí)時(shí)性的要求。因此有研究者提出使用GPU對(duì)SPH求解進(jìn)行加速的方法。比如Harada等[2]將SPH計(jì)算過(guò)程全部在GPU上實(shí)現(xiàn)，并使用桶紋理結(jié)構(gòu)實(shí)現(xiàn)了高效的鄰域搜索算法。Junior等[3]不僅實(shí)現(xiàn)了GPU上的數(shù)值求解，還使用空間哈希函數(shù)改進(jìn)了鄰域搜索算法，實(shí)現(xiàn)了無(wú)邊界的流體模擬。

本文基于CUDA架構(gòu)實(shí)現(xiàn)了并行化的沙粒模擬過(guò)程，將Alduan等人[1]的PCISPH數(shù)值求解和摩擦力模型完全基于GPU進(jìn)行并行化處理，避免了CPU和GPU之間不必要的數(shù)據(jù)傳輸，顯著提升了沙粒模擬的運(yùn)算效率。

2 基于CUDA架構(gòu)實(shí)現(xiàn)沙粒模擬

2.1 PCISPH方法

采用標(biāo)準(zhǔn)SPH方法模擬的沙粒具有可壓縮性，嚴(yán)重影響模擬的效果。本文使用PCISPH方法實(shí)現(xiàn)沙粒的不可壓縮模擬，每個(gè)時(shí)間步的算法步驟如下：

（1）首先針對(duì)每一個(gè)粒子搜索一定范圍內(nèi)的鄰域粒子，用于后續(xù)各種物理量的計(jì)算;

（2）使用SPH插值方法計(jì)算出所有非壓強(qiáng)力，比如粒子間粘性力、重力、所受外力等。同時(shí)將壓強(qiáng)和壓強(qiáng)力初始化為0;

（3）使用預(yù)測(cè)-矯正策略循環(huán)更新壓強(qiáng)和壓強(qiáng)力的值。迭代之初，使用粒子所受合力預(yù)測(cè)粒子臨時(shí)速度和臨時(shí)位置，并基于預(yù)測(cè)的位置信息預(yù)估出粒子的密度信息。該密度值用于矯正壓強(qiáng)和壓強(qiáng)力，并作為是否繼續(xù)下一迭代的依據(jù)。若當(dāng)前預(yù)測(cè)的密度值過(guò)大，說(shuō)明粒子之間過(guò)于緊密，需要按照如下公式增大粒子所受壓強(qiáng)Pi：

Pi+=δρerr

其中δ為預(yù)計(jì)算的常量，ρerr為與預(yù)測(cè)的密度相關(guān)的變量。所得壓強(qiáng)值用于計(jì)算出新的壓強(qiáng)力，重新進(jìn)入循環(huán)更新粒子的速度和位置。經(jīng)過(guò)若干次迭代，粒子間距便可滿足預(yù)先設(shè)定的壓縮率要求;

（4）基于粒子所有更新后的力進(jìn)行時(shí)間積分，得到下一時(shí)間步的粒子速度和位置。

上面步驟中第（3）步包含若干次迭代計(jì)算，是整個(gè)算法最耗時(shí)的部分。這種迭代式的方式使得每個(gè)時(shí)間步的計(jì)算比標(biāo)準(zhǔn)SPH方法花費(fèi)更多的時(shí)間。然而該方法允許設(shè)置較大的時(shí)間步長(zhǎng)，這反而極大地提升了沙粒模擬整體的運(yùn)行效率。

2.2 添加摩擦力模型

與單純的流體模擬不同，沙粒模擬還需要添加粒子間摩擦力，從而阻止粒子間的相對(duì)滑動(dòng)。摩擦力模型基于最大耗散原則實(shí)現(xiàn)，這又可通過(guò)最小化應(yīng)變速率來(lái)表達(dá)。為此，使用反推法基于粒子間相對(duì)速度得到粒子應(yīng)變速率值，然后迭代地計(jì)算出摩擦應(yīng)力張量值s，從而計(jì)算出粒子所受摩擦力大小：

其中θ為預(yù)定義的靜止角，用來(lái)決定沙粒穩(wěn)定堆積的坡度。由于該摩擦力模型本質(zhì)上也是基于預(yù)測(cè)-矯正策略，因此可以很容易地集成到前述PCISPH方法中。

2.3 基于CUDA的算法實(shí)現(xiàn)

SPH作為一種典型的拉格朗日方法，存在大量粒子間無(wú)關(guān)的數(shù)值計(jì)算過(guò)程，因此可以充分借助GPU的性能優(yōu)勢(shì)實(shí)現(xiàn)高度并行化。算法初始化時(shí)，所有沙粒粒子的屬性數(shù)據(jù)及邊界信息從CPU端內(nèi)存拷貝到GPU端顯存。后續(xù)的數(shù)值計(jì)算完全在GPU端進(jìn)行，避免了CPU和GPU之間不必要的數(shù)據(jù)傳輸。

沙粒模擬的每一幀都需要計(jì)算出所有粒子更新后的屬性值，比如密度、壓強(qiáng)、應(yīng)變速率等，以及各種受力，比如壓強(qiáng)力、摩擦力等。在CUDA編程模型中，可以并行執(zhí)行的代碼稱為核函數(shù)。因此針對(duì)沙粒模擬中的鄰域搜索及各種物理量的計(jì)算，本文使用核函數(shù)創(chuàng)建的輕量級(jí)線程與每個(gè)粒子一一對(duì)應(yīng)，實(shí)現(xiàn)對(duì)所有粒子的并行化處理。

2.4 實(shí)例驗(yàn)證

為了驗(yàn)證上述算法，設(shè)計(jì)了沙粒模擬實(shí)驗(yàn)，如圖1所示。初始時(shí)由沙粒顆粒組成的長(zhǎng)方體沙堆位于長(zhǎng)方體盒子包圍下的角落位置。相比于流體模擬，沙粒模擬由于粒子間摩擦力的作用，在流動(dòng)過(guò)程中逐漸呈現(xiàn)出明顯的堆積效果，驗(yàn)證了本文沙粒模型的準(zhǔn)確性。為了對(duì)本文基于CUDA的并行算法效果進(jìn)行驗(yàn)證，針對(duì)沙粒模擬場(chǎng)景分別采用CPU串行和GPU并行的方式實(shí)現(xiàn)，每幀平均運(yùn)行時(shí)間對(duì)比如表1所示。可以看到，相比于串行方式，并行化使得算法在運(yùn)行效率上達(dá)到了8倍的加速比，充分驗(yàn)證了本文并行化算法的有效性。

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

本文基于CUDA框架實(shí)現(xiàn)了并行化的沙粒模擬過(guò)程，將PCISPH流體模擬的方法和摩擦力模型完全基于GPU進(jìn)行數(shù)值求解，充分挖掘了GPU并行計(jì)算的性能優(yōu)勢(shì)，提升了沙粒模擬的運(yùn)行效率。今后的研究將在本文研究的基礎(chǔ)上繼續(xù)探索CUDA優(yōu)化技術(shù)在大規(guī)模沙粒模擬場(chǎng)景中的效果。

參考文獻(xiàn)：

[1]Alduán I， Otaduy M A. SPH granular flow with friction and cohesion[C]//Proceedings of the 2011 ACM SIGGRAPH/Eurographics symposium on computer animation， 2011： 25-32.

[2]Harada T ， Koshizuka S ， Kawaguchi Y . Smoothed Particle Hydrodynamics on GPUs[J]. Proc.computer Graphics Int.rio De Janeiro Brazil May Jun， 2007， 4（4）：671-691.

[3]Junior J R S， Joselli M， Zamith M， et al. An architecture for real time fluid simulation using multiple GPUs[J]. XI SBGames， Brasília， 2012：8.