基于GPU并行算法的水動力數(shù)學(xué)模型建立及其效率分析

趙旭東，梁書秀＊，孫昭晨，劉忠波，2，韓松林，任喜峰

（1.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.大連海事大學(xué) 交通管理學(xué)院，遼寧 大連 116026）

0 引 言

海洋水動力模型的發(fā)展，一方面要求實(shí)現(xiàn)不同動力因素的耦合，如目前風(fēng)－浪－潮耦合模型的需求，這將導(dǎo)致計(jì)算區(qū)域加大；另一方面需要對所關(guān)心區(qū)域不同尺度物理現(xiàn)象的精確模擬，這要求局部網(wǎng)格較小.這兩方面的要求使采用非結(jié)構(gòu)化網(wǎng)格模擬計(jì)算大范圍海域的水動力模型成為一種趨勢.但是網(wǎng)格密度和計(jì)算區(qū)域的增大將導(dǎo)致網(wǎng)格數(shù)量變得十分龐大，在沒有集群機(jī)或在集群機(jī)計(jì)算條件受限情況下，單機(jī)計(jì)算時(shí)間過長，難以在較短時(shí)間內(nèi)獲得計(jì)算結(jié)果，在一定程度上限制了非結(jié)構(gòu)化網(wǎng)格水動力模型在工程領(lǐng)域上的應(yīng)用.

GPU（graphic processing unit，圖形處理器）與CPU（central processing unit，中央處理器）的計(jì)算能力的差別主要源于二者的計(jì)算架構(gòu)不同，一般的處理器由控制單元、邏輯單元和存儲單元3個(gè)部分組成.在CPU 中，控制單元和存儲單元占的份額較大，而GPU 中擁有較多的邏輯運(yùn)算單元，因此相對于CPU，GPU 擁有更強(qiáng)的浮點(diǎn)運(yùn)算能力.以本文使用的GTX460 顯卡為例，其擁有336個(gè)流處理器，核心頻率為700 MHz，顯存容量為1GB，浮點(diǎn)運(yùn)算能力達(dá)到0.961TFLOPS.

早期的GPU 僅僅用來完成圖形處理工作，但隨著技術(shù)的進(jìn)步，對GPU 的可編程性增強(qiáng)，其能力已經(jīng)超出了原先設(shè)計(jì)的功能，尤其是近年來，圖形處理器硬件水平得到高速發(fā)展［1］，更多的開發(fā)者將其強(qiáng)大的并行計(jì)算能力應(yīng)用到通用計(jì)算領(lǐng)域.目前GPU 已經(jīng)應(yīng)用于n－body模擬、醫(yī)學(xué)CT圖像重建、計(jì)算流體力學(xué)、分子動力學(xué)模擬、金融計(jì)算和油／氣藏模擬等方面.2003年Li等［2］利用API環(huán)境實(shí)現(xiàn)了格子玻爾茲曼方法在GPU 上的計(jì)算.2011年Aaron把GPU 應(yīng)用在利用大渦模擬解決基于結(jié)構(gòu)化網(wǎng)格的湍流問題上，選用的CPU 是AMD Phenom 2.6GHz四核處理器，在Linux x64操作系統(tǒng)上運(yùn)行，最多將計(jì)算速度提高了15 倍左右.2011年朱小松［3］實(shí)現(xiàn)了基于GPU 并行加速的粒子法解決液體晃蕩問題，加速能力最多提高到26倍.

總的來說，目前基于GPU 的流體計(jì)算，多采用粒子法求解拉格朗日型的N－S方程，而對于求解歐拉型N－S方程，尤其是非結(jié)構(gòu)網(wǎng)格下，利用GPU 解決流體動力學(xué)問題的研究很少.針對以上問題，本文采用CUDA 編程模型，基于GPU 架構(gòu)并行算法建立非結(jié)構(gòu)化網(wǎng)格水動力模型，并對所建模型的計(jì)算效率與集群機(jī)的性能進(jìn)行對比分析.

