朱煉華，郭照立

（華中科技大學(xué)煤燃燒國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430074）

文章編號(hào)：1001?246X（2015）01?0020?07

基于格子Boltzmann方法的多孔介質(zhì)流動(dòng)模擬GPU加速

朱煉華，郭照立

（華中科技大學(xué)煤燃燒國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430074）

利用NVIDIA CUDA平臺(tái)，在GPU上結(jié)合稀疏存貯算法實(shí)現(xiàn)基于格子Boltzmann方法的孔隙尺度多孔介質(zhì)流動(dòng)模擬加速，測(cè)試該算法相對(duì)基本算法的性能.比較該算法在不同GPU上使用LBGK和MRT兩種碰撞模型及單、雙精度計(jì)算時(shí)的性能差異.測(cè)試結(jié)果表明在GPU環(huán)境下采用稀疏存貯算法相對(duì)基本算法能大幅提高計(jì)算速度并節(jié)省顯存，相對(duì)于串行CPU程序加速比達(dá)到兩個(gè)量級(jí).使用較新構(gòu)架的GPU時(shí)，MRT和LBGK碰撞模型在單、雙浮點(diǎn)數(shù)精度下計(jì)算速度相同.而在較上一代的GPU上，計(jì)算精度對(duì)MRT碰撞模型計(jì)算速度影響較大.

多孔介質(zhì)；GPU；格子Boltzmann方法；并行計(jì)算

0 引言

多孔介質(zhì)內(nèi)的流動(dòng)現(xiàn)象廣泛存在于自然界及油氣藏開(kāi)采、化工過(guò)程、環(huán)境保護(hù)等諸多工業(yè)領(lǐng)域.多孔介質(zhì)流動(dòng)是一類典型的跨尺度現(xiàn)象，在孔隙尺度上利用計(jì)算流體力學(xué)方法對(duì)這類現(xiàn)象進(jìn)行直接模擬是研究其微觀流動(dòng)機(jī)理的重要手段.由于多孔介質(zhì)的孔隙結(jié)構(gòu)十分復(fù)雜，使用常規(guī)的基于直接離散求解Navier?Stoke方程的方法（有限體積法、有限差分法、有限元法）在孔隙尺度模擬時(shí)流固邊界難以處理.格子Boltzmann方法（LBM）作為一類基于微觀分子動(dòng)力學(xué)的介觀數(shù)值方法，復(fù)雜流固邊界容易處理，并且并行效率高，十分適合孔隙尺度多孔介質(zhì)流動(dòng)的并行數(shù)值模擬［1－2］.

模擬孔隙尺度多孔介質(zhì)流動(dòng)時(shí)，刻畫高度非規(guī)則的流固邊界需要很高網(wǎng)格分辨率，因此計(jì)算量很大，通常需要采用并行計(jì)算加速.最近幾年流行起來(lái)的通用GPU（GPGPU）高性能技術(shù)是一種非常有潛力的并行計(jì)算技術(shù).現(xiàn)代GPU提供的計(jì)算能力和存儲(chǔ)帶寬遠(yuǎn)遠(yuǎn)超過(guò)同期主流的CPU.NVIDIA公司于2007年推出CUDA GPU計(jì)算構(gòu)架后，GPGPU編程難度降低，諸多科學(xué)計(jì)算應(yīng)用被移植到GPU平臺(tái)［3－4］.LBM算法空間局部性優(yōu)良并且計(jì)算瓶頸在訪存帶寬，因此特別適合于GPU并行計(jì)算.T?lke等人最早實(shí)現(xiàn)了LBM的CUDA程序，在單塊GPU加速比達(dá)到一個(gè)量級(jí)以上［5］.國(guó)內(nèi)黃昌盛等人系統(tǒng)分析了LBM的GPU優(yōu)化方法［6］，李博等人做了關(guān)于多GPU的LBM實(shí)現(xiàn)方面的工作［7］.

