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

吳蔚

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

關鍵詞：沙粒模擬;SPH;并行計算;CUDA

中圖分類號：TP391.7 ? ? ? 文獻標志碼：A ? ? ? ? ? ? ?文章編號：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 概述

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

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

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

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

2 基于CUDA架構實現(xiàn)沙粒模擬

2.1 PCISPH方法

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

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

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

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

Pi+=δρerr

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

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

上面步驟中第（3）步包含若干次迭代計算，是整個算法最耗時的部分。這種迭代式的方式使得每個時間步的計算比標準SPH方法花費更多的時間。然而該方法允許設置較大的時間步長，這反而極大地提升了沙粒模擬整體的運行效率。

2.2 添加摩擦力模型

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

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

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

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

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

2.4 實例驗證

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

3 結束語

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

參考文獻：

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