基于GPU加速的Boussinesq類波浪傳播變形數(shù)值模型

孫家文，朱 桐，房克照，劉忠波

(1. 大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室 DUT-UWA海洋工程聯(lián)合研究中心，遼寧 大連 116024； 2. 國家海洋環(huán)境監(jiān)測中心，遼寧 大連 116023； 3. 大連海事大學(xué) 交通運(yùn)輸工程學(xué)院，遼寧 大連 116026)

作為一種典型的水波方程，Boussinesq類方程同時(shí)考慮色散性與非線性特征，能較準(zhǔn)確地描述近岸海域波浪的傳播演變過程，再現(xiàn)波浪在傳播過程中出現(xiàn)的折射、繞射、散射、反射、淺水變形等物理現(xiàn)象。因此，Boussinesq方程被廣泛應(yīng)用于近海區(qū)域波浪的數(shù)值計(jì)算[1-2]。Peregrine[3]于1967年推導(dǎo)了適應(yīng)變水深的二維Boussinesq方程，該方程被稱為經(jīng)典Boussinesq方程。但是該方程僅含有弱非線性和弱色散性項(xiàng)，難以適應(yīng)大水深及非線性較強(qiáng)的問題。針對這類問題，Madsen和S?rensen[4]于1992年推導(dǎo)了二維擴(kuò)展型Boussinesq水波方程(后稱MS92方程)。該方程在動(dòng)量方程中加入具有特定系數(shù)的三階偏導(dǎo)數(shù)項(xiàng)，使得方程的色散性提高了一階。因?yàn)镸S92方程形式簡單，又具有較好的色散性，得到廣泛應(yīng)用。

以Boussinesq方程建立的時(shí)域數(shù)值模型屬于相位解析模型，為保證計(jì)算精度則需要較高的時(shí)間和空間分辨率。研究表明，單純使用有限差分方法求解這類方程易產(chǎn)生高頻數(shù)值震蕩，影響計(jì)算的穩(wěn)定性，通常需要采用具有耗散性質(zhì)的光滑器以滿足穩(wěn)定性的需求[5-6]。利用有限體積法求解Boussinesq方程一定程度上可解決上述問題[7-8]。求解時(shí)，采用有限體積方法計(jì)算對流項(xiàng)，剩余項(xiàng)則采用有限差分方法求解，這一方式經(jīng)過驗(yàn)證已經(jīng)取得了較好的效果[9-10]。但是利用有限體積法離散對流項(xiàng)時(shí)，需要對界面左右變量進(jìn)行狀態(tài)重構(gòu)和數(shù)值通量求解，過程較為復(fù)雜，耗時(shí)較高，當(dāng)網(wǎng)格數(shù)目較多時(shí)更會(huì)大幅度增加計(jì)算時(shí)間。若想同時(shí)保證計(jì)算的精確性與高效性，單個(gè)主機(jī)的計(jì)算性能已經(jīng)很難滿足要求。采用主機(jī)與圖形處理器相結(jié)合的方式搭建高性能計(jì)算平臺(tái)是解決辦法之一。迄今，眾多CPU上的算法已成功移植到GPU上，計(jì)算效率得到了明顯提升[11-12]。近兩年，GPU已經(jīng)開始使用到Boussinesq類模型的數(shù)值計(jì)算當(dāng)中。Tavakkol[13]開發(fā)了一種基于GPU加速的Boussinesq類波浪模擬軟件，不僅能提高計(jì)算效率，而且交互性強(qiáng)，提供了包含真實(shí)化渲染等很多可視性選項(xiàng)。Kim等[14]使用CUDA FORTRAN語言實(shí)現(xiàn)了Boussinesq方程模型的GPU并行化，計(jì)算效率相比CPU模型明顯提升，最大加速比(CPU時(shí)間/GPU時(shí)間)可達(dá)20以上。GPU具有高并行性與強(qiáng)大的數(shù)據(jù)處理能力，能很好地解決當(dāng)下數(shù)值模擬中網(wǎng)格數(shù)量多和計(jì)算效率低的問題。將來，采用GPU的海洋數(shù)值模型是海洋研究發(fā)展的重要方向之一。

