SPH-GPU并行計算在風沙流中的應用

梁嵐博，金阿芳，聞騰騰

新疆大學機械工程學院，烏魯木齊 830047

風沙運動為主要標志的荒漠化及其帶來的風沙環境是全球最突出的環境問題之一。風沙流活動對農業、牧業、沙漠公路交通、沙漠石油開采等造成了極大的危害，沙害造成的經濟損失高達數百億元人民幣。而風沙兩相流是一門既古老又年輕的科學，是認識風沙運動規律的重要基礎，同時也是開展沙漠化治理和促進沙漠化區域可持續發展的科學依據，因此風沙兩相流的研究具有重要的理論和實踐意義[1]。

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

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

近年來，隨著科學技術的迅速發展，科學研究領域對計算機性能的要求也在不斷增長，更是對計算能力提出了極高的要求[9]。傳統計算機的單核處理器CPU，已經遠遠不能滿足計算要求，即使是目前廣泛使用的多核處理器CPU，在龐大的計算量壓力下，也捉襟見肘[10]。為了滿足巨大計算需求，眾核架構圖形處理器（GPU）的出現成為了必然。目前，由于GPU 超強的計算能力且適合計算密集型和易于并行的程序，以廣泛應用到SPH 水體仿真算法中[11-12]，但應用到風沙運動中還比較匱乏。

因此本文引入了GPU 眾核架構，并且通過對SPH串行算法的改進，提出了一種SPH-GPU并行計算方法，用以在計算機上快速的模擬風沙運動。

1 風沙流的SPH描述和離散

1.1 SPH基本思想

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

1.2 SPH基本方程

SPH方法的構造按兩個關鍵步驟進行：第一步為積分表示法，又稱場函數核近似法；第二步為粒子近似法[15]。場函數核近似法如下式：

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

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

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

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

其中，N是粒子i的支持域內的粒子總數；xi、xj分別為i和j粒子的坐標向量；ρj、mj為j粒子的密度和質量。

1.3 N-S方程的SPH離散

1.3.1 連續密度法

對連續性方程進行速度散度SPH近似：

1.3.2 動量方程的粒子近似法

對動量方程的梯度項直接應用SPH 粒子近似法進行變換可得：

1.4 鄰近粒子的搜索

在模擬風沙運動的過程中，氣沙兩相粒子之間的位置是不斷變化的，所以需要在每個時間步內都要重新計算每個粒子作用域中所包含的粒子，可以說鄰近粒子搜索法的選擇是非常重要的一步。

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

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

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

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

2 計算條件

文中主要是基于氣沙兩相耦合之間的受力機制來模擬沙粒的運動過程。模擬過程中涉及的各類參數及運行環境分別為表1和表2所示。由于在邊界上的粒子積分時會被邊界截斷，會導致求解結果出錯，所以在計算區域上部和下部采用固體邊界，進出口采用周期性邊界。在計算區域頂部和沙粒的底部分別設置一層虛粒子，從而阻止粒子非物理穿透邊界如圖2所示。初始床面上的風速廓線公式：

圖2 計算區域初始狀態示意圖Fig.2 Schematic diagram of initial state of calculation area

表1 沙相、氣相物性和計算參數Tabel 1 Sand，gas phase physical properties and calculation parameters

表2 計算環境Tabel 2 Computing environment

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

3 SPH方法在CUDA計算架構下的實現

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

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

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

4 CPU-GPU并行計算結果與分析

4.1 SPH串行算法熱點分析

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

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

圖4 SPH串行熱點程序Fig.4 SPH serial hotspot program

4.2 不同搜索方法的串行SPH計算時間

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

表3 不同搜索法的SPH計算時間Tabel 3 SPH calculation time for different search methods

4.3 GPU并行計算模型驗證

4.3.1 沙粒群運動的時空變化

本算例重點討論流體起動，不考慮碰撞和沖擊產生的瞬時沖擊力的作用。

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

圖5 沙粒群運動的時空變化Fig.5 Temporal and spatial changes of movement of sand particles

4.3.2 沙粒的運動軌跡

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

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

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

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

4.4 GPU-CPU計算效率的對比

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

圖8 CPU-GPU計算時間Fig.8 CPU-GPU calculation time

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

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

5 結論

本文首先對串行SPH程序進行分析和修改，最后采用SPH-GPU大規模并行加速技術對風沙流動過程進行數值模擬，結果表明：

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

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

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