一種改進GPU加速策略在光滑粒子流體動力學方法中的應用

管延敏, 楊彩虹, 康 莊, 周 利

(1. 江蘇科技大學 船舶與海洋工程學院,江蘇 鎮江 212003;2. 哈爾濱工程大學 船舶工程學院,哈爾濱 150001)

光滑粒子流體動力學(Smoothed Particle Hydrodynamics,SPH)對模擬自由表面流動問題具有先天優勢,但隨著粒子數量的增加,尤其是三維流動問題,計算效率急劇降低,極大限制了該方法在大規模計算中的適用性.利用圖形處理器(GPU)眾核架構開展多核并行計算的計算機統一設備架構(CUDA)技術,其強大的并行特性非常適合解決大規模的高級計算問題[1],尤其是計算密集型問題.SPH計算能力密集的特性使其能容易地在GPU上實現并行運算,因此一些學者開始嘗試將GPU技術應用于SPH方法.Harada等[2]首次描述了SPH方法在GPU上的實現,Crespo等[3]將GPU加速的SPH方法應用于模擬復雜的自由表面流動.他們的研究均表明在單個GPU上的SPH模擬比在單核中央處理器(CPU)上進行的SPH模擬要快兩個數量級.目前SPH方法在GPU平臺下的計算研究發展迅速,GPU加速的SPH模型已應用于顆粒流體流動的模擬和淺水方程的求解[4-5],并通過結合自適應粒子技術[6]實現了高效計算.國內,徐鋒[7]在實現了基于GPU眾核架構的并行SPH算法的同時,結合所使用的GPU硬件上的特點,對并行算法進行了優化,將計算效率提高了20倍以上;金善勤等[8]提出一種基于粒子對的并行計算方法并將其與改進的SPH方法結合,實現了超過10的加速比;楊志國等[9]將GPU算法應用于二維潰壩模擬中也獲得了數量級的加速效果.

但是,GPU的系統設計理念與CPU正好相反,GPU面對的是類型一致、互不關聯的大量數據,其顯存的讀寫對算法求解效率影響很大.隨著流場的不斷演變,SPH粒子的無序化很容易導致多個線程同時讀寫同一地址而引起訪問沖突,極大地影響了計算效率的穩定性.針對此問題,并行索引排序方法[3, 10]、Z索引排序方法[11]及其改進方法[12]相繼被提出,這些方法通過優化SPH粒子索引存儲方式使同單元粒子在GPU顯存中盡可能相鄰,一定程度上改善了GPU顯存訪問的不連續問題,但沒有從實質上改變單元粒子的無序化.為此,本文提出了一種粒子重新編號技術,實現了SPH粒子的有序排列和GPU顯存的連續訪問,并將該方法應用于三維帶障礙物潰壩模擬,通過與實驗數據比較以及不同硬件設施上不同算法求解效率的對比,驗證了本文方法的精確性和高效性.

1 SPH的控制方程

SPH控制方程的離散形式可表示為

(1)

(2)

式中:ρi、pi、ui、Vj分別為流體粒子i的密度、壓力、速度和體積;t為時間;mj和ρj為粒子i支持域中粒子j的質量和密度;g為重力加速度;Wij為核函數,本文采用Wendland[13]的C2核函數:

(3)

(0≤q≤2)

式中:r為兩粒子間距離;h為光滑長度;q為粒子間相對距離,q=r/h;αd為常數,對于三維問題,αd=21/8πh3.

人工黏性為

(4)

式中:

(5)

假定流體弱可壓縮,采用Monaghan等[14]提出的人工壓縮法求解壓強:

(6)

2 CUDA-GPU架構下的SPH模型

2.1 并行算法的實現流程

針對現代CUDA-GPU架構設計了一種高效的SPH流體數據,以提高SPH流體模擬在GPU上的數據處理速度.并行代碼使用C#語言和CUDA編程語言進行開發,大多數源代碼都是CPU和GPU所共有的,可以在CPU或GPU上執行,也可以在沒有啟用GPU的工作站上運行,只使用CPU實現.并行算法實現的流程如圖1所示,粒子信息的初始數據存儲在CPU中,將數據傳輸到GPU中之后,后續的所有計算都在GPU中進行,即所有涉及粒子循環的任務執行都是由GPU的并行結構來完成的.在完成控制方程求解和粒子物理量更新之后,將更新后的粒子物理量從GPU再次傳輸到CPU進行保存和輸出.

