楊志國，黃興，鄭興，段文洋

（哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001）

GPU在SPH方法模擬潰壩問題的應(yīng)用研究

楊志國，黃興，鄭興，段文洋

（哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001）

SPH方法是一種無網(wǎng)格的粒子方法，對于求解強非線性水動力學(xué)問題具有重要意義。隨著粒子數(shù)增加，該方法的計算效率成為限制其大規(guī)模工程應(yīng)用的重大瓶頸?？蓪⒋笠?guī)模并行計算引入SPH方法中，以得到良好的計算加速效果。采用將GPU運用于SPH方法并行計算的技術(shù)，借助CUDA硬件計算架構(gòu)，研究SPH方法的并行計算通用性問題。以二維潰壩問題作為數(shù)值算例，對GPU計算結(jié)果的穩(wěn)定性和收斂性進(jìn)行驗證，比較CPU與GPU的計算效率。通過計算，驗證了GPU在SPH方法并行計算應(yīng)用中的可靠性、可行性以及高效性，為提高SPH方法的計算效率提供一種重要的參考途徑。

GPU；并行計算；CUDA；SPH方法；潰壩；水動力學(xué)；數(shù)值計算；

SPH方法是一種基于拉格朗日力學(xué)的無網(wǎng)格粒子法，該方法通過將連續(xù)的流場離散成一系列的流體粒子，而后用離散的流體粒子通過核近似的方式代替連續(xù)流場，并求解離散后的流場動力學(xué)方程。由于SPH方法的無網(wǎng)格特性以及帶有相應(yīng)物理量的粒子的自由運動的特點，使得SPH方法在研究帶有強非線性自由表面的流場問題等領(lǐng)域中有著獨特的優(yōu)勢。SPH方法最早是在求解天體物理問題［1-2］時提出，后經(jīng)一系列研究發(fā)展，其研究領(lǐng)域得到不斷拓展，如應(yīng)用于可壓縮非耗散流體問題［3］、彈性和斷裂問題［4］、弱可壓縮流體自由表面流動問題［5］、波浪破碎問題［6-7］。

隨著計算機科學(xué)的不斷發(fā)展以及科學(xué)計算中對計算機的性能要求不斷提高，圖形處理器（graphics processing unit，GPU）作為新型的計算處理單元，以其強大的并行計算能力，滿足了諸多領(lǐng)域的計算需求，近年來得到迅猛的發(fā)展。其功能已由最初的圖形處理加速的簡單需求，轉(zhuǎn)變?yōu)樾滦偷拿嫦蛴嬎銠C通用計算平臺的通用計算圖形處理器（general purpose GPU，GPGPU）。目前，GPU在物理［8］、生物［9］、化學(xué)［10］等諸多科學(xué)領(lǐng)域已經(jīng)得到較為廣泛的應(yīng)用。

針對SPH方法計算過程中的效率問題，將GPU應(yīng)用于SPH方法計算中，是提高計算實效的有效途徑。目前，SPH的GPU并行研究主要集中于流體實時模擬問題［11］、近海海浪仿真［12］等問題。

由于GPU平臺下SPH方法并行計算結(jié)果可靠性問題研究目前尚不充分，本文首先介紹了基于GPU平臺的SPH并行計算方法，之后著重對GPU模擬二維潰壩問題數(shù)值結(jié)果的收斂性與穩(wěn)定性進(jìn)行詳細(xì)研究，探討GPU并行架構(gòu)在SPH方法求解潰壩問題中的可行性。之后，根據(jù)GPU硬件特性和數(shù)值求解結(jié)果，對CPU-GPU計算效率以及GPU內(nèi)線程數(shù)量對計算效率的影響進(jìn)行重點分析，從而為GPU上實現(xiàn)SPH方法并行計算，驗證求解結(jié)果可靠性，提高SPH方法計算效率等問題提供一定參考。

1 SPH方法

式中：Ω為包含x的積分體積，f為三維坐標(biāo)向量x的函數(shù)；δ( x-x′，h )為δ函數(shù)，性質(zhì)如下：

在SPH方法中f x()函數(shù)的積分表達(dá)式定義為

經(jīng)離散化后N-S方程質(zhì)量守恒，動量守恒的SPH表達(dá)式可寫為