CUDA平臺使用帶有關(guān)鍵詞擴(kuò)展的標(biāo)準(zhǔn)C開發(fā)語言，可以通過類C語言編寫在顯卡芯片上執(zhí)行的程序.CUDA 編程模型中，將CPU 作為主機(jī)（Host），將GPU 作為協(xié)處理器或設(shè)備（Device）［4］，二者共同協(xié)作完成整個(gè)程序.其中，Host端負(fù)責(zé)執(zhí)行串行部分的程序，Device端負(fù)責(zé)并行部分，Device端的程序又稱為“kernel”.GPU 執(zhí)行時(shí)的最小單位是thread，CUDA 架構(gòu)可以擁有海量的thread，非常適合大規(guī)模數(shù)據(jù)并行計(jì)算［5］.

GPU 并行計(jì)算相對于傳統(tǒng)的CPU 有以下幾個(gè)優(yōu)勢：成本相對較低，一個(gè)普通計(jì)算機(jī)的價(jià)格，計(jì)算性能達(dá)到了TFLOPS，相當(dāng)于一個(gè)高性能計(jì)算集群系統(tǒng)［6］；具有可擴(kuò)展性，可以多個(gè)GPU 協(xié)同并行計(jì)算，支持高效線程，使得應(yīng)用程序可以提供比現(xiàn)有硬件執(zhí)行資源更高的并行度［7］；節(jié)點(diǎn)之間可以共享存儲器數(shù)據(jù)，顯存帶寬是內(nèi)存的10倍左右，可以高速傳遞數(shù)據(jù)；有眾多的核心，完全可以使用胖節(jié)點(diǎn)式的并行算法，而CPU 的大型的并行，都是在集群層面的，而不是在單機(jī)上的；CUDA 這一類的GPU 并行編程語言，并行層面在線程級別，相對于CPU 要方便許多.

1 非結(jié)構(gòu)化網(wǎng)格水動力模型

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

從N－S方程出發(fā)，假設(shè)海水為不可壓縮黏性流體，密度恒定，以連續(xù)性方程和動量方程作為二維水動力數(shù)值模型的控制方程：

底部剪切力

固邊界條件

其中（u，v）為水平面流速，（τsx，τsy）為表面風(fēng)應(yīng)力，（τbx，τby）為底部剪切應(yīng)力，f為科氏力系數(shù)，η為水面高程，M為曼寧系數(shù)，C為謝才系數(shù).根據(jù)Smagorinsky渦黏參數(shù)法［8］，水平黏性項(xiàng)可以由下式近似表示：

1.2 空間離散方法

本文使用的模型選用非結(jié)構(gòu)化三角形網(wǎng)格，采用有限體積法進(jìn)行數(shù)值求解，在求解過程中將整個(gè)數(shù)值計(jì)算區(qū)域劃分為相互不重疊的非結(jié)構(gòu)化三角形網(wǎng)格，如圖1所示.

圖1 數(shù)值求解離散化非結(jié)構(gòu)化三角形網(wǎng)格示意圖Fig.1 Illustration of numerically solving the discretized unstructured triangular grid

三角形單元τ中的變量為線性分布，滿足二階空間精度，如式（11）～（15）所示：

1.3 時(shí)間離散方法

在本文的數(shù)學(xué)模型中對連續(xù)性方程（1）求面積積分，通過修正過的四階R－K 時(shí)間積分進(jìn)行求解運(yùn)算，如式（16）～（18）所示：

2 基于GPU并行算法的實(shí)現(xiàn)及優(yōu)化

目前CUDA 的應(yīng)用程序大多為基于C 語言進(jìn)行開發(fā)的，程序?qū)崿F(xiàn)相對容易.而現(xiàn)在廣泛使用的海洋模型如POM、ECOM、FVCOM 等均為使用FORTRAN 語言進(jìn)行開發(fā)的，雖然可以在原有的模型中調(diào)用CUDA 計(jì)算核心，但是并行計(jì)算部分的設(shè)計(jì)，與重新建立模型的工作量相差不大，計(jì)算效率也略低于C程序，因此本文以C 語言為開發(fā)語言，在CUDA 平臺下建立基于GPU 并行計(jì)算的水動力模型.

