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

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

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

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

1 SPH的控制方程

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

(1)

(2)

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

(3)

(0≤q≤2)

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

人工黏性為

(4)

式中:

(5)

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

(6)

2 CUDA-GPU架構(gòu)下的SPH模型

2.1 并行算法的實現(xiàn)流程

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

圖1 SPH并行算法的實現(xiàn)流程Fig.1 Parallel arithmetic mode of SPH

2.2 粒子搜索

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

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

2.3 粒子重新編號技術(shù)

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

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

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

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

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

(3) 清空數(shù)組IDCount;

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

3 結(jié)果驗證與分析

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

3.1 潰壩模型

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

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

3.2 結(jié)果對比

數(shù)值模擬選取初始粒子間距為0.01 m、流場實粒子總數(shù)為 676 500,時間步長為0.1 ms.為驗證本文粒子重新編號方法的準確性,對不同位置處的水深展開比較研究.圖5所示為 H1、H2和H3位置處不同時刻的水深粒子重新編號前后計算值和實驗值的比較.圖中:hw為水柱高度.可見粒子重新編號前后計算值基本一致,不會對精度造成影響;在初始水柱中間(H3)、箱體中部(H2)、障礙物前緣附近(H1),本文計算值與實驗值的趨勢和大小都很接近.圖中H3探頭清晰地再現(xiàn)了潰壩的整個過程,在最初的2 s內(nèi)水柱坍塌,相應(yīng)地,這段時間的水位也不斷下降,而在其他探頭中,水流依次到達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 效率測試結(jié)果與分析

4.1 測試平臺

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

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

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

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

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

4.3 算法求解效率對比

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

表2 各硬件核數(shù)和效率比

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

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

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

5 結(jié)論

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

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

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

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