本文選取Monaghan［13］提出的三次樣條函數(shù)作為核函數(shù)。邊界粒子布置采用鏡像法［14］，鏡像粒子滿足方向不可穿透，鏡像可滑移條件。此外，為確保數(shù)值結(jié)果的穩(wěn)定性，傳統(tǒng)的弱可壓縮SPH方法需考慮人工壓縮性、人工粘性［5］等問題，對于人工粘性問題，特別是對粘性系數(shù)選取問題，許多學(xué)者進(jìn)行了深入研究［15］。

2 CUDA架構(gòu)

GPU一般存在于計算機的顯卡中，擁有多處理器結(jié)構(gòu)。隨著計算機科學(xué)的不斷發(fā)展，GPU已經(jīng)從原本的圖形處理器逐步轉(zhuǎn)變?yōu)檫m合處理具有較高計算密度和簡單邏輯分支任務(wù)的通用處理器。

目前，市場上出售的GPU主要是NVIDIA和AMD／ATI兩家公司的產(chǎn)品。本文適用NVIDIA公司GTX660產(chǎn)品。

2.1 CUDA及存儲模型

統(tǒng)一計算設(shè)備架構(gòu)（compute unified device architecture，CUDA）是NVIDIA公司于2007年6月推出的計算架構(gòu)體系，它為GPU應(yīng)用于數(shù)據(jù)并行計算提供了一套完備的硬件架構(gòu)和軟件體系。CUDA編程模型中，分為主機（Host）與設(shè)備（Device）或者主機（Host）與協(xié)處理器（co-processor）［16］。主機是CPU在CUDA中的代號，而設(shè)備、協(xié)處理器則是GPU的代號，之所以將二者分開來，主要是為了便于在調(diào)用CPU或GPU時區(qū)分兩者功能。GPU上每個完整的處理核心稱為一個core。

設(shè)備端內(nèi)存種類較多，以其在執(zhí)行過程中作用不同，可分為全局存儲器、常數(shù)存儲器、紋理存儲器、局部存儲器、共享寄存器和寄存器［16］。

2.2 執(zhí)行方式及執(zhí)行單位

CUDA執(zhí)行以主機端控制設(shè)備端，設(shè)備端返回給主機端的執(zhí)行結(jié)果的方式進(jìn)行。當(dāng)執(zhí)行任務(wù)時，先由主機端程序開辟任務(wù)所需的內(nèi)存或顯存，然后將主機端數(shù)據(jù)拷入顯存中，設(shè)備端接收到主機端的執(zhí)行指令后執(zhí)行相應(yīng)的任務(wù)要求，執(zhí)行結(jié)束后，將執(zhí)行結(jié)果傳遞給主機端，主機端釋放開辟的內(nèi)存或顯存，執(zhí)行結(jié)束，執(zhí)行過程如圖1所示。

圖1 GPU-CPU執(zhí)行模式Fig.1 GPU-CPU execution model

CUDA線程是設(shè)備端運行的基本單位。線程是以Grid、Block、Thread為層次機構(gòu)執(zhí)行的。每個Block由若干個Thread組成一維、二維或者三維形式，而每個Grid則由若干個Block組成一維、二維或者三維形式。在設(shè)備端執(zhí)行的過程中，線程作為最小執(zhí)行形式被分派到Block中，而Block又被分派到Grid中，控制器將Block發(fā)送到流多處理器中，進(jìn)而將線程發(fā)送到流處理器中計算相應(yīng)指令，執(zhí)行結(jié)束后完成設(shè)備端計算任務(wù)，返回主機端。

2.3 SPH并行模式

本文在SPH方法的每步迭代過程中，以單線程控制單個流場粒子，以多線程同時執(zhí)行的方式計算不同粒子的相應(yīng)流場物理量。線程數(shù)量的多少取決于GPU硬件的配置，而并行效率則與線程數(shù)量和GPU計算能力有很大的關(guān)系。

初始時刻，先將粒子信息考入顯存當(dāng)中，在GPU上完成所有計算任務(wù)，之后將所得結(jié)果返回CPU，減少了由于CPU-GPU數(shù)據(jù)通訊花費的計算時間。每步計算過程中，采用鏈表查找法對所有粒子進(jìn)行排序，采用改進(jìn)歐拉迭代法，具體流程如圖2所示。

圖2 SPH并行方式Fig.2 SPH parallel manner