利用CUDA C語言搭建基于GPU加速的Boussinesq類波浪傳播變形數(shù)值模型，并模擬一系列經(jīng)典工況，把所得的結(jié)果與解析解及CPU模型的模擬結(jié)果相比較，以驗(yàn)證模型的計(jì)算精度。除此之外，還比較了同種工況采用GPU模擬與CPU模擬的運(yùn)算效率，得到了計(jì)算效率與網(wǎng)格數(shù)目之間的關(guān)系。

1 數(shù)學(xué)模型

1.1 控制方程

MS92方程的守恒形式如下[15]：

Ut+F(U)x+G(U)y=S(U)

(1)

式中：U為變量矢量，F(xiàn)和G為x與y方向上的通量矢量，定義為：

(2)

其中，

(3)

式中：d=η+h為總水深，η表示自由波面，h為靜水水深，P和Q分別表示x和y方向上的速度通量(P=du和Q=dv，u、v分別表示x、y方向的水深平均速度)，B為色散性參數(shù)，取1/15。為實(shí)現(xiàn)海床變化和海岸水-陸動(dòng)邊界存在時(shí)高精度有限體積方法的和諧解，引入水位ζ=zb+d。式(1)中S為源項(xiàng)，可將其分成三個(gè)組成部分，分別為水底坡度項(xiàng)Sb、水底摩擦項(xiàng)Sf和色散項(xiàng)Sd，如下：

(4)

式中：色散項(xiàng)為：

(5)

水底摩擦項(xiàng)：

(6)

式中：f為底摩擦系數(shù)，取值范圍為0.000 1～0.01，本文所有工況計(jì)算均未考慮海床摩擦。

1.2 方程空間離散、時(shí)間積分與模型邊界條件

采用矩形網(wǎng)格單元，將計(jì)算域在空間、時(shí)間上均勻劃分單元網(wǎng)格，波面η和速度u均定義在控制體中心。為了提高空間精度，采用四階高精度狀態(tài)差值方法(MUSCL)對界面左右變量進(jìn)行重構(gòu)，應(yīng)用重構(gòu)后的界面變量進(jìn)行通量計(jì)算。模型的時(shí)間積分采用具有TVD性質(zhì)的三階龍格-庫塔方法[16]。為保證計(jì)算收斂，時(shí)間計(jì)算步長由Courant-Friedrichs-Lewy(CFL)穩(wěn)定性限制條件確定：

(7)

式中：ν為常數(shù)，這里取ν= 0.25。

對于文中計(jì)算工況，波浪通過給定初始條件或內(nèi)部造波方法產(chǎn)生，計(jì)算域末端均為固壁邊界，視情況不同在固壁前設(shè)置海綿層以吸收波浪能，以有效降低數(shù)值計(jì)算中的二次反射[15]。

2 GPU并行計(jì)算簡介

在CUDA的并行計(jì)算模型中，CPU作為主機(jī)端，負(fù)責(zé)運(yùn)行主程序并在運(yùn)行過程中調(diào)用核函數(shù)。GPU作為設(shè)備端，主要進(jìn)行線性化的并行計(jì)算。在設(shè)備端運(yùn)行的函數(shù)被稱為核函數(shù)。每個(gè)核函數(shù)對應(yīng)的線程均以線程網(wǎng)格(grid)的形式存在，線程網(wǎng)格又被分成多個(gè)線程塊(block)，每個(gè)線程塊又由大量的線程(thread)組成。文中模型采用每個(gè)線程塊1 024個(gè)線程的分配方案。