在基于CPU 并行計(jì)算水動力模型中，由于CPU 核心的數(shù)量遠(yuǎn)小于網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)，傳統(tǒng)的方法是利用美國密西根大學(xué)Karypis等編寫的用于圖的分區(qū)和稀疏矩陣排序的串行包——METIS庫，將計(jì)算區(qū)域分為N個(gè)子區(qū)域，分配給N個(gè)CPU 進(jìn)行計(jì)算，把子區(qū)域的初始流場信息、幾何信息分別存儲到對應(yīng)的CPU 內(nèi)存中，在每一個(gè)CPU中啟動計(jì)算進(jìn)程，由主進(jìn)程調(diào)度各CPU 的計(jì)算［9］.在這一模式下，計(jì)算速度受到了多個(gè)方面因素的影響，包括CPU 的計(jì)算速度、內(nèi)存帶寬、在分區(qū)的邊界處須保持CPU間數(shù)據(jù)交換過程等，因此并行計(jì)算的效率并不隨CPU 核心的個(gè)數(shù)呈線性增加.

由于GPU 的并行計(jì)算架構(gòu)與CPU 的并行計(jì)算架構(gòu)有較大的區(qū)別，在并行化的實(shí)現(xiàn)過程中也有很大的不同.GPU 中的計(jì)算核心稱為流處理器（SP），以本文中使用的GTX460 顯卡為例，包含了330個(gè)流處理器.GPU 通過同時(shí)往每個(gè)流處理器上并發(fā)執(zhí)行線程來實(shí)現(xiàn)加速計(jì)算.若像傳統(tǒng)的并行方式一樣利用METIS庫將計(jì)算區(qū)域分塊計(jì)算，一方面如果分塊過多會導(dǎo)致計(jì)算過程過于復(fù)雜，另一方面并行程序會因?yàn)槊總€(gè)分塊之間進(jìn)行數(shù)據(jù)交換而導(dǎo)致并行計(jì)算效率的下降，因此在GPU 并行程序設(shè)計(jì)中盡可能多地分配線程，每個(gè)線程執(zhí)行對一個(gè)網(wǎng)格或網(wǎng)格節(jié)點(diǎn)的計(jì)算.

2.1 GPU 并行算法的實(shí)現(xiàn)

在GPU 上實(shí)現(xiàn)并行算法需要以下幾個(gè)步驟：

（1）執(zhí)行GPU 并行計(jì)算前檢測設(shè)備，根據(jù)不同的硬件設(shè)備，設(shè)置執(zhí)行環(huán)境，并在CPU 上完成程序的前處理.

（2）在顯卡上分配存儲空間，用cudaMemcpy函數(shù)將在文件中讀取的初始數(shù)據(jù)發(fā)送到顯存中.

（3）將串行計(jì)算程序中需要每次對網(wǎng)格循環(huán)計(jì)算的程序段用對應(yīng)的GPU 并行子程序?qū)崿F(xiàn)，在每個(gè)并行子程序中設(shè)定在GPU 并行計(jì)算中的Grid和Block的維度.

（4）輸出計(jì)算結(jié)果，用cudaMemcpy函數(shù)將顯存中的數(shù)據(jù)返回到內(nèi)存中，再輸出到結(jié)果文件.

2.2 對GPU 并行算法的優(yōu)化

在GPU 上實(shí)現(xiàn)并行計(jì)算，核心技術(shù)是算法的并行化和程序優(yōu)化.在本文中按照計(jì)算模型的特點(diǎn)進(jìn)行了適當(dāng)?shù)膬?yōu)化.

在GPU 并行程序中處理的對象是空間離散后的每個(gè)網(wǎng)格以及每個(gè)網(wǎng)格節(jié)點(diǎn)，這樣就避免了由于分區(qū)不均衡引起的計(jì)算效率降低，同時(shí)還節(jié)省了分區(qū)邊界處數(shù)據(jù)交換的時(shí)間.因此在原串行程序中的每次對網(wǎng)格數(shù)組或節(jié)點(diǎn)數(shù)組進(jìn)行操作中，將其每次循環(huán)需要計(jì)算的內(nèi)容對應(yīng)地分配到GPU 的各個(gè)線程中，以海量的線程個(gè)數(shù)隱藏了數(shù)據(jù)訪問延遲的時(shí)間.