目前關(guān)于GPU加速LBM計(jì)算的文章中，流場(chǎng)大多為較規(guī)則的空腔體，而LBM的典型應(yīng)用領(lǐng)域則是具有復(fù)雜流固邊界的流動(dòng)模擬.針對(duì)規(guī)則流場(chǎng)的常規(guī)LBM GPU算法在模擬多孔介質(zhì)這種高度非規(guī)則流場(chǎng)時(shí)，存在顯存浪費(fèi)、GPU執(zhí)行效率不高等問(wèn)題.另一方面，在CPU平臺(tái)，Pan等人在2004年提出了一種針對(duì)多孔介質(zhì)流動(dòng)的稀疏存貯LBM算法，可以大幅降低內(nèi)存耗用［8］，但該工作使用的是傳統(tǒng)的基于消息傳遞接口（MPI）的并行方式.鑒于目前GPU上搭載的顯存容量非常有限，在GPU上應(yīng)用該算法顯得尤為必要.由于GPU構(gòu)架與CPU相差較大，該算法的GPU版本性能受更多因素影響.本文將上述稀疏存貯算法應(yīng)用于多孔介質(zhì)流動(dòng)模擬的GPU LBM程序，測(cè)試并分析該算法相對(duì)基本算法的性能優(yōu)勢(shì)，然后比較兩種常見(jiàn)碰撞模型（MRT、LBGK）以及所使用的浮點(diǎn)數(shù)精度對(duì)該算法性能的影響.

1 計(jì)算模型

模擬外力驅(qū)動(dòng)下多孔介質(zhì)流動(dòng)問(wèn)題，使用標(biāo)準(zhǔn)D3Q19模型，其分布函數(shù)（DF）fi的演化方程為

其中下標(biāo)i表示離散速度方向，x為空間位置，t為時(shí)間，δt為時(shí)間步長(zhǎng)，ci為格子離散速度，其各分量為（2）

方程（1）右端第一項(xiàng)Ωi為碰撞算子，第二項(xiàng)Fi為離散外力項(xiàng).本文使用的He［9］等人提出的外力離散格式，即

其中cs＝1/3為格子聲速， 是無(wú)量綱松弛時(shí)間，與運(yùn)動(dòng)學(xué)粘性ν有關(guān)，ν＝c2s（－0.5）δt，F(xiàn)為外力大小，ρ和u分別為宏觀密度和速度，（x，t）為平衡態(tài)分布函數(shù)，其定義為

其中ωi為權(quán)系數(shù)，ω0＝1/3，ω1，…，6＝1/18，ω7，…，18＝1/36．流體宏觀密度ρ和速度u由密度分布函數(shù)得到

使用兩種碰撞算子模型，即單松弛時(shí)間（LBGK）模型和多松弛時(shí)間（MRT）模型.LBGK模型碰撞算子形式為

MRT模型的碰撞過(guò)程在密度分布函數(shù)的速度矩空間完成，其碰撞算子形式為

其中M是變換矩陣，向量meq是速度矩的平衡態(tài)，M和meq的具體形式可參考文獻(xiàn)［10］．矩陣S為對(duì)角矩陣，其對(duì)角元素分別為各個(gè)矩的松弛因子，即

其中sν和sε與運(yùn)動(dòng)學(xué)粘度ν和體粘度ζ相關(guān)：

sν由運(yùn)動(dòng)學(xué)粘度ν控制，其它松弛因子按Ginzburg等人提出的雙松弛時(shí)間（TRT）模型［11］取值，即

2 LBM算法的GPU實(shí)現(xiàn)及分析

2.1 CUDA構(gòu)架

CUDA構(gòu)架通過(guò)一種Grid?Block?Thread的三層次結(jié)構(gòu)將計(jì)算任務(wù)按數(shù)據(jù)并行的方式映射到GPU上的大量處理單元上并行執(zhí)行，同一個(gè)Block中的Thread之間可以通過(guò)共享內(nèi)存通信.并行計(jì)算任務(wù)被編寫為稱之為Kernel的函數(shù).為充分利用GPU提供的帶寬，訪問(wèn)顯存時(shí)需要滿足連續(xù)、合并訪問(wèn).為高效利用GPU上的大量處理單元，需要合理選擇Block大小并控制Kernel函數(shù)中使用的寄存器個(gè)數(shù)，并且要避免Kernel函數(shù)中出現(xiàn)分支判斷語(yǔ)句.

