SPH-GPU并行計(jì)算在風(fēng)沙流中的應(yīng)用

梁嵐博，金阿芳，聞騰騰

新疆大學(xué)機(jī)械工程學(xué)院，烏魯木齊 830047

風(fēng)沙運(yùn)動(dòng)為主要標(biāo)志的荒漠化及其帶來的風(fēng)沙環(huán)境是全球最突出的環(huán)境問題之一。風(fēng)沙流活動(dòng)對(duì)農(nóng)業(yè)、牧業(yè)、沙漠公路交通、沙漠石油開采等造成了極大的危害，沙害造成的經(jīng)濟(jì)損失高達(dá)數(shù)百億元人民幣。而風(fēng)沙兩相流是一門既古老又年輕的科學(xué)，是認(rèn)識(shí)風(fēng)沙運(yùn)動(dòng)規(guī)律的重要基礎(chǔ)，同時(shí)也是開展沙漠化治理和促進(jìn)沙漠化區(qū)域可持續(xù)發(fā)展的科學(xué)依據(jù)，因此風(fēng)沙兩相流的研究具有重要的理論和實(shí)踐意義[1]。

風(fēng)沙兩相流的研究方法主要有風(fēng)洞模擬、野外觀測(cè)和數(shù)值模擬。風(fēng)洞試驗(yàn)觀測(cè)是研究風(fēng)沙運(yùn)動(dòng)最常用的手段，通常采用非接觸測(cè)量?jī)x器如粒子圖像測(cè)速儀（particle image velocimeter，PIV）、相位多普勒粒子分析儀（phase Doppler particle analyzer，PDPA）、高速攝影儀等，可以對(duì)沙粒的微細(xì)運(yùn)動(dòng)和風(fēng)洞輸沙的時(shí)空變化進(jìn)行分析。野外觀測(cè)更多采用風(fēng)速儀和集沙儀以及各種類型的輸沙傳感器等對(duì)風(fēng)沙運(yùn)動(dòng)及輸沙量進(jìn)行測(cè)量[2]。隨著計(jì)算機(jī)硬件和軟件以及計(jì)算方法的發(fā)展，數(shù)值模擬近年來在風(fēng)沙兩相流研究中的作用越來越重要。數(shù)值模擬方法便于對(duì)風(fēng)沙運(yùn)動(dòng)進(jìn)行定量模擬和預(yù)測(cè)，易于控制參數(shù)，相對(duì)于風(fēng)洞實(shí)驗(yàn)和野外觀測(cè)有成本較低和便于實(shí)現(xiàn)的優(yōu)點(diǎn)。