圖1 SPH并行算法的實現流程Fig.1 Parallel arithmetic mode of SPH

2.2 粒子搜索

SPH方法由于核函數的緊支性,只有粒子支持域內的相鄰粒子才會相互作用,因粒子間的相對位置不斷變化,每個時刻都需要進行最近相鄰粒子的搜索,所以搜索算法極大地影響著SPH的計算效率.本文采用Monahan[15]的鏈表搜索算法,通過使用控制網格搜索粒子的相鄰質點.其主要思想為:在計算域建立臨時網格,網格大小和粒子支持域大小一致,然后建立每個粒子與所處網格的關聯形成關聯鏈表.對于三維問題,計算域被劃分為立方單元,每個單元與其周圍的單元(共27個)組成緊鄰單元,如圖2所示.粒子根據它們在區域中的位置被放置在單元格中,在搜索粒子時只需搜索粒子所在的緊鄰單元即可,形成最近相鄰粒子對.

圖2 三維相鄰粒子鏈表搜索法示意圖Fig.2 Sketch for three-dimensional neighbor list search method

2.3 粒子重新編號技術

鄰近粒子搜索的GPU并行算法的構建與在CPU上使用的過程有所不同,在GPU上需要建立線程與粒子之間的關聯,并為每個粒子單獨分配一個計算線程.所有粒子按編號排入一個連續的標簽數組,并給固定數目的連續粒子分配一個線程塊.在流場模擬過程中,隨著時間步不斷推進,SPH粒子分布與初始粒子分布差異越來越大,鏈表搜索法中SPH粒子不再是順序分布,如圖3(a)所示,這就導致線程結束并行訪問內存時,極易導致多個線程同時對同一內存地址進行讀寫操作而產生訪問沖突,降低并行算法的求解效率.為保證算法效率穩定性,提出了一種粒子重新編號方法,圖3顯示了用于生成按單元格重新編號法的粒子標簽數組簡單示意圖.

圖3 粒子重新編號算法示例Fig.3 Example of neighbor list reorder

該方法主要包含以下步驟:

(1) 遍歷所有SPH粒子,計算每個粒子所在網格索引,統計該索引所在網格存儲SPH粒子數量,存儲于數組IDCount中.

(2) 創建數組IDBegin,記錄新粒子編號中每個網格中首粒子編號,對于網格m有IDBegin[m] = IDBegin[m-1] + IDCount[m-1];

(3) 清空數組IDCount;

(4) 遍歷所有SPH粒子,重新統計該索引所在網格存儲SPH粒子數量,根據粒子所在網格索引,對其重新編號,如對于某SPH粒子,其對于網格索引為n,其新的粒子編號為IDBegin[n] + IDCount[n],其中IDCount[n]為統計至該粒子時網格n中存儲的SPH粒子數量,最終得到的重新粒子編號效果如圖3(b)所示.

3 結果驗證與分析

以SPH法驗證自由表面流動問題的基準測試案例——潰壩問題為研究對象,對三維帶障礙物的潰壩進行數值模擬,通過與實驗數據進行對比驗證本文并行算法的可靠性和有效性.

3.1 潰壩模型

荷蘭海洋研究所(Marin)的潰壩試驗[16]被廣泛認為是驗證SPH自由表面流動的基準算例.試驗包含一個潰壩流與障礙物的碰撞,如圖4所示.水箱長3.22 m、寬1 m、高1 m,水柱被儲存在水箱一端,長1.228 m、寬1 m、高0.55 m,并在試驗開始瞬間釋放.障礙物設置在水流下游,隨著擋水墻的拆除,由于重力作用,流體逐漸淹沒水箱干床并與障礙物發生碰撞.試驗通過設置3個垂直高度探頭(H1、H2和H3)測量不同位置的水深,具體位置如圖4所示.

圖4 試驗配置和試驗數據的測量位置(m)Fig.4 Experimental configuration and measurement position of experimental data (m)

3.2 結果對比