2.2 基本LBM算法的GPU實(shí)現(xiàn)

LBM算法是一種顯式時(shí)間推進(jìn)算法，在每個(gè)時(shí)間步內(nèi)，依次執(zhí)行碰撞、遷移和邊界處理三個(gè)步驟.多孔介質(zhì)中有大量不參與計(jì)算的固體區(qū)域，可以用一個(gè)三維Flag數(shù)組來(lái)標(biāo)記格子中格點(diǎn)類型，演化時(shí)只有對(duì)流體格點(diǎn)執(zhí)行碰撞和遷移步驟.為減少全局顯存訪問(wèn)次數(shù)，碰撞和遷移可以放在一步執(zhí)行.為了實(shí)現(xiàn)連續(xù)、合并訪問(wèn)全局顯存，應(yīng)保證相鄰格點(diǎn)的同一離散速度的DF在顯存中的位置是連續(xù)的.以C語(yǔ)言為例，開(kāi)辟的DF數(shù)組形式應(yīng)為f［Q］［NZ］［NY］［NX］，這里Q表示離散速度個(gè)數(shù)，NX、NY、NZ表示格子大?。捎谶w移步每個(gè)格點(diǎn)需要訪問(wèn)周圍格點(diǎn)DF，當(dāng)遷移的速度方向有X方向分量時(shí)會(huì)出現(xiàn)內(nèi)存訪問(wèn)不對(duì)齊的情況，T?lke最早提出利用GPU上的共享內(nèi)存輔助遷移過(guò)程來(lái)克服這個(gè)問(wèn)題［5］.隨后Obrecht［12］等人指出這種方式在實(shí)現(xiàn)復(fù)雜Kernel函數(shù)時(shí)會(huì)因共享內(nèi)存數(shù)量不足而限制GPU執(zhí)行效率，而且在較新的GPU上，對(duì)顯存訪問(wèn)地址對(duì)齊的要求大為降低.本文在實(shí)現(xiàn)時(shí)采用Obrecht的方法，即直接訪問(wèn)周圍格點(diǎn)DF.

在本文所考慮的多孔介質(zhì)流動(dòng)模擬中，所有外邊界都是用周期邊界條件，而孔隙壁面則用標(biāo)準(zhǔn)反彈格式處理，這兩種邊界條件都可以在遷移過(guò)程中隱式處理，不需要增加單獨(dú)的Kernel函數(shù)來(lái)特殊處理.總結(jié)以上討論，基本算法在每其一個(gè)演化步的計(jì)算流程為

1）通過(guò)Flag數(shù)組判斷當(dāng)前格點(diǎn)類型，若當(dāng)前為固體格點(diǎn)，則不作任何處理，否則進(jìn)入下一步；

2）讀取周圍格點(diǎn)DF，若某一方向相鄰格點(diǎn)為固體格點(diǎn)，則從當(dāng)前格點(diǎn)反方向讀取DF；

3）統(tǒng)計(jì)宏觀量；

4）執(zhí)行碰撞操作，將新的DF寫入顯存.

2.3 稀疏存儲(chǔ)算法

上一小節(jié)描述的通過(guò)Flag數(shù)組來(lái)刻畫流場(chǎng)結(jié)構(gòu)的方法雖然實(shí)現(xiàn)簡(jiǎn)單，但在GPU上實(shí)現(xiàn)時(shí)存在一些問(wèn)題.首先，在開(kāi)辟DF數(shù)組時(shí)，不需要演化固體格點(diǎn)也占據(jù)了存儲(chǔ)空間和帶寬.自然界中的多孔介質(zhì)孔隙率通常在0.3～0.4之間［8］，相當(dāng)于有多達(dá)70%的顯存被浪費(fèi)了.目前主流的GPU所搭載的顯存容量有限，并且其存儲(chǔ)帶寬是其性能瓶頸，所以在GPU上這一問(wèn)題也更突出；其次，由于演化時(shí)對(duì)固體格點(diǎn)和流體格點(diǎn)處理方式不同，導(dǎo)致線程執(zhí)行過(guò)程分支，GPU執(zhí)行效率大為降低.