3 二維潰壩模擬

本節(jié)以二維潰壩問題為研究對象，驗證GPU并行方法的可靠性與實用性問題。計算平臺為Intel i5-3570 CPU，8 G內(nèi)存，GTX660顯卡。

3.1 潰壩模型

潰壩模型如圖3所示。

圖3 潰壩模型Fig.3 Dam breaking model

水柱寬度B＝0.2 m，高度H＝0.4 m。矩形體寬度為1 m，光滑長度h＝0.01 m，粒子總數(shù)為800，初始時刻粒子間距與光滑長度相等。ρ0＝1 000 kg／ m3，采用弱可壓縮性假設(shè)，時間步長取0.1 ms。

3.2 計算結(jié)果及分析

圖4給出了T為0.1、0.4、0.7、0.9、1.05 s的SPH方法的GPU計算結(jié)果。計算結(jié)果較好地描述了從潰壩開始，潰壩前端觸壁、爬升、翻卷，二次入水等強非線性特征的典型模擬過程。

圖4 不同時刻計算結(jié)果（800粒子）Fig.4 Calculation result at different time（800 particles）

比較水柱前緣位置是驗證計算結(jié)果準(zhǔn)確性的重要參考。取時間步長為0.1 ms，將計算結(jié)果與Martin［17］等的實驗結(jié)果、Lee［18］等的VOF方法、徐卜男［19］SPH方法的數(shù)值結(jié)果進(jìn)行比較，如圖5所示。本文計算值與Lee、徐等數(shù)值模擬結(jié)果吻合較好，但與實驗值相比，數(shù)值結(jié)果普遍呈現(xiàn)出偏快的特點。

圖5 水柱前緣位置比較Fig.5 Comparison of water column front edge position

3.3 收斂性分析

收斂性分析作為驗證數(shù)值結(jié)構(gòu)的有效性與可靠性的重要手段，一直都作為研究論證的重點。本文采取相同粒子數(shù)量下不同時間步長以及相同時間步長下不同的粒子數(shù)量2種方式驗證收斂性。

3.3.1 時間步長影響

考察相同粒子數(shù)下不同時間步長的影響，本文以粒子數(shù)為800時為例，時間步長分別取為0.1、0.05、0.01、0.005 ms。

圖6 不同時間步長比較Fig.6 Comparison of different time steps

圖6 表明，水柱前緣隨著計算時間步長的改變，變化并不顯著，結(jié)果表明數(shù)值計算結(jié)果趨于穩(wěn)定，GPU并行計算取得良好的收斂性。

3.3.2 粒子數(shù)目對潰壩問題的影響

研究流場中粒子數(shù)量多少對數(shù)值結(jié)果的影響，是收斂性分析的重要工作之一。本小結(jié)保持計算域的大小不變，時間步長取為0.01 ms，光滑長度為0.01、0.005、0.004、0.002和0.001 m所對應(yīng)的粒子數(shù)分別為N＝20×40、N＝40×80、N＝50×100、N＝100×200、N＝200×400。

圖7分別取T＝0.3 s、T＝0.6 s、T＝1.0 s、T＝1.1 s 4個時刻自由液面數(shù)值結(jié)果。結(jié)果表明，在潰壩開始的最初階段，粒子數(shù)的數(shù)值結(jié)果自由表面形狀吻合較好，至前緣爬升階段仍較為一致。出現(xiàn)翻卷時，自由表面開始顯現(xiàn)出差異，尤其是水舌前部重新入水以及出水后的形狀有些許差異。但這種差異隨著粒子數(shù)增加而減小，計算結(jié)果趨于穩(wěn)定。

圖7 自由液面不同粒子數(shù)比較Fig.7 Free surface comparison of different particle number

圖8 不同粒子數(shù)水柱前緣對比Fig.8 Front edge contrast of different particle number

從圖8中可以看出，當(dāng)粒子數(shù)增加時，水柱前緣位置呈先慢后快的趨勢，但差距不超過1%。

3.4 計算時間