目前，流體數(shù)值模擬研究包括基于網(wǎng)格的歐拉法和基于粒子的拉格朗日法。在基于歐拉網(wǎng)格的數(shù)值方法中如FDM，首先需要在問題域上生成網(wǎng)格，要在不規(guī)則或者很復(fù)雜的幾何形狀上構(gòu)造規(guī)律的網(wǎng)格是一件非常困難的事情，在流體模擬的整個(gè)過程中還會(huì)出現(xiàn)很多如大變形、變形邊界、自由表面和物質(zhì)交界面的問題，給基于網(wǎng)格的數(shù)值模擬方法帶來了很大的挑戰(zhàn)。拉格朗日無網(wǎng)格法主要有SPH、MPS和DEM等方法，以上粒子之間不需要進(jìn)行網(wǎng)格劃分（不是預(yù)先定義好的網(wǎng)格系統(tǒng)），是一系列任意分布的離散點(diǎn)，因此可以克服流體模擬中網(wǎng)格畸變和網(wǎng)格移動(dòng)等引起的各種問題。目前，SPH方法被廣泛應(yīng)用于各種領(lǐng)域如潰壩、輪胎劃水、飛機(jī)開溝、血流血管和船的撞擊等[3-4]，且該方法更加穩(wěn)定和真實(shí)，但SPH方法運(yùn)用在風(fēng)沙兩相流中還較為缺乏。文獻(xiàn)[5]首次將SPH 方法運(yùn)用在風(fēng)沙流中，建立了基于SPH 方法的風(fēng)沙流控制方程，提出了模擬過程中一些關(guān)鍵參數(shù)的選擇，并使用虛粒子法阻止粒子穿透邊界。文獻(xiàn)[6]利用文獻(xiàn)[5]的算法，對(duì)沙粒的躍移運(yùn)動(dòng)進(jìn)行了模擬，并在模擬過程中引入了起風(fēng)和停風(fēng)模塊，分析了在風(fēng)的加載和卸載時(shí)，沙粒的運(yùn)動(dòng)軌跡。文獻(xiàn)[7]利用文獻(xiàn)[5]的算法，建立了具有坡面形態(tài)的計(jì)算模型，得到了坡面不同位置沙粒的躍移高度和水平距離，驗(yàn)證了沙粒躍移的平均距離是單個(gè)沙波紋長(zhǎng)度的1.16倍。文獻(xiàn)[8]針對(duì)歐拉-歐拉流體計(jì)算模型在求解氣沙兩項(xiàng)耦合流動(dòng)問題中存在的缺陷，提出了基于SPH-FVM 耦合的方法模擬風(fēng)沙運(yùn)動(dòng)，用SPH 法離散沙粒相進(jìn)行求解，用有限體積法（finite volume method，F(xiàn)VM）求解氣相粒子，氣沙兩相通過拖曳力以及體積分?jǐn)?shù)等參數(shù)進(jìn)行耦合，模擬了自由來流風(fēng)作用下沙粒的運(yùn)動(dòng)過程以及沙丘緩慢向前運(yùn)動(dòng)的過程。由于以上模擬過程中需要將計(jì)算區(qū)域離散成數(shù)量龐大的單個(gè)粒子，其計(jì)算時(shí)間長(zhǎng)和模擬效率低一直是該方法的瓶頸。

近年來，隨著科學(xué)技術(shù)的迅速發(fā)展，科學(xué)研究領(lǐng)域?qū)τ?jì)算機(jī)性能的要求也在不斷增長(zhǎng)，更是對(duì)計(jì)算能力提出了極高的要求[9]。傳統(tǒng)計(jì)算機(jī)的單核處理器CPU，已經(jīng)遠(yuǎn)遠(yuǎn)不能滿足計(jì)算要求，即使是目前廣泛使用的多核處理器CPU，在龐大的計(jì)算量壓力下，也捉襟見肘[10]。為了滿足巨大計(jì)算需求，眾核架構(gòu)圖形處理器（GPU）的出現(xiàn)成為了必然。目前，由于GPU 超強(qiáng)的計(jì)算能力且適合計(jì)算密集型和易于并行的程序，以廣泛應(yīng)用到SPH 水體仿真算法中[11-12]，但應(yīng)用到風(fēng)沙運(yùn)動(dòng)中還比較匱乏。

因此本文引入了GPU 眾核架構(gòu)，并且通過對(duì)SPH串行算法的改進(jìn)，提出了一種SPH-GPU并行計(jì)算方法，用以在計(jì)算機(jī)上快速的模擬風(fēng)沙運(yùn)動(dòng)。

1 風(fēng)沙流的SPH描述和離散

1.1 SPH基本思想

SPH 方法是在1977 年由Gingold 和Monaghan 等人[13-14]提出，最初用于解決天體物理學(xué)問題，后來經(jīng)過其他學(xué)者的不斷改進(jìn)和完善，推廣到流體力學(xué)各個(gè)方面。該無網(wǎng)格方法是將連續(xù)的流體離散成一系列任意分布的SPH 粒子，每個(gè)粒子攜帶質(zhì)量、速度、密度等物理信息，通過求解SPH控制方程來描述整個(gè)系統(tǒng)的演變過程。

1.2 SPH基本方程