由于在本文中使用的是非結(jié)構(gòu)化網(wǎng)格，無法按照矩陣方式分配在每個(gè)線程中計(jì)算的網(wǎng)格，故在GPU 并行程序的線程設(shè)計(jì)中將Block設(shè)置為一維，并且每個(gè)Block中包含64個(gè)線程，即兩個(gè)線程束.則整體的Block的個(gè)數(shù)為N／64＋1（其中N為計(jì)算單元個(gè)數(shù)），這樣就將GPU 劃分為（N／64＋1）×64個(gè)線程，這種線程分配方式最大限度地降低了分支流程和空載線程的個(gè)數(shù)，實(shí)現(xiàn)了GPU 核心的最高利用效率.根據(jù)式（1）～（16），每個(gè)GPU 的線程計(jì)算一個(gè)網(wǎng)格或網(wǎng)格節(jié)點(diǎn)的數(shù)據(jù)，將所需要的數(shù)據(jù)和需要計(jì)算的結(jié)果作為kernel函數(shù)的參數(shù)，并在所有節(jié)點(diǎn)的數(shù)據(jù)計(jì)算完成后，將計(jì)算結(jié)果返回到顯存的數(shù)組中.

為了減少顯存和內(nèi)存之間數(shù)據(jù)傳輸速度的限制，在整個(gè)計(jì)算過程中，完成初始化后將變量傳遞到顯存中，在迭代計(jì)算過程中，所有數(shù)據(jù)均在顯存中進(jìn)行交換，直至需要將中間結(jié)果輸出，再將顯存中的數(shù)據(jù)返回到內(nèi)存的數(shù)組中，這樣減少了由于數(shù)據(jù)交換引起的時(shí)間損耗.

3 模型驗(yàn)證

數(shù)值實(shí)驗(yàn)的模型水深D＝10 m，長L＝3.7km，寬B＝1.5km，整個(gè)流場的右側(cè)為開邊界，有作用振幅為0.5m 的M2分潮.為了避免受到邊界影響，在該數(shù)值實(shí)驗(yàn)結(jié)果對比中，驗(yàn)證點(diǎn)選擇整個(gè)流場的中間點(diǎn)，距離左側(cè)邊界1.85 km［10］.

根據(jù)Ippen 的理論［11］，對于此種水域情況，有如下解析解：

其中A為振幅，ω為角頻率，x為該點(diǎn)距離左側(cè)邊界 的 距 離，h為 水 深，T為 潮 周 期，η為 潮 波 面 升高，u為縱向流速.

在數(shù)值模型中將計(jì)算區(qū)域剖分成1 620個(gè)三角形網(wǎng)格，驗(yàn)證點(diǎn)潮位和縱向流速計(jì)算結(jié)果歷時(shí)曲線如圖2所示.從圖中可以看出，數(shù)值計(jì)算結(jié)果與解析解吻合度較好.

圖2 驗(yàn)證點(diǎn)潮位及流速歷時(shí)曲線計(jì)算結(jié)果對比Fig.2 The contrast of elevation and current speed duration curve result in the verified points

4 GPU 算法的并行效率分析