文章通過將GPU計算耗時長度與CPU計算耗時進(jìn)行對比，研究GPU的加速效率問題。GTX660顯卡擁有960 core，2G顯存，GPU主頻頻率為1 033 MHz，每個線程塊中理論上允許同時使用的最大線程數(shù)量為1 024。CPU采用i5 3570型號，主頻為3.4 GHz。本節(jié)取粒子數(shù)分別為800、3 200、5 000、20 000、80 000，時間步長取0.01 ms，計算150 000步。分析GPU線程對計算速度的影響以及比較CPU-GPU間計算效率問題。

3.4.1 GPU線程影響

本小節(jié)分別取單個Block內(nèi)線程數(shù)量為128、256、512、768，比較分析不同線程數(shù)量對計算效率的影響。圖9顯示了粒子數(shù)相同情況下，取不同線程所需的平均單步計算時間，從圖中可以看出，GPU的計算效率并不是隨著線程數(shù)量的增加而呈現(xiàn)出簡單的遞增關(guān)系，線程數(shù)量越多，計算時間未必越少，在本算例中單塊內(nèi)的線程數(shù)取256時，不同粒子數(shù)下計算時間均為最短，而當(dāng)線程數(shù)增加或減小時，計算所需時間反而增加，且粒子數(shù)不同，增加的趨勢也不盡相同。

圖9 線程時間比較Fig.9 Time comparison of thread spend

3.4.2 CPU-GPU比較

分析比較CPU與GPU計算時間的差異，是本文工作的重點之一。圖10顯示了2種不同計算核心的計算時間對比。

當(dāng)粒子數(shù)量為800時，CPU運算速度較GPU有14%的增速，其原因在于GPU的單線程計算能力不及CPU，在計算量相對較小的情況下，CPU在計算時間上顯現(xiàn)出一定的優(yōu)勢。但隨著粒子數(shù)的逐漸增加，GPU多線程優(yōu)勢開始發(fā)揮，從2.35倍的加速效果遞增至3.89倍，同時，隨著粒子數(shù)的增加，GPU并行計算加速效果更加明顯。

圖10 CPU-GPU時間比較Fig.10 Time comparison of CPU and GPU

4 結(jié)論

本文以GPU為計算平臺，研究SPH方法潰壩問題并行計算的數(shù)值模擬方法，結(jié)果表明：

1）GPU用于SPH方法并行計算研究的可靠性。本文以SPH方法理論和GPU計算結(jié)果為依據(jù)，將數(shù)值結(jié)果與實驗值比較，通過取不同粒子數(shù)以及不同時間步長進(jìn)行收斂性驗證，表明GPU并行計算求解水動力學(xué)問題的可靠性與適用性。

2）GPU并行加速SPH方法的有效性。比較GPU與CPU計算時間，單步GPU的計算加速最高時達(dá)到CPU計算時間的3.89倍，隨著計算量的增加，加速效果逐漸凸顯，由此可見，GPU在研究SPH方法并行計算中存在很大的適用性。

3）在SPH方法計算過程中，GPU線程數(shù)量對計算時間影響較大。通過相同粒子數(shù)下不同線程數(shù)量對計算時間比較可知，線程數(shù)量對計算速度而言有顯著影響，但并不成正比關(guān)系。使用GPU時，綜合考慮硬件特點以及研究問題的特殊性，結(jié)合合理的算法，是取得良好加速效果的重要保證。

［1］GINGOLD R A，MONAGHAN J J.Smoothed particle hydrodynamics：theory and application to non-spherical stars mon［J］.Monthly Notices of the Royal Astronomical Society，1977，181：375-389.

［2］LUCY L B.A numerical approach to the testing of the fission hypothesis［J］.Astronomical Journal，1977，82：1013-1024.

［3］LIBERSKY L，PETSCHEK A G.Smooth particle hydrodynamics with strength of materials advances in free Lagrange methods［J］.Lecture Notes in Physics，1991，395：248-257.

［4］GINGOLD R A，MONAGHAN J J.Kernel estimates as a basis for general particle methods in hydrodynamics［J］.Jour-nal of Computational Physics，1982，46：429-453.

［5］MONAGHAN J J.Simulating free surface flows with SPH［J］.Computational Physics，1994，110：399-406.

［6］MONAGHAN J J.Energy distribution in a particle alpha model［J］.Journal of Turbulence，2004，5：22.

［7］鄭興，馬慶位，段文洋.K2_SPH方法及二維破碎波的模擬［J］.計算物理，2012，29（3）：317-325.ZHENG Xing，MA Qingwei，DUAN Wenyang.K2_SPH method and simulation of 2D breaking waves［J］.Journal of Computational Physics，2012，29（3）：317-325.