SPH方法的構(gòu)造按兩個(gè)關(guān)鍵步驟進(jìn)行：第一步為積分表示法，又稱場(chǎng)函數(shù)核近似法；第二步為粒子近似法[15]。場(chǎng)函數(shù)核近似法如下式：

式中，Ω是支持域，f(x)為三維坐標(biāo)向量x的函數(shù)，δ(x-x′)是狄拉克函數(shù)，性質(zhì)如下：

用光滑函數(shù)W(x-x′)取代f(x)積分表示式中的δ(x-x′)則：

式中，角括號(hào)＜＞在SPH的習(xí)慣用法中是核近似算子的標(biāo)記；W(x-x′)即光滑函數(shù)；h是光滑長(zhǎng)度，用來定義光滑函數(shù)影響區(qū)域的范圍。

式（3）可轉(zhuǎn)化為支持域內(nèi)所有粒子疊加求和的離散化形式，稱為粒子近似法。在粒子i處的粒子近似形式可寫為：

其中，N是粒子i的支持域內(nèi)的粒子總數(shù)；xi、xj分別為i和j粒子的坐標(biāo)向量；ρj、mj為j粒子的密度和質(zhì)量。

1.3 N-S方程的SPH離散

1.3.1 連續(xù)密度法

對(duì)連續(xù)性方程進(jìn)行速度散度SPH近似：

1.3.2 動(dòng)量方程的粒子近似法

對(duì)動(dòng)量方程的梯度項(xiàng)直接應(yīng)用SPH 粒子近似法進(jìn)行變換可得：

1.4 鄰近粒子的搜索

在模擬風(fēng)沙運(yùn)動(dòng)的過程中，氣沙兩相粒子之間的位置是不斷變化的，所以需要在每個(gè)時(shí)間步內(nèi)都要重新計(jì)算每個(gè)粒子作用域中所包含的粒子，可以說鄰近粒子搜索法的選擇是非常重要的一步。

