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

朱煉華，郭照立

（華中科技大學(xué)煤燃燒國家重點實驗室，武漢 430074）

文章編號：1001?246X（2015）01?0020?07

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

朱煉華，郭照立

（華中科技大學(xué)煤燃燒國家重點實驗室，武漢 430074）

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

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

0 引言

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

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

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

1 計算模型

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

其中下標i表示離散速度方向，x為空間位置，t為時間，δt為時間步長，ci為格子離散速度，其各分量為（2）

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

其中cs＝1/3為格子聲速， 是無量綱松弛時間，與運動學(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ù)得到

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

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

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

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

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

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

2.1 CUDA構(gòu)架

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

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

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

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

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

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

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

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

2.3 稀疏存儲算法

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

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

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

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

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

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

3 數(shù)值算例

3.1 驗證算例

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

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

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

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

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

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

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

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

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

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

3.2 性能測試

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

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

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

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

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

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

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

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

對比圖5和圖6可以發(fā)現(xiàn)，K20 GPU計算速度約為Tesla C1060的2.5到3倍，這是因為LBM算法的性能瓶頸正是存儲帶寬，而K20顯存帶寬為C1060顯存帶寬的2.8倍.

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

圖5 Tesla C1060 GPU計算速度Fig.5 Computing speed with Tesla C1060 GPU

圖6 Tesla K20 GPU計算速度Fig.6 Computing speed with Tesla K20 GPU

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

4 結(jié)論

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

國家自然科學(xué)基金（51125024）資助項目

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

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