［8］JIA X，PAGANETTI H S J，JIANG S B.GPU-based fast Monte Carlo does calculation for proton therapy［J］.Physics in Medicine and Biology，2012，57（23）：7783-7797.

［9］JONES E A，Van ZEIJL R J M，ANDREN P E，et al.High speed data processing for imaging MS-based molecular histology using graphical processing units［J］.Journal of the A-merican Society for Mass Spectrometry，2012，23（4）：745-752.

［10］CHENG Ling，BENKRID K.Design and implementation of a CUDA-Compatible GPU-based core for gapped BLAST algorithm［J］.Procedia Computer Science，2010，1（1）：495-504.

［11］郭秋雷，唐逸之，劉詩秋，等.一個SPH流體實時模擬的全GPU實現(xiàn)框架［J］.計算機應(yīng)用與軟件，2011，28（11）：69-72.GUO Qiulei，TANG Yizhi，LIU Shiqiu，et al.A full GPU implementation framework of SPH fluid real-time simulation［J］.Computer Applications and Software，2011，28（11）：69-72.

［12］陳俊.近海海浪的仿真研究［D］.武漢：武漢理工大學(xué)，2011：13-14.CHEN Jun.Research on simulation of waves in shallowwater［D］.Wuhan：Wuhan University of Technology，2011：13-14.

［13］FLUCK D A，QUINN D W.An analysis of 1-D smoothed particle hydrodynamics Kernels［J］.Journal of Computational Physics，1996，126：699-709.

［14］LIU G R，LIU M B.Smoothed particle hydrodynamics：a meshfree particle method［M］.Singapore：World Scientific Publishing Co.Pte.Ltd，2003：138-141.

［15］ZHENG Xing，DUAN Wenyang.Numerical simulation of dam breaking using smoothed particle hydrodynamics and viscosity behavior［J］.Journal of Marine Science，2010，9：34-41.

［16］張舒，褚艷利，趙開勇.GPU高性能運算之CUDA［M］.北京：中國水利水電出版社，2009：14-16.

［17］MARTIN J C，MOYCE W J.An experimental study of the collapse of liquid columns on a rigid horizontal plane［J］.Philosophical Transactions of the Royal Society of London A，1952，244（882）：312-324.

［18］LEE B H，PARK J C，KIM M H，et al.Numerical simulation of impact loads using a particle method［J］.Ocean Engineering，2010，37（2／3）：164-173.

［19］徐卜男.SPH方法模擬液艙晃蕩及量綱分析［D］.大連：大連理工大學(xué)，2011：45-51.XU Pu′nan.SPH simulation and parametric analysis on sloshing tank［D］.Dalian：Dalian University of Technology，2011：45-51.

The application research of GPU in the SPH method to simulate the dam breaking problem

YANG Zhiguo，HUANG Xing，ZHENG Xing，DUAN Wenyang
（School of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China）

As a mesh-free particle method，the SPH method is very important for dealing with the hydrodynamic problems with strong nonlinear problems.However，following the increase of the particle number，the calculation efficiency becomes a bottleneck for applying the method to the engineering practice.Good acceleration performance can be attained by applying massively parallel computing into the SPH method.In order to study the general parallel purpose of the SPH method，the parallel technology through the use of GPU on the compute unified device architecture platform has been used.The convergence and the stability of the computing results generated by GPU have been verified and the computational efficiency of GPU and CPU has been compared according to the computing results by using a numerical example of the two-dimensional dam-breaking problem.The capability，feasibility and high efficiency of GPU computing have been proved by the results which provide important reference channels to increase the calculation efficiency of the SPH method.

graphics processing unit（GPU）；parallel computing；compute unified device architecture（CUDA）；SPH method；dam breaking；hydrodynamic；numerical simulation

10.3969／j.issn.1006-7043.201305062

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1006-7043.201305062.html

O352

A

1006-7043（2014）06-0661-06

2013-05-23.網(wǎng)絡(luò)出版時間：2014-05-14 15：49：41.

國家自然科學(xué)基金資助項目（51009034，51279041）.

楊志國（1964-），男，副研究員.

楊志國，E-mail：yangzhiguo＠hrbeu.edu.cn.