為了更好地說明GPU 在水動力模型并行計(jì)算過程中起到的作用，數(shù)值模型實(shí)驗(yàn)采用在上一章模型驗(yàn)證中的計(jì)算水域，并設(shè)計(jì)了4種不同網(wǎng)格密度的工況算例，分別使用CPU 單線程（Intel CoreTMi7 930 2.8GHz）、集群機(jī)（使用30核心、24核心、20核心，最大計(jì)算節(jié)點(diǎn)數(shù)為12，CPU 型號為Intel Xeon E5430，主頻2.66 GHz），以及GTX460顯卡（330 個(gè)流處理器，核心頻率700 MHz）進(jìn)行計(jì)算，并對計(jì)算耗用時(shí)間進(jìn)行對比分析.在集群機(jī)計(jì)算過程中發(fā)現(xiàn)，在工況1～3中，網(wǎng)格數(shù)不多時(shí)，總的計(jì)算核心數(shù)保持不變，集群機(jī)節(jié)點(diǎn)數(shù)以及每個(gè)節(jié)點(diǎn)使用核心數(shù)變化對計(jì)算時(shí)間影響不大，計(jì)算所需時(shí)間基本不變；在工況4中，網(wǎng)格數(shù)超過20萬，受集群機(jī)節(jié)點(diǎn)數(shù)以及每個(gè)節(jié)點(diǎn)使用核心數(shù)影響明顯，不同節(jié)點(diǎn)的分配方式得到的加速比R如圖3所示，其中加速比為單線程計(jì)算時(shí)間與集群機(jī)計(jì)算使用時(shí)間的比值，橫坐標(biāo)為節(jié)點(diǎn)數(shù)×每個(gè)節(jié)點(diǎn)使用核心數(shù).

由圖3可以看出，在核心總數(shù)相同的條件下，使用的節(jié)點(diǎn)數(shù)越多，計(jì)算效率越高，相應(yīng)的加速比就越高.因此集群機(jī)的計(jì)算效率受計(jì)算核心總數(shù)，以及使用節(jié)點(diǎn)個(gè)數(shù)的影響比較明顯.

采用集群機(jī)的最優(yōu)化的節(jié)點(diǎn)分配方案，并與GPU（Nvidia GTX460）并行計(jì)算的加速能力進(jìn)行對比如表1和圖4所示.

圖3 集群機(jī)在不同節(jié)點(diǎn)分配條件下計(jì)算工況4得到的加速比Fig.3 The speedup ratio in different distributions of compute nodes when using cluster to compute Case 4

表1 集群機(jī)和GPU 不同網(wǎng)格數(shù)條件下所用時(shí)間Tab.1 The time of cluster and GPU cost in different number of grids

圖4 不同網(wǎng)格數(shù)條件下并行加速能力對比Fig.4 The contrast of speedup capacity at different number of grids

表1中結(jié)果為分別使用單線程CPU 串行算法、3種不同計(jì)算核心數(shù)的MPI并行算法，以及GPU 并行算法，計(jì)算4種不同網(wǎng)格數(shù)下算例所使用時(shí)間.圖4為表1中各并行算法所對應(yīng)的加速比.在本文中最多只使用了12個(gè)計(jì)算節(jié)點(diǎn)，由圖4可以看出，隨著網(wǎng)格個(gè)數(shù)的增加，GPU 并行計(jì)算的加速比能夠?qū)崿F(xiàn)相對穩(wěn)定的提高.計(jì)算網(wǎng)格數(shù)在工況1～3 時(shí)，加速比提高的速度低于集群機(jī)，網(wǎng)格數(shù)在工況3～4時(shí)，由于節(jié)點(diǎn)使用個(gè)數(shù)的限制，使用30 核和20 核的加速比略有下降，GPU 的加速比仍穩(wěn)定地提高，其加速比提高的速度基本與24核集群機(jī)的平行.由加速比隨網(wǎng)格數(shù)變化曲線的趨勢可以看出，若網(wǎng)格數(shù)繼續(xù)增加，集群機(jī)的加速能力會受到節(jié)點(diǎn)分配方案和等待其他分區(qū)計(jì)算完成時(shí)間的影響而有所降低，而GPU加速能力仍會相對穩(wěn)定地提高.

5 結(jié) 論

（1）對不同計(jì)算核心分配方案的MPI并行計(jì)算與單CPU 串行計(jì)算性能進(jìn)行對比分析，得到結(jié)論：當(dāng)網(wǎng)格節(jié)點(diǎn)數(shù)在工況1～3（即網(wǎng)格數(shù)較少）時(shí)，不同的節(jié)點(diǎn)分配方案對集群機(jī)的加速能力影響不大；在網(wǎng)格數(shù)量較多時(shí)，集群機(jī)加速能力受計(jì)算節(jié)點(diǎn)分配方式影響顯著.