為解決上述問(wèn)題，可以借鑒稀疏矩陣存貯方法.稀疏矩陣中大量元素都是0，存貯時(shí)可以只儲(chǔ)非零元素及其位置.同樣，在處理多孔介質(zhì)時(shí)可以只存貯流體格點(diǎn)的DF，即事先將流體格點(diǎn)從多孔介質(zhì)結(jié)構(gòu)中提取出來(lái)，然后開(kāi)辟形式為F［fluid_num］［Q］的數(shù)組存儲(chǔ)其DF，其中fluid_num為多孔結(jié)構(gòu)中流體格點(diǎn)的個(gè)數(shù)．采用這種存儲(chǔ)方案時(shí)，格點(diǎn)之間的相對(duì)位置關(guān)系丟失了，所以還需要構(gòu)造一個(gè)輔助數(shù)組node_map［fluid_num］［Q－1］來(lái)記錄每個(gè)格點(diǎn)的相鄰Q－1個(gè)格點(diǎn)在數(shù)組F中的索引位置．在格點(diǎn)i執(zhí)行k方向遷移操作前，先通過(guò)這個(gè)輔助數(shù)組查詢其k方向的相鄰格點(diǎn)在數(shù)組F中的位置（node_map［i］［k］），再執(zhí)行寫入操作，完成遷移過(guò)程.最后，在模擬計(jì)算結(jié)束并統(tǒng)計(jì)出宏觀量后，再將其還原到原始的三維矩陣結(jié)構(gòu)用于后處理.

可以分析對(duì)比基本算法和稀疏存貯算法的顯存耗用情況.假設(shè)計(jì)算區(qū)域流固格點(diǎn)總數(shù)為N，多孔結(jié)構(gòu)的孔隙率為ε，DF使用單精度浮點(diǎn)數(shù)．node_map元素使4字節(jié)無(wú)符號(hào)整形變量存儲(chǔ).若不采用稀疏存儲(chǔ)模式，忽略Flag數(shù)組所占空間，則需要消耗顯存

若采用稀疏存儲(chǔ)模式，則需要消耗顯存

可以發(fā)現(xiàn)當(dāng)孔隙率ε＜0.75時(shí)，M2＜M1．因此在大部分情況下，采用稀疏存貯均能節(jié)省顯存開(kāi)銷.

采用稀疏存儲(chǔ)方案時(shí)，格點(diǎn)之間的拓?fù)潢P(guān)系已經(jīng)存儲(chǔ)在node_map數(shù)組里面了，在構(gòu)造這個(gè)數(shù)組時(shí)，周期邊界格點(diǎn)和孔壁處反彈格點(diǎn)的遷移規(guī)則均可反映在這個(gè)數(shù)組里，即不用顯式判斷并分別處理外邊界格點(diǎn)和孔隙壁面反彈格點(diǎn)，演化的Kernel函數(shù)里面不會(huì)出現(xiàn)判斷分支語(yǔ)句.因此，可以預(yù)測(cè)采用稀疏存貯方案后GPU執(zhí)行效率會(huì)有顯著提高.

3 數(shù)值算例

3.1 驗(yàn)證算例

本節(jié)先通過(guò)模擬計(jì)算得出一種理想多孔介質(zhì)模型的絕對(duì)滲透率并與解析解對(duì)比來(lái)驗(yàn)證GPU程序?qū)崿F(xiàn)的正確性，然后在此基礎(chǔ)上，測(cè)試稀疏存儲(chǔ)算法與基本算法的計(jì)算速度，最后測(cè)試分析影響計(jì)算性能的其它因素，如LBGK/MRT模型、單/雙精度.

計(jì)算使用的多孔介質(zhì)模型——體心立方（BCC）結(jié)構(gòu)的最小單元如圖1所示，位于中心的球體和立方體頂點(diǎn)上的1/8半球體半徑均為r，立方體邊長(zhǎng)為L(zhǎng).固體區(qū)域所占比例和孔隙率分別為

其中c的最大值為cmax＝3π/8．BCC結(jié)構(gòu)的絕對(duì)滲透率近似解析解形式為

其中d?為單個(gè)球體所受的無(wú)量綱阻力，在Re＜1的時(shí)，其級(jí)數(shù)解析解為