圖1（a）全配對(duì)搜索法是最簡(jiǎn)單和直接的搜索方法。對(duì)于目標(biāo)粒子i（實(shí)心紅色圓點(diǎn)），要計(jì)算它到問題域中每一個(gè)粒子j(j=1,2,…,N,N是問題域中的粒子總數(shù)，圖中用藍(lán)色實(shí)心圓點(diǎn)表示）之間的距離，若距離小于粒子的支持域的半徑kh，則粒子j為粒子i的鄰近粒子（所有在紅色圓圈之內(nèi)的點(diǎn)）。在本文中因考慮對(duì)稱光滑長(zhǎng)度，則粒子i和粒子j為一組相互作用的粒子對(duì)，即粒子i也在粒子j的支持域內(nèi)，那么，在搜索粒子j的鄰近粒子時(shí)就不需要再考慮粒子i，因此可以減少一半的搜索計(jì)算量。該搜索方法對(duì)計(jì)算域內(nèi)的所有粒子進(jìn)行，雖然這個(gè)方法的應(yīng)用較簡(jiǎn)單，但顯然，這種方法的時(shí)間復(fù)雜度是O(N2)。

圖1 鄰近粒子的搜索Fig.1 Search for nearby particles

在運(yùn)用鏈表搜索法時(shí)，在問題域上要鋪設(shè)一臨時(shí)網(wǎng)格，如圖1（b）所示。將每個(gè)粒子都分布在網(wǎng)格單元內(nèi)，并通過簡(jiǎn)單的存儲(chǔ)規(guī)則將每個(gè)網(wǎng)格內(nèi)的所有粒子連接起來。網(wǎng)格單元的大小與粒子支持域的大小相同，此網(wǎng)格只是用來劃分區(qū)域，不參與計(jì)算。即，若粒子i的支持域大小為kh，則網(wǎng)格單元的大小也為kh。那么，對(duì)于給定的粒子i，其鄰近粒子只能在同一網(wǎng)格或相鄰網(wǎng)格里，由于本文的問題域是二維的，所以相鄰網(wǎng)格的數(shù)量為8。若每個(gè)單元內(nèi)的平均粒子數(shù)量很小，則鏈表搜索法的復(fù)雜度為O(N)。

上述算法在構(gòu)建相鄰粒子對(duì)時(shí)所使用到的輸入數(shù)據(jù)僅與當(dāng)前粒子（沙粒、氣體以及虛粒子）的空間位置有關(guān)，與其他變量無關(guān)。每次計(jì)算輸出的結(jié)果是粒子對(duì)的信息，并不改變當(dāng)前時(shí)刻的粒子位置數(shù)據(jù)。在構(gòu)建粒子對(duì)循環(huán)函數(shù)體中計(jì)算出的結(jié)果不影響輸入數(shù)據(jù)，那么函數(shù)體中每次迭代的過程都不相互依賴。因此任一時(shí)刻某個(gè)粒子在尋找相鄰粒子對(duì)時(shí)都是相互獨(dú)立的，也就是說讓GPU中的每一個(gè)線程攜帶粒子空間位置信息進(jìn)行計(jì)算，每個(gè)線程執(zhí)行時(shí)的區(qū)別僅是輸入位置數(shù)據(jù)的不同，而線程的計(jì)算指令保持一致，符合GPU大規(guī)模并行計(jì)算的要求。

2 計(jì)算條件

文中主要是基于氣沙兩相耦合之間的受力機(jī)制來模擬沙粒的運(yùn)動(dòng)過程。模擬過程中涉及的各類參數(shù)及運(yùn)行環(huán)境分別為表1和表2所示。由于在邊界上的粒子積分時(shí)會(huì)被邊界截?cái)啵瑫?huì)導(dǎo)致求解結(jié)果出錯(cuò)，所以在計(jì)算區(qū)域上部和下部采用固體邊界，進(jìn)出口采用周期性邊界。在計(jì)算區(qū)域頂部和沙粒的底部分別設(shè)置一層虛粒子，從而阻止粒子非物理穿透邊界如圖2所示。初始床面上的風(fēng)速廓線公式：

圖2 計(jì)算區(qū)域初始狀態(tài)示意圖Fig.2 Schematic diagram of initial state of calculation area

表1 沙相、氣相物性和計(jì)算參數(shù)Tabel 1 Sand，gas phase physical properties and calculation parameters

表2 計(jì)算環(huán)境Tabel 2 Computing environment

式中，μz為床面高程z處的風(fēng)速；z0為沙床面粗糙度，z0=Ds/30。

3 SPH方法在CUDA計(jì)算架構(gòu)下的實(shí)現(xiàn)

為了實(shí)現(xiàn)并行計(jì)算，需要對(duì)SPH 串行程序進(jìn)行分析和修改，首先尋找程序中的熱點(diǎn)函數(shù)（程序中耗時(shí)最多的函數(shù)），此外還要考慮數(shù)據(jù)的依賴性（前面的計(jì)算結(jié)果緊接著被后面的函數(shù)用到）分析哪部分程序適合串行和并行計(jì)算等等。最后用C語言編寫主流程，采用CUDA C 語言調(diào)用GPU 設(shè)備實(shí)現(xiàn)并行計(jì)算，使改進(jìn)后的程序，能使CPU與GPU協(xié)同合作完成任務(wù)，最大限度的發(fā)揮計(jì)算機(jī)的并行處理能力，從而提高SPH數(shù)值計(jì)算效率。