（2）通過對不同計(jì)算核心數(shù)的MPI并行算法以及GPU 并行算法的加速比的對比得出：使用GPU 并行算法，能夠在保證較高計(jì)算精度的基礎(chǔ)上，大幅提高方程求解的速度，尤其是隨著網(wǎng)格數(shù)的增加，基于GPU 求解方程的計(jì)算加速比得到顯著的提高.當(dāng)網(wǎng)格個(gè)數(shù)增加到20萬個(gè)以上時(shí)，GPU 的計(jì)算能力已經(jīng)接近使用20核集群機(jī)的計(jì)算能力，加速比提高了10倍以上，并且加速比的變化率已經(jīng)超過20核和30核的集群機(jī)，并行加速效果十分顯著.若考慮其他方面的計(jì)算成本，如集群機(jī)節(jié)點(diǎn)的利用率、能耗等方面，使用GPU 進(jìn)行并行計(jì)算的優(yōu)勢更加明顯.由此可見，GPU 在大范圍海域的水動力模型的數(shù)值計(jì)算方面具有較好的發(fā)展前景.

［1］ NVIDIA.NVIDIA CUDA C programming guide［EB／OL］.［2013－08－19］.https：／／www.nvidia.com.

［2］ LI Wei，WEI Xiao－ming，Kaufman A.Implementing lattice Boltzmann computation on graphics hardware ［J］.Visual Computer，2003，19（7－8）：444－456.

［3］ 朱小松.粒子法的并行加速及在液體晃蕩研究中的應(yīng)用［D］.大連：大連理工大學(xué)，2011.ZHU Xiao－song.Parallel acceleration of particle method and its application on the study of liquid sloshing ［D］.Dalian：Dalian University of Technology，2011.（in Chinese）

［4］ 張 舒，褚艷麗.GPU 高性能運(yùn)算之CUDA［M］.北京：中國水利水電出版社，2009.ZHANG Shu，CHU Yan－li.GPU High Performance Computing：CUDA［M］.Beijing：China Water Power Press，2009.（in Chinese）

［5］ 王英俊，王啟富，王 鋼，等.CUDA 架構(gòu)下的三維彈性靜力學(xué)邊界元并行計(jì)算［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2012，24（1）：112－119.WANG Ying－jun，WANG Qi－fu，WANG Gang，etal.CUDA based parallel computation of BEM for 3Delastostatics problems［J］.Journal of Computer－Aided Design ＆Computer Graphics，2012，24（1）：112－119.（in Chinese）

［6］ 多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室.基于GPU 的多尺度離散模擬并行計(jì)算［M］.北京：科學(xué)出版社，2009.State Key Laboratory of Multiphase Complex Systems.Multi－scale Discrete Simulation Parallel Computing Based on GPU ［M］.Beijing：Science Press，2009.（in Chinese）

［7］ Kirk D B，Hwu Wen－mei W.Programming Massively Parallel Processors［M］.Beijing：Tsinghua University Press，2010.

［8］ 陳長勝.海洋生態(tài)系統(tǒng)動力學(xué)與模型［M］.北京：高等教育出版社，2004.CHEN Chang－sheng.Ocean Ecosystem Dynamics and the Model ［M］.Beijing：Higher Education Press，2004.（in Chinese）

［9］ 王 敏，王 明，楊 明.小浪底水庫三維數(shù)學(xué)模型并行計(jì)算研究［J］.人民黃河，2012，34（5）：25－27.WANG Min，WANG Ming，YANG Ming.Parallel computing research of Xiaolangdi Reservoir threedimensional mathematical model［J］.Yellow River，2012，34（5）：25－27.（in Chinese）

［10］ 張瑞瑾.三維潮流數(shù)值模型及其應(yīng)用［D］.大連：大連理工大學(xué)，2002.ZHANG Rui－jin.Establish and the application of a three dimensional numerical tidal model ［D］.Dalian：Dalian University of Technology，2002.（in Chinese）

［11］ Ippen A T.Estuary and Coastline Hydrodynamics［M］.New York：McGraw－Hill Book Co，1966.