圖1 BCC多孔介質(zhì)模型Fig.1 BCC porousmediamodel

數(shù)值模擬時(shí)，通過(guò)給定外力g，在流場(chǎng)達(dá)到穩(wěn)定后，統(tǒng)計(jì)流場(chǎng)平均流速ud，通過(guò)達(dá)西滲流定律可確定多孔介質(zhì)的滲透率

其中ud為流場(chǎng)平均流速，▽pz和ρg分別為壓力梯度和外力大小.

圖2 使用不同碰撞模型的計(jì)算結(jié)果Fig.2 LBGK and MRT results

Pan等人指出LBGK模型結(jié)合標(biāo)準(zhǔn)反彈格式模擬多孔介質(zhì)流動(dòng)時(shí)，計(jì)算所得的滲透率與流體運(yùn)動(dòng)粘度相關(guān)聯(lián)，這是非物理的現(xiàn)象，而若使用MRT模型，則可顯著改善這一問(wèn)題［10］.本文分別使用MRT和LBGK模型計(jì)算了χ＝0.8在不同粘性時(shí)對(duì)應(yīng)的滲透率.計(jì)算區(qū)域格點(diǎn)數(shù)為1283，根據(jù)解析解確定外力大小使得流動(dòng)雷諾數(shù)Re固定為0.1.計(jì)算使用雙精度浮點(diǎn)數(shù)，收斂標(biāo)準(zhǔn)為相鄰1 000個(gè)演化步計(jì)算所得滲透率相對(duì)變化小于10－5.使用LBGK和MRT模型的計(jì)算結(jié)果與解析解得比值如圖2所示.從圖中可以看出，在同一孔隙率情況下使用不同松弛時(shí)間（即不同的運(yùn)動(dòng)粘度ν），LBGK模型模擬的結(jié)果與解析解相差明顯，而使用MRT模型則結(jié)果較理想.由此可見(jiàn)，本文模擬結(jié)果驗(yàn)證了Pan等人結(jié)論，驗(yàn)證了GPU程序的正確性.

為考察GPU浮點(diǎn)數(shù)精度對(duì)計(jì)算結(jié)果的影響，測(cè)試了χ＝0.8，＝1.0時(shí)，使用單精度和雙精度計(jì)算時(shí)K/K?的收斂過(guò)程，結(jié)果如圖3所示，單精度計(jì)算演化9 300步后滲透率收斂到0.986 2，而雙精度計(jì)算演化9 400步收斂到0.986 4，單精度相對(duì)于雙精度誤差為0.02%.說(shuō)明在GPU上使用單精度浮點(diǎn)數(shù)計(jì)算滲流問(wèn)題，由于計(jì)算精度引起的誤差較小.下文計(jì)算性能測(cè)試中如無(wú)單獨(dú)說(shuō)明，均使用單精度計(jì)算.

3.2 性能測(cè)試

性能測(cè)試計(jì)算環(huán)境為Intex Xeor X5550 CPU，如無(wú)特殊說(shuō)明，GPU均使用Nvida Tesla C1060.CUDA toolkit版本為3.2，操作系統(tǒng)為 CentOS 5.8 64bit版，程序使用CUDA擴(kuò)展的C語(yǔ)言.計(jì)算速度衡量標(biāo)準(zhǔn)為每秒百萬(wàn)格點(diǎn)更新速率（MFLUPS），其計(jì)算方法為

圖3 使用不同浮點(diǎn)數(shù)精度計(jì)算的收斂過(guò)程Fig.3 Convergence processwith different float precisions

其中Tevol為演化步數(shù)，在本工作中，Tevol均取10 000步，Nf為流體格點(diǎn)個(gè)數(shù)，tcomp為總的計(jì)算時(shí)間.