如圖3 所示模擬過程分為4 個(gè)部分：讀取需要模擬的數(shù)據(jù)、搜尋鄰居粒子、計(jì)算粒子間作用力以及粒子位置的更新。CUDA 程序開始計(jì)算時(shí)，CPU 程序占主導(dǎo)地位，主機(jī)端創(chuàng)建輸入和輸出數(shù)組，為輸入數(shù)據(jù)和結(jié)果提供存儲(chǔ)空間，將已有的粒子信息文件（初始的位置、速度、密度等屬性值）讀取到CPU 內(nèi)存里；隨后CPU加載cudaMalloc函數(shù)分配GPU顯存和cudaMemcpy函數(shù)把主機(jī)內(nèi)存里的粒子信息拷貝至GPU 顯存里，接下來CPU 加載核函數(shù)給GPU 做計(jì)算，GPU 計(jì)算時(shí)CPU阻塞不執(zhí)行其他任務(wù)，核函數(shù)在設(shè)備端進(jìn)行相鄰粒子搜索，搜索之后計(jì)算粒子屬性值，在計(jì)算粒子屬性值時(shí)要考慮數(shù)據(jù)之間的依賴性，具有依賴性的數(shù)據(jù)不適合進(jìn)行并行計(jì)算，否則會(huì)導(dǎo)致計(jì)算結(jié)果不準(zhǔn)確甚至程序崩潰；最后進(jìn)行粒子位置更新，更新完成后，跳至鄰居粒子搜索進(jìn)行循環(huán)執(zhí)行，在完成所有計(jì)算任務(wù)之后，才將所得結(jié)果拷貝到內(nèi)存中，這樣就減少了內(nèi)存和顯存數(shù)據(jù)傳輸所花費(fèi)的時(shí)間，提高了并行計(jì)算的效率。這個(gè)模擬過程中還要考慮合理利用設(shè)備資源，比如使SM 始終保持計(jì)算狀態(tài)，合理的分配線程，選擇合適的存儲(chǔ)器等。

圖3 GPU加速的SPH計(jì)算流程Fig.3 Accelerated SPH calculation process through GPU

4 CPU-GPU并行計(jì)算結(jié)果與分析

4.1 SPH串行算法熱點(diǎn)分析

并行算法實(shí)現(xiàn)之前，要對(duì)整個(gè)串行程序進(jìn)行熱點(diǎn)分析，在進(jìn)行并行加速時(shí)，要考慮串行程序中計(jì)算時(shí)間很長(zhǎng)的程序段，也就是熱點(diǎn)程序。如果首先對(duì)熱點(diǎn)程序進(jìn)行并行實(shí)現(xiàn)，那么對(duì)整個(gè)程序計(jì)算效率的提升是非常巨大的。上結(jié)所說模擬過程中耗時(shí)的部分主要有3個(gè)：搜尋鄰居粒子、計(jì)算粒子間作用力以及粒子更新。

如圖4是模擬25 000步得到的結(jié)果，從圖中可以很明顯地看出，整個(gè)程序耗時(shí)最大的部分在于搜尋鄰居粒子上，這部分占據(jù)了整個(gè)程序一半以上的時(shí)間，是串行程序的主要熱點(diǎn)所在。計(jì)算粒子間作用力，花費(fèi)了接近40%的時(shí)間，也是熱點(diǎn)程序的部分。相比較而言，粒子信息更新在總時(shí)間中所占據(jù)的比例僅僅2%。以上可得，在進(jìn)行并行算法實(shí)現(xiàn)時(shí)，應(yīng)該首先分析搜尋鄰居粒子這段程序上，尋找更好的粒子搜尋方法，其次在優(yōu)化計(jì)算粒子間相互作用力的程序，最后對(duì)于粒子更新可以選擇保留串行計(jì)算，也可以進(jìn)行并行計(jì)算。

圖4 SPH串行熱點(diǎn)程序Fig.4 SPH serial hotspot program

4.2 不同搜索方法的串行SPH計(jì)算時(shí)間

由上一節(jié)得出在搜尋鄰近粒子上的耗時(shí)占用了整個(gè)程序的大部分時(shí)間，因此表3比較了全配對(duì)法和鏈表法所占的耗時(shí)，由于時(shí)間復(fù)雜度由O(N2)降低到O(N)，可見鏈表法是全配對(duì)搜索法計(jì)算速度的2.3倍。