由于GPU無法直接讀寫內(nèi)存中的數(shù)據(jù)，因此在核函數(shù)運(yùn)行先需要初始化數(shù)據(jù)并將其由主機(jī)端傳遞到設(shè)備端的全局存儲(chǔ)器中；之后，各個(gè)線程塊會(huì)分配到不同的GPU大核 (Streaming Multiprocessors，簡稱SM)上，線程塊上的不同線程交給GPU大核上不同的CUDA核心 (Streaming Processor，簡稱SP)來執(zhí)行，從而達(dá)到并行運(yùn)算的效果，如圖1所示。為充分發(fā)揮GPU計(jì)算效率，在程序設(shè)計(jì)時(shí)，將主循環(huán)中計(jì)算耗時(shí)較重的模塊，如狀態(tài)重構(gòu)、通量計(jì)算、時(shí)間積分等設(shè)計(jì)為核函數(shù)，在設(shè)備端運(yùn)行。模型的初始化參數(shù)以及輸入輸出等模塊在主機(jī)端部署。設(shè)備端和主機(jī)端之間通過CUDA C提供的內(nèi)置函數(shù)(memcpy)進(jìn)行數(shù)據(jù)傳遞，模型計(jì)算流程見圖2。

圖1 CUDA編程模型Fig. 1 CUDA programming model

3 數(shù)值驗(yàn)證

針對所建立的GPU模型，通過幾個(gè)典型工況的數(shù)值模擬分別對其進(jìn)行驗(yàn)證。其中工況一為常水深孤立波的長距離傳播問題；工況二為規(guī)則波跨越潛堤傳播；工況三為封閉方形水池中的水面晃動(dòng)；工況四為波浪在橢圓形淺灘上的傳播。

3.1 孤立波在常水深水槽中的傳播

在本工況中，x方向計(jì)算域長度為409.6 m，靜水水深h0為1.0 m，在計(jì)算域內(nèi)給出孤立波解析解作為初始條件，波高H=0.6 m，孤立波從計(jì)算域的左側(cè)向右側(cè)傳播，初始時(shí)刻波峰位于x0=25 m處，網(wǎng)格尺寸Δx=0.2 m，忽略底摩擦的影響，且計(jì)算域邊界處均采用固壁邊界條件。圖3給出了模擬的孤立波在t=0，10，20，30，40，50，60，70和80 s時(shí)的波面狀態(tài)。從圖中可見，經(jīng)過長距離的傳播，孤立波的波形和波高值均保持穩(wěn)定，幾乎未發(fā)生變化，說明數(shù)值離散正確，未引入額外的數(shù)值耗散和數(shù)值彌散。圖4為時(shí)間t=80 s時(shí)，模型計(jì)算結(jié)果與解析解的比較，可以看出兩者高度吻合，說明模型具有較好的精確性。

圖4 t= 80 s時(shí)數(shù)值模擬結(jié)果與解析解的對比Fig. 4 Comparison of simulated and analytical surface profiles at t=80 s

3.2 規(guī)則波在潛堤海床上的傳播

規(guī)則波跨越水下潛堤傳播是非常復(fù)雜的波浪演化過程，涉及波浪的變淺、高次諧波產(chǎn)生釋放等，該類試驗(yàn)被用來檢驗(yàn)各類數(shù)值模型的色散性和非線性等綜合性能。采用Dingemans[17]的規(guī)則波跨越潛堤傳播的試驗(yàn)來對數(shù)值模型進(jìn)行檢驗(yàn)，試驗(yàn)海床設(shè)置如圖5中所示，海床兩端設(shè)置1.5倍波長海綿層吸收波浪，造波機(jī)在x= 0 m處，平底處靜水水深為0.4 m，梯形壩的坡腳位于x=6.0 m處，梯形壩頂端長為2.0 m，距離靜水面0.1 m，梯形潛堤前坡坡度為1∶20，后坡坡度為1∶10，根據(jù)Dingemans[17]的試驗(yàn)，在計(jì)算域內(nèi)設(shè)置10個(gè)浪高儀，位置分別為x= 2，4，10.5，12.5，13.5，14.5，15.7，17.3，19和21 m。入射波波高H=0.02 m，周期T=2.02 s，計(jì)算網(wǎng)格長度為Δx=0.025 m。

圖5 規(guī)則波跨越潛堤傳播試驗(yàn)布置Fig. 5 Experimental setup of regular waves propagation over submerged bar