數值模擬選取初始粒子間距為0.01 m、流場實粒子總數為 676 500,時間步長為0.1 ms.為驗證本文粒子重新編號方法的準確性,對不同位置處的水深展開比較研究.圖5所示為 H1、H2和H3位置處不同時刻的水深粒子重新編號前后計算值和實驗值的比較.圖中:hw為水柱高度.可見粒子重新編號前后計算值基本一致,不會對精度造成影響;在初始水柱中間(H3)、箱體中部(H2)、障礙物前緣附近(H1),本文計算值與實驗值的趨勢和大小都很接近.圖中H3探頭清晰地再現了潰壩的整個過程,在最初的2 s內水柱坍塌,相應地,這段時間的水位也不斷下降,而在其他探頭中,水流依次到達H2和H1.在1.75 s后,反射的水波撞擊左墻后反向移動,第2次撞擊障礙物以及右墻,同時,水深達到第2個峰值(H1的水深峰值時刻大概為4.8 s,H2為 4.6 s,H3為3.8 s).本文計算的水深與實驗水深隨時間的變化趨勢大致相同,表明本文方法具有良好的計算精度.

圖5 實驗測量和本文模擬的探測點處的垂直水柱高度比較Fig.5 Comparation of vertical water heights at the detection point measured experimentally and simulated in this paper

4 效率測試結果與分析

4.1 測試平臺

通過比較CPU串行、CPU并行、GPU并行SPH算法三維帶障礙物潰壩數值模擬計算耗時,驗證本文所采用GPU算法求解效率.為保證計算結果的適用性,分別在不同的CPU、GPU硬件上運行SPH算法,具體配置如表1所示.

表1 CPU和GPU配置Tab.1 Hardware configurations of CPU and GPU

4.2 GPU算法粒子重新編號效果對比

圖6為實粒子總數676 500時不同硬件條件下粒子重新編號和未重新編號GPU加速SPH算法單步運行時間(ts)對比.圖中:sn為計算步數.由圖可見,本文采用的粒子重新編號算法在不同硬件上都獲得了穩定的單步運行時間,而粒子未重新編號時,隨著流場中粒子的無序化導致GPU顯存訪問沖突,其單步運行時間呈對數增長,算例表明本文所采用的粒子重新編號方法可以保證穩定的單步運行時間,是有效的.

圖6 粒子重新編號和未重新編號單步運行時間對比Fig.6 Comparison of step running time between reorder and non-reorder method

4.3 算法求解效率對比

圖7為實粒子數為676 500時CPU并行和GPU并行算法迭代60 000步單步運行時間對比,可見GPU并行算法都有良好的計算效率,而求解效率穩定性弱于CPU并行算法.受SPH方法部分算法、函數間存在串行關系影響,計算效率未能隨核數增加而線性增加,以Intel Core i9-10900F為參照,各硬件核數、效率比如表2所示,可見相對計算成本(核數效率比)隨核數的增加而增大.圖8為不同實粒子數下CPU并行、GPU并行SPH算法單步平均用時(tm)對比.圖中:np為實粒子數.可見隨著粒子數量的增加,CPU并行算法運行時間顯著增加, GPU并行算法大幅縮短計算時間的優勢愈發明顯.

表2 各硬件核數和效率比

圖7 CPU并行和GPU算法單步運行時間對比Fig.7 Comparison of step running time between CPU parallel and GPU parallel

圖8 不同實粒子總數下CPU與GPU運行時間的對比Fig.8 Comparison of running time between CPU and GPU at different number of particles

為進一步驗證本文粒子重新編號算法的有效性,對不同實粒子數下GPU算法并行效率與開源軟件DualSPHysics進行了對比,如圖8(b)所示,可見同實粒子總數、同硬件設備條件下,本文方法平均單步運行時間均小于DualSPHysics軟件,算例表明本文粒子重新編號方法具有良好的效率優勢.

5 結論

運用粒子重新編號技術開發了一套高效GPU-SPH并行算法,將該算法應用于三維帶障礙物潰壩問題,并對算法求解效率進行了比較研究,得到以下結論:

(1) 粒子重新編號前后計算值基本一致,不會對精度造成影響,與試驗值的對比表明本文所采用的方法精確有效.

(2) 粒子重新編號技術能夠有效解決GPU-SPH算法中的顯存訪問沖突問題.

(3) GPU并行算法能夠大幅提高SPH方法求解效率,隨著粒子數量的增加,其大幅縮短計算時間的優勢愈發明顯.