圖4所示為采用基本算法（完全存貯算法）和稀疏存貯算法在不同孔隙率下的計(jì)算速度對(duì)比，兩種算法均使用LBGK模型，網(wǎng)格規(guī)模為1283.從圖中可以發(fā)現(xiàn)，基本算法計(jì)算速度隨孔隙率的減小近似線性下降，當(dāng)孔隙率從0.95減為0.3左右時(shí)，計(jì)算速度從180 MFLUPS降到了80 MFLUPS以下，與第二節(jié)中的分析相符，即這是由于無(wú)效的固體格點(diǎn)占據(jù)了GPU流處理單元的執(zhí)行時(shí)間.而作為對(duì)比，稀疏存貯算法計(jì)算速度基本不受孔隙率影響，在所考察的孔隙率范圍內(nèi)基本能維持在200MFLUPS以上.另外可以發(fā)現(xiàn)在孔隙率接近1即計(jì)算區(qū)域幾乎不含固體格點(diǎn)時(shí)，稀疏存貯算法計(jì)算速度仍比基本算法高，這是因?yàn)樵谶@種情況下，雖然基本算法的Kernel中判斷語(yǔ)句幾乎沒(méi)有分支，但是其總體指令數(shù)量較稀疏存貯算法大得多，并且需要顯式處理流場(chǎng)外邊界.后文分析影響性能的其它因素均使用稀疏存儲(chǔ)算法.

圖4 稀疏存貯算法與基本算法的計(jì)算速度Fig.4 Performance of different algorithms

值得指出的是，同樣使用稀疏存貯算法的CPU串行程序，用4.6版本的gcc編譯器結(jié)合O2級(jí)優(yōu)化選項(xiàng)編譯時(shí)，在各種孔隙率下計(jì)算速度均低于2MFLUPS.可見(jiàn)使用GPU加速的稀疏存儲(chǔ)算法相對(duì)于CPU平臺(tái)加速比能達(dá)到兩個(gè)量級(jí).

為分析碰撞模型（LBGK/MRT）和浮點(diǎn)數(shù)精度對(duì)計(jì)算速度的影響，測(cè)試了在不同網(wǎng)格規(guī)模下，使用不同碰撞模型和浮點(diǎn)數(shù)精度時(shí)的計(jì)算速度.測(cè)試時(shí)，孔隙率固定為0.32（BCC結(jié)果所能達(dá)到的最小孔隙率），計(jì)算網(wǎng)格規(guī)模從403變化到3203.為使性能規(guī)律更具一般性，在更新一代的Tesla K20 GPU也做了相同測(cè)試.兩種GPU上的測(cè)試結(jié)果如圖5和圖6所示.

從圖5和圖6中均可發(fā)現(xiàn)當(dāng)Nf較小時(shí)，計(jì)算速度隨Nf增大而增大，而當(dāng)Nf增大到一定程度時(shí)，計(jì)算速度基本保持不變.從圖中還可發(fā)現(xiàn)在某些特定的網(wǎng)格規(guī)模時(shí)，計(jì)算速度呈現(xiàn)突然的降低（圖中曲線的V字形區(qū)域）.在使用單精度浮點(diǎn)數(shù)或C1060 GPU計(jì)算時(shí)，這一現(xiàn)象更為明顯.文獻(xiàn)［14］也觀察到了這一現(xiàn)象，盡管該文沒(méi)有采用本文使用的稀疏存貯算法.出現(xiàn)該現(xiàn)象是因?yàn)槊總€(gè)格點(diǎn)讀取周圍格點(diǎn)DF時(shí)不滿足合并訪問(wèn)條件，這會(huì)導(dǎo)致GPU流處理器效率在某些特殊的計(jì)算規(guī)模處出現(xiàn)很大幅度的下降［14］，相比較而言，在較新構(gòu)架的K20 GPU上，由于非合并訪問(wèn)造成的的性能下降較小，所以上述現(xiàn)象相對(duì)不明顯.

對(duì)比圖5和圖6可以發(fā)現(xiàn)，K20 GPU計(jì)算速度約為Tesla C1060的2.5到3倍，這是因?yàn)長(zhǎng)BM算法的性能瓶頸正是存儲(chǔ)帶寬，而K20顯存帶寬為C1060顯存帶寬的2.8倍.