圖6給出了10個(gè)浪高儀位置上數(shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)和CPU模型模擬結(jié)果的對比。如圖6所示，規(guī)則波跨越堤頂傳播過程中，水深變淺，波浪的非線性增強(qiáng)，高次諧波產(chǎn)生，波形發(fā)生明顯變化，潛堤前以及堤頂位置處波面的數(shù)值計(jì)算結(jié)果與Dingemans的試驗(yàn)數(shù)據(jù)吻合程度高，說明GPU模型能夠準(zhǔn)確模擬規(guī)則波的傳播及跨越潛堤時(shí)的波浪狀態(tài)。堤后伴隨水深的增加，各次諧波變成自由波，計(jì)算結(jié)果與實(shí)測值存在一定誤差，這主要是方程捕捉高次諧波的色散性和非線性均不夠造成的。整體來看，GPU模型與CPU模型的計(jì)算結(jié)果高度一致，驗(yàn)證了GPU模型的正確性。

圖6 規(guī)則波跨越潛堤傳播各浪高儀處的波面時(shí)間序列Fig. 6 Time series of surface elevations for regular waves passing over a submerged bar

3.3 封閉方形水池中的水面晃動(dòng)

為了驗(yàn)證模型精度以及對全反射邊界的處理能力，對封閉方形水池中液體晃動(dòng)進(jìn)行了模擬。數(shù)值水池為正方形，尺寸為7.5 m×7.5 m，水深h= 0.45 m。水池四邊為全反射直墻且與x、y軸平行，坐標(biāo)原點(diǎn)位于水池左下角。網(wǎng)格尺寸Δx=Δy=0.058 593 75 m。初始波面為Gaussian形，即：

η0(x,y)=H0exp{-2[(x-a)2+(y-b)2]}

(8)

式中：a=b=3.75 m，H0=0.045 m。

圖7給出了波面在一段時(shí)間后的分布情況。如圖7所示，在t= 2 s時(shí)，波面關(guān)于x軸和y軸完全對稱。經(jīng)過比較長時(shí)間的傳播后，在t=20 s時(shí)，對稱性仍能很好地保持，說明模型離散和數(shù)值計(jì)算正確。

圖7 t=2 s與t=20 s時(shí)波面分布Fig. 7 Free-surface contours at time t=2 s and t=20 s

圖8給出了H0= 0.045 m時(shí)水池角點(diǎn)處和中心點(diǎn)處自由表面計(jì)算結(jié)果和線性理論結(jié)果的比較。波浪傳播過程中前20 s內(nèi)的非線性相對較弱，因而基于GPU加速的弱非線性Boussinesq數(shù)值模型與解析解結(jié)果基本一致。在后20 s內(nèi)，伴隨邊界反射疊加，波浪的非線性增強(qiáng)，組成波成份增多，超出了方程本身非線性和色散性使用范圍，導(dǎo)致二者存在一定的誤差。

圖8 H0=0.045 m時(shí)波面時(shí)間過程與線性理論計(jì)算結(jié)果比較Fig. 8 Comparison of time series of free surface elevation between numerical and theoretical results for H0=0.045 m

3.4 波浪在橢圓形淺灘上的傳播

為了檢驗(yàn)?zāi)Ｐ湍M復(fù)雜地形上波浪傳播變形的能力，對Berkhoff等[18]所做的橢圓形淺灘上的波浪折、繞射物理模型試驗(yàn)進(jìn)行了數(shù)值模擬。在模擬中，計(jì)算區(qū)域長30 m，寬20 m，具體試驗(yàn)地形如圖9所示。

模擬波浪入射波高為H0=0.046 4 m，周期T=1 s，在x=3 m處沿x軸正向入射。上下兩邊界(y=0 m和y=20 m)為全反射邊界，左側(cè)海綿層寬度為2 m，右側(cè)海綿層寬度為3 m。網(wǎng)格數(shù)目為512(x方向)×256(y方向)=131 072。模型運(yùn)行時(shí)間為40 s。