表3 不同搜索法的SPH計(jì)算時(shí)間Tabel 3 SPH calculation time for different search methods

4.3 GPU并行計(jì)算模型驗(yàn)證

4.3.1 沙粒群運(yùn)動(dòng)的時(shí)空變化

本算例重點(diǎn)討論流體起動(dòng)，不考慮碰撞和沖擊產(chǎn)生的瞬時(shí)沖擊力的作用。

如圖5 給出了在SPH-GPU 并行計(jì)算下時(shí)間步為1 500步圖（a），5 100步圖（b），8 400步圖（c），9 900步圖（d），19 200 步圖（e），21 600 步圖（f）的計(jì)算結(jié)果。從圖中可知，沙粒的起動(dòng)是在重力、氣流拖曳力等各種力共同作用的結(jié)果，計(jì)算結(jié)果較好地描述了沙粒群的運(yùn)動(dòng)規(guī)律。如圖（a）、（b）所示，在氣流的作用下，平坦沙床上層沙粒產(chǎn)生加速度，當(dāng)風(fēng)速足夠大時(shí)，沙粒從地面跳起開始運(yùn)動(dòng)。如圖（c）、（d）中橢圓區(qū)域所示，隨著時(shí)間的推移，當(dāng)時(shí)間步進(jìn)行到8 400步時(shí)，沙床此時(shí)開始發(fā)生松動(dòng)或堆積，沙床中有些沙粒會(huì)發(fā)生離開沙床起躍的趨勢(shì)，也有些沙粒其速度朝下，這些粒子會(huì)造成沙床更加緊實(shí)，沙床的松動(dòng)和緊實(shí)是造成沙床產(chǎn)生沙波紋的直接原因之一。如圖（e）、（f）所示，起跳沙粒在氣流中不斷的加速，又在重力作用下有了明顯的下落。當(dāng)沙床面上起跳的沙粒數(shù)與落回床面的沙粒數(shù)相差不大時(shí)，風(fēng)沙運(yùn)動(dòng)進(jìn)入了比較穩(wěn)定的狀態(tài)。

圖5 沙粒群運(yùn)動(dòng)的時(shí)空變化Fig.5 Temporal and spatial changes of movement of sand particles

4.3.2 沙粒的運(yùn)動(dòng)軌跡

如圖6為床面上典型沙粒的躍移運(yùn)動(dòng)軌跡，從圖中可以看到，沙粒的運(yùn)動(dòng)軌跡具有上升段比較急劇，降落段比較平緩的特點(diǎn)，都呈現(xiàn)出一種非對(duì)稱的拋物線形狀，觀察到靠近右周期邊界的躍移軌跡的下落段并沒有顯示完整，這是由于進(jìn)出口采用周期性邊界，右邊界出去的粒子從左邊界進(jìn)來，以上沙粒運(yùn)動(dòng)軌跡與文獻(xiàn)[16]風(fēng)洞實(shí)驗(yàn)中高速攝影記錄的沙粒躍移運(yùn)動(dòng)相符。分析這是由于沙粒在躍移過程中始終受到氣流拖曳力的作用，導(dǎo)致了在水平方向的分速度越來越大，當(dāng)沙粒飛到一定高度，又在氣流阻力和重力的作用下，迅速向下運(yùn)動(dòng)，因此導(dǎo)致降落段的水平位移較大。

圖6 典型沙粒躍移軌跡Fig.6 Typical sand grain trajectories

如圖7 表示沙粒運(yùn)動(dòng)中具有拋物線形狀但同時(shí)又有些變異的軌跡，其躍移軌跡并不是非常光滑的拋物線形狀，而是在轉(zhuǎn)折處出現(xiàn)尖角，即沙粒會(huì)急劇降落，造成此現(xiàn)象的原因可能是局部氣流縱向脈動(dòng)較大，也可能是沙粒在空中與其他沙粒發(fā)生碰撞后獲得了一個(gè)向下的沖量。