從圖5和圖6中還可以看出，對(duì)于同一種碰撞模型，雙精度（DP）浮點(diǎn)數(shù)計(jì)算比單精度（SP）計(jì)算大約慢50%，這是因?yàn)槭遣捎秒p精度浮點(diǎn)數(shù)計(jì)算時(shí)顯存數(shù)據(jù)傳輸量增大，并且寄存器文件使用增多導(dǎo)致GPU執(zhí)行效率降低.從圖中還可發(fā)現(xiàn)，在K20 GPU上，MRT和LBGK碰撞模型的計(jì)算速度幾乎相同，而在C1060 GPU上，使用雙單精度時(shí)，MRT模型計(jì)算速度明顯比LBGK慢.理論上，MRT模型計(jì)算量只比LBGK模型多15%～20%［1］，而訪存量沒(méi)有變化，按照LBM算法的性能瓶頸在顯存帶寬這一點(diǎn)來(lái)看，兩種碰撞模型的計(jì)算速度應(yīng)該沒(méi)有顯著差異.但由于MRT模型的Kernel函數(shù)中需要的寄存器文件比LBGK多，在C1060 GPU上寄存器容量較為有限，使用雙精度浮點(diǎn)數(shù)時(shí)，寄存器容量限制了其運(yùn)行效率.而在構(gòu)架較新的K20 GPU上，寄存器容量相對(duì)充裕，所以無(wú)論是單精度還是雙精度，LBGK和MRT模型計(jì)算速度幾乎相同.

圖5 Tesla C1060 GPU計(jì)算速度Fig.5 Computing speed with Tesla C1060 GPU

圖6 Tesla K20 GPU計(jì)算速度Fig.6 Computing speed with Tesla K20 GPU

最后，對(duì)于稀疏存儲(chǔ)算法，流體格點(diǎn)存貯順序可能會(huì)影響對(duì)計(jì)算速度.在上文計(jì)算中，提取三維孔隙結(jié)構(gòu)中流體格點(diǎn)構(gòu)造稀疏存貯數(shù)組時(shí)，均是依次按照Z(yǔ)、Y、X的順序遍歷.這里我們測(cè)試了按其他遍歷次序時(shí)的計(jì)算性能，但發(fā)現(xiàn)計(jì)算速度沒(méi)有顯著差別.如果相應(yīng)地重新排布離散速度順序，則遍歷次序可能會(huì)對(duì)計(jì)算速度有明顯影響［15］，這一問(wèn)題還有待進(jìn)一步研究.

4 結(jié)論

將Pan等人［8］提出的稀疏存貯LBM算法運(yùn)用到了GPU加速的多孔介質(zhì)孔隙尺度模擬.性能測(cè)試結(jié)果表明，在GPU上，使用這種算法能大幅提高加速比并節(jié)省顯存，并且孔隙率越低，效果越明顯.測(cè)試了GPU浮點(diǎn)數(shù)精度對(duì)滲透率計(jì)算結(jié)果的影響，結(jié)果表明單雙計(jì)算精度計(jì)算結(jié)果相差極小，滿足工程應(yīng)用要求.分析了在C1060和較新的K20 GPU上，碰撞模型、浮點(diǎn)數(shù)精度以及計(jì)算規(guī)模對(duì)該算法計(jì)算速度的影響，測(cè)試結(jié)果表明，使用雙精度浮點(diǎn)數(shù)要比單精度性能下降約一倍.在較新的K20 GPU上，無(wú)論使用單精度還是雙精度浮點(diǎn)數(shù)，MRT和LBGK碰撞模型性能相同，而在上一代的C1060 GPU上，使用單精度浮點(diǎn)數(shù)時(shí)兩種碰撞模型計(jì)算速度相同，但使用雙精度浮點(diǎn)數(shù)時(shí)，MRT模型明顯比LBGK模型慢.值得指出的是，以上工作都是基于標(biāo)準(zhǔn)均勻格子的并行化方法，在下一步的研究中，我們將考慮在非常規(guī)網(wǎng)格如樹(shù)形網(wǎng)格［16］上結(jié)合GPU實(shí)現(xiàn)并行化.

［1］ Guo Z L，Zheng C G.Theory and application of lattice Boltzmann method［M］.Beijing：Science Press，2009：202－204.

［2］ Deng Y Q，Tang Z，Dong Y H.Lattice Boltzmann method for simulating propagating acoustic waves［J］.Chinese Journal of Computational Physics，2013，30（6）：808－814.