將GPU模型計(jì)算結(jié)果與Tonelli和Petti[8]計(jì)算結(jié)果進(jìn)行了對比，見圖10?？梢姡灸Ｐ偷臄?shù)值計(jì)算結(jié)果與Tonelli和Petti[8]所得結(jié)果基本一致，說明本模型可以較好地模擬波浪的折射和繞射等現(xiàn)象。

圖9 試驗(yàn)地形圖(Berkhoff等[18])Fig. 9 Topographic map of experiment (Berkhoff, et al [18])

圖10 GPU模型計(jì)算波高與試驗(yàn)數(shù)據(jù)和他人結(jié)果比較Fig. 10 Comparisons of normalized wave height from GPU model, data and other simulation

4 GPU并行模型加速性能分析

本文計(jì)算使用專業(yè)圖形加速卡Tesla k20進(jìn)行，為了對比也使用CPU模型進(jìn)行了計(jì)算。

具體設(shè)備參數(shù)和開發(fā)平臺(tái)如下：

CPU：Intel i5 8500，主頻3.0 GHz，含6個(gè)CPU核心。

GPU：專業(yè)圖形加速卡Tesla k20，含2 688個(gè)流處理器，主頻732 MHz，單精度浮點(diǎn)性能3.95 Teraflops，顯存5.2 GHz 6GB，帶寬250 GB/s。

開發(fā)平臺(tái)：Intel Complier Visual studio 2015，CUDA 8.0。

利用同一臺(tái)設(shè)備分別使用GPU模型與CPU模型進(jìn)行計(jì)算，統(tǒng)計(jì)計(jì)算時(shí)間，并計(jì)算出加速比(CPU模型計(jì)算時(shí)間TC與GPU模型計(jì)算時(shí)間TG之比)，具體結(jié)果見表1。從表1可知，隨著網(wǎng)格數(shù)增加，加速比穩(wěn)定提高。同時(shí)，不同的工況中的加速比一般不同，其中封閉方形水池的水面晃動(dòng)所達(dá)到的加速比最大，計(jì)算效率的提升最為明顯，最大加速比可達(dá)12.72。綜合來看，GPU并行模型的確對計(jì)算效率有很明顯的提升。

Kim等[14]發(fā)展的GPU類Boussinesq模型加速比可達(dá)20，高于本文建立模型，可能是如下原因所致：1)圖形處理器不同，本文選用的圖形處理器是NVIDIA GeForce Tesla k20，Kim等選用的是NVIDIA GeForce GTX TITAN Z；2)編程框架不同，本文的編程語言是CUDA C，而Kim等選擇的是CUDA FORTRAN；3)計(jì)算工況不同，工況的差異會(huì)影響加速比，這在本文中已經(jīng)得到了證明；4)算法不同，不同的算法由于數(shù)據(jù)讀取方式的不同，會(huì)對計(jì)算效率有著較大影響。

表1 各工況CPU模型與GPU模型計(jì)算效率比較Tab. 1 Computational efficiency comparison between CPU model and GPU model

5 結(jié) 語

利用CUDA C語言與圖形處理器(GPU)實(shí)現(xiàn)了Boussinesq模型的并行化。采用有限差分與有限體積的混合格式來求解方程。在計(jì)算過程中，通過把較為復(fù)雜的計(jì)算任務(wù)分配到GPU上分塊并行處理，從而實(shí)現(xiàn)大大提升計(jì)算效率的目的。將本模型的計(jì)算結(jié)果與CPU模型上的計(jì)算結(jié)果和模型本身的解析解相比較，所得結(jié)果基本一致，驗(yàn)證了本模型計(jì)算的正確性。同時(shí)對本模型與CPU模型的計(jì)算效率進(jìn)行了對比。通過比較發(fā)現(xiàn)，本模型對計(jì)算效率有很明顯的提升，文中工況的最大加速比可達(dá)12以上。

進(jìn)一步的研究將圍繞不同圖形處理器設(shè)備加速性能、并行優(yōu)化以及實(shí)際應(yīng)用方面展開。