圖7 變異的尖角軌跡Fig.7 Mutated sharp trajectories

4.4 GPU-CPU計(jì)算效率的對(duì)比

如圖8 顯示了粒子數(shù)分別為2 500、10 000、40 000來計(jì)算不同粒子數(shù)下GPU與CPU的效率，可以看出，使用不同的鄰近粒子搜索法在GPU中的計(jì)算效率遠(yuǎn)高于在CPU 中的計(jì)算效率。GPU 并行鏈表搜索法相較于CPU 鏈表搜索法最高可以獲得3.1 倍的加速比，而相對(duì)于CPU全配對(duì)搜索法最高可以獲得20倍的加速比。當(dāng)粒子數(shù)目為2 500 時(shí)，GPU 并行計(jì)算的加速效果并不明顯，其原因是GPU的單線程計(jì)算能力遠(yuǎn)弱于CPU，隨著粒子數(shù)目的增加GPU多線程并行計(jì)算能力的優(yōu)勢(shì)就發(fā)揮了出來，能夠帶來更大的加速比。

圖8 CPU-GPU計(jì)算時(shí)間Fig.8 CPU-GPU calculation time

為了更好地理解CPU 與GPU 的加速機(jī)制，以10 000 粒子為例，統(tǒng)計(jì)了搜尋鄰居粒子、計(jì)算粒子間作用力以及粒子更新在計(jì)算時(shí)占總時(shí)間的比例。如圖9所示，GPU并行計(jì)算相較于CPU計(jì)算時(shí)，搜尋鄰居粒子和粒子更新所占時(shí)間比例都減少了，然而計(jì)算粒子間作用的時(shí)間比例增加了。綜合以上，可以推測(cè)出搜尋鄰居粒子和粒子更新相對(duì)于計(jì)算粒子間作用力有更好的并行計(jì)算效率。

圖9 SPH各部分占總時(shí)間的比例Fig.9 Proportion of each parts of SPH in total time

5 結(jié)論

本文首先對(duì)串行SPH程序進(jìn)行分析和修改，最后采用SPH-GPU大規(guī)模并行加速技術(shù)對(duì)風(fēng)沙流動(dòng)過程進(jìn)行數(shù)值模擬，結(jié)果表明：

（1）對(duì)SPH 串行算法代碼進(jìn)行了熱點(diǎn)程序分析，得到在搜尋鄰居粒子和計(jì)算粒子間作用力時(shí)耗時(shí)最長(zhǎng)，因此優(yōu)先優(yōu)化這兩部分程序段。由于搜尋鄰居粒子是最耗時(shí)的部分，比較了全配對(duì)搜索法和鏈表搜索法的計(jì)算效率，得到鏈表法是全配對(duì)搜索法計(jì)算速度的2.3 倍。

（2）對(duì)SPH-GPU 并行計(jì)算模型進(jìn)行驗(yàn)證。宏觀上得到了沙粒群在各種力的作用下有起跳、上升和回落三個(gè)階段的時(shí)空變化規(guī)律以及沙波紋形成過程，微觀上得到了典型的沙粒躍移軌跡和變異的尖角軌跡。

（3）比較了CPU 與GPU 的計(jì)算效率，得到GPU 大規(guī)模并行計(jì)算的執(zhí)行時(shí)間遠(yuǎn)遠(yuǎn)低于CPU 的執(zhí)行時(shí)間。在粒子數(shù)較少時(shí)GPU 加速效果并不大，隨著粒子數(shù)的增加，加速效果越加明顯，最高可以獲得20 倍的加速比。為了更好地理解CPU 與GPU 的加速機(jī)制，比較了搜尋鄰居粒子、計(jì)算粒子間作用力以及粒子更新在計(jì)算時(shí)占總時(shí)間的比例，推測(cè)出搜尋鄰居粒子和粒子更新相對(duì)于計(jì)算粒子間作用力有更好的并行計(jì)算效率。