［3］ Kirk D B，Hu W M.Programmingmassively parallel processors［M］.Beijing：Tsinghua University Press，2010：31－32.

［4］ Huang S，Xiao JY，Hu Y C，et al.GPU?accelerated boundary elementmethod for Burton?Miller equation in acoustics［J］. Chinese Journal of Computational Physics，2011，28（4）：481－487.

［5］ Jonas T.Implementation of a lattice Boltzmann kernel using the compute unified device architecture developed by nVIDIA ［J］.Computing and Visualization in Science，2010，13（1）：29－39.

［6］ Huang C S，Zhang W H，Hou ZM，et al.CUDA based lattice Boltzmann method：Algorithm design and program optimization ［J］.Chinese Science Bulletin，2011，56（28）：2434－2444.

［7］ Xiong Q G，Li B，Xu J，et al.Efficient parallel implementation of the lattice Boltzmann method on large clusters of graphic processing units［J］.Chinese Science Bulletin，2012，57（7）：707－715.

［8］ Pan C X，Prins J F，Miller C T.A high?performance lattice Boltzmann implementation to model flow in porous media［J］. Computer Physics Communications，2004，158：89－105.

［9］ He X Y，Shan XW，Doolen G D.Discrete Boltzmann equation model for nonideal gases［J］.Physical Review E，1998，57 （1）：13.

［10］ Pan C X，Luo L S，Miller C T.An evaluation of lattice Boltzmann schemes for porousmedium flow simulaion［J］.Computers ＆Fluids，35（8－9）：898－909.

［11］ Ginzburg I，Carlier JP，Kao C.Lattice Boltzmann approach to Richards'equation［C］∥Proceedings of the XVth International Conference on Computational Methods in Water Resources，Chapel Hill：Elsevier，2004：15－23.

［12］ Obrecht C，Kuznik F，Tourancheau B，etal.A new approach to the lattice Boltzmannmethod for graphics processing units［J］. Computers＆Mathematicswith Applications，2011，61（12）：3628－3638.

［13］ Sangani A S，Acrivos A.Slow flow through a periodic array of spheres［J］.International Journal of Multihase Flow，1982，8 （4）：343－360.

［14］ McClure JE，Prinsy JF，Miller C T.Comparison of CPU and GPU implementations of the lattice Boltzmann method［C］∥XVIIth International Conference on Computational Methods in Water Resources，Barelona：CIMNE，2010：1－8.

［15］ Mattila K，Hyv?luoma J，Timonen，et al.Comparison of implementations of the lattice?Boltzmann method［J］.Computers＆Mathematics with Applications，2007，55（7）：1514－1524.

［16］ Chen Y，Xia Z H，Cai Q D.Lattice Boltzmann method with tree?structured mesh and treatment of curved boundaries［J］. Chinese Journal of Computational Physics，2010，27（1）：23－30.

GPU Accelerated Lattice Boltzmann Simulation of Flow in Porous M edia

ZHU Lianhua，GUO Zhaoli （State Key Laboratory ofCoal Combustion，Huazhong Uniersity of Science and Technology，Wuhan 430074，China）

A sparse lattice representation lattice Boltzmannmethod algorithm is implemented on Graphics Processing Units（GPU）to accelerate pore scale flow simuation.Prefomance testing shows that sparse lattice representation approach grately reduces memory requirement and maintains performance under low porosity compared with basic algorithm.Overall speedup reaches two orders of magnitude compared with serial code.Various factors including collisionmodel，floatnumber precision，and GPU thataffect computing speed of the algorithm are invesgated independently.It indicates that MRT model runs as fast as LBGK model on new generation of GPU cards.While on old GPU cards，MRTmodel's computing speed matchs LBGK only when using single precision float.

porousmedia；GPU；lattice Boltzmannmethod；parallel computing

TQ021.1

A

2013－12－10；

2014－04－02

國(guó)家自然科學(xué)基金（51125024）資助項(xiàng)目

朱煉華（1992－），男，湖北孝感，碩士生，主要從事多相滲流數(shù)值模擬研究，Email：lhzhu＠hust.edu.cn

Received date： 2013－12－10；Revised date： 2014－04－02