摘  要： 針對LBM在應(yīng)用過程中計(jì)算規(guī)模過大的問題，提出了一種優(yōu)化算法。該算法明顯改善了計(jì)算時(shí)間過長的問題；對訪問的邊界數(shù)據(jù)進(jìn)行特殊處理，減少了分支判斷邏輯；根據(jù)三維巖石微觀圖像的特點(diǎn)進(jìn)行了空間結(jié)構(gòu)優(yōu)化，避免了無效格子點(diǎn)在線程中的運(yùn)行，提高了GPU運(yùn)算效率。實(shí)驗(yàn)結(jié)果表明，優(yōu)化后的計(jì)算速度有明顯提升。
關(guān)鍵詞： Boltzmann；并行計(jì)算；GPU；CUDA

    研究巖石的滲透特性、孔隙特性以及變化規(guī)律是預(yù)測油氣形成、分布規(guī)律與探索提高采收率的主要方法[1]。滲透特性反映流體通過巖石空間的能力，以往大多通過實(shí)驗(yàn)來定性研究，但近年來，為從微觀尺度對滲流特性進(jìn)行定量分析，國內(nèi)外學(xué)者在該領(lǐng)域展開了大量研究。為了使研究的滲流結(jié)果具有實(shí)際意義，建立一個(gè)能夠真實(shí)反映巖心孔隙特征的數(shù)字巖心往往是開展研究的第一步。目前，國內(nèi)外學(xué)者能夠借助數(shù)值重建方法或物理實(shí)驗(yàn)手段建立數(shù)字巖心[2]。利用數(shù)字巖心進(jìn)行滲透率分析的方法主要有兩種：孔隙網(wǎng)絡(luò)模型和格子Boltzman方法LBM（Lattice Boltzmann Method）?？紫毒W(wǎng)絡(luò)模型是以孔隙和喉道為基本組成單元的模型，它對不規(guī)則的孔隙喉道進(jìn)行規(guī)則化抽象處理，將其等效為具有任意橫截面的管道，并且采用形狀因子記錄真實(shí)的形狀信息，使幾何計(jì)算得以簡化，從而使得更多的物理公式能直接用于滲流模擬分析?？紫毒W(wǎng)絡(luò)模型是根據(jù)研究的問題定義具體的規(guī)則進(jìn)行流動(dòng)模擬的[3]。LBM把模擬介質(zhì)看作由細(xì)微的“格子”狀的基本單元組成，通過模擬微觀粒子的運(yùn)動(dòng)過程和分布規(guī)律來模擬滲流過程。LBM是一種成熟的算法，被廣泛應(yīng)用于多種領(lǐng)域的流體運(yùn)動(dòng)模擬[4-5]。但是傳統(tǒng)算法的計(jì)算規(guī)模過大，運(yùn)算時(shí)間太長，在實(shí)際應(yīng)用中極為耗時(shí)。為了解決模擬時(shí)間過長的問題，國內(nèi)外學(xué)者都在研究一些基于并行技術(shù)的加速方案，這些方案多采用集群并行計(jì)算機(jī)來加速LBM模擬的方法。
　本文主要以D3Q19為模型來模擬流體，并針對LBM的算法特點(diǎn)對訪問的邊界數(shù)據(jù)進(jìn)行特殊處理，減少了訪問邊界數(shù)據(jù)時(shí)的分支跳轉(zhuǎn)；同時(shí)，根據(jù)三維微觀巖石的空間結(jié)構(gòu)特點(diǎn)，剔除了無效點(diǎn)，減少了運(yùn)算時(shí)間。在國家自然科學(xué)基金項(xiàng)目中，參考文獻(xiàn)[5]在計(jì)算結(jié)構(gòu)和粒子空間上進(jìn)行了優(yōu)化，本文在此基礎(chǔ)上進(jìn)一步優(yōu)化，取得了一定的效果。
1 LBM和CUDA平臺
1.1 LBM的基本方法
　在LBM中，復(fù)雜流體運(yùn)動(dòng)被轉(zhuǎn)化為規(guī)則、離散的格子空間中的運(yùn)動(dòng)，一般采用DdQm模型來表示粒子在離散空間中的運(yùn)動(dòng)模型。其中，d代表空間的維數(shù)，m代表離散速度的方向數(shù)。與Navier-Stockes方程求解非線性偏微分方程不同[6]，LBM僅用離散格子的遷移與碰撞來表述流體的運(yùn)動(dòng)。
模擬三維空間中流體的運(yùn)動(dòng)，一般采用D3Q13、D3Q15、D3Q19或D3Q27模型。在DdQm模型中，m越大時(shí)，LBM的精度越高，收斂速度越慢。
1.2 CUDA平臺
　在20世紀(jì)90年代，NVIDIA公司提出了圖形處理器GPU（Graphics Processing Unit），并且實(shí)現(xiàn)了硬件加速的三維坐標(biāo)變換和光源計(jì)算T&L（Transformation and Lighting）技術(shù)。在數(shù)字圖像應(yīng)用（游戲、工程制圖等）的強(qiáng)烈刺激下，GPU的性能不斷進(jìn)步。在2001年，NVIDIA推出了一款支持可編程圖形流水線（Programmable Graphics Pipeline）的標(biāo)志性產(chǎn)品GeForce 3。隨后，人們開始注意到GPU的高度并行化架構(gòu)的優(yōu)勢，通用計(jì)算圖形處理器GPGPU技術(shù)逐步發(fā)展起來，并且在科學(xué)計(jì)算領(lǐng)域得到了應(yīng)用[7]。在2006年，NVIDIA公司發(fā)布了GeForce 8系列GPU，同時(shí)發(fā)布了GPU運(yùn)用于通用計(jì)算的統(tǒng)一計(jì)算設(shè)備架構(gòu)CUDA（Compute Unified Device Architecture）[8]。從此，GPU技術(shù)的發(fā)展進(jìn)入了全新的時(shí)代。目前，GPU通用計(jì)算技術(shù)已廣泛應(yīng)用于流體滲流模擬、天文計(jì)算、實(shí)時(shí)數(shù)字圖像處理與醫(yī)療診斷等眾多科學(xué)計(jì)算領(lǐng)域。在全球超級計(jì)算機(jī)前10排行榜中，有3個(gè)采用了CPU加GPU的異構(gòu)方式[9]。
2 CUDA平臺上LBM并行算法優(yōu)化
　本文主要從兩個(gè)方面對CUDA平臺上LBM并行算法進(jìn)行優(yōu)化：一是針對邊界數(shù)據(jù)進(jìn)行處理，減少訪問存儲器次數(shù)；二是針對巖石的多孔介質(zhì)的空間特點(diǎn)，從實(shí)質(zhì)上提高運(yùn)算過程中Wrap內(nèi)有效點(diǎn)的比例，從而達(dá)到提高運(yùn)算速度的目的。首先，在CUDA平臺上對LBM進(jìn)行層次化劃分。CUDA線程劃分方式如圖1所示。
本文選擇每個(gè)CUDA線程處理一個(gè)LBM的離散格子。GPU中采用網(wǎng)格分層的方式來組織線程，合理地劃分Grid與Block才能充分發(fā)揮GPU的計(jì)算能力。
2.1 邊界數(shù)據(jù)處理
　在圖1所示的CUDA線程劃分方式下，分析粒子空間中格子數(shù)據(jù)的具體處理過程?？紤]x與y方向的邊界格子中粒子的碰撞與遷移過程。因?yàn)樵诹Ｗ友刈鴺?biāo)軸正方向遷移時(shí)，粒子分布函數(shù)將會(huì)離開格子空間發(fā)生越界，所以在這里對邊界數(shù)據(jù)進(jìn)行特殊處理。

2.2 空間結(jié)構(gòu)優(yōu)化
　對于三維巖石微觀圖像而言，其孔隙空間具有多孔介質(zhì)形態(tài)的特點(diǎn)，并且大多數(shù)情況下其孔隙度不超過20%。在滲流分析中，LBM的離散格子空間中有超過80%的粒子點(diǎn)不參與計(jì)算，如果這些格子點(diǎn)在參與計(jì)算，那么GPU中的Wrap有效率會(huì)大大降低。針對這種情況，可以在進(jìn)行GPU執(zhí)行語句前，在CPU程序中通過分支語句跳過。然而，引入分支條件一般會(huì)降低程序的運(yùn)行效率（Wrap內(nèi)的分支條件會(huì)被串行執(zhí)行）。
在GPU中，計(jì)算與訪存均以Wrap（即32個(gè)線程）為單位。經(jīng)過統(tǒng)計(jì)，在格子空間規(guī)模為128×128×128，孔隙度為12.42%的巖石孔隙結(jié)構(gòu)中，有效點(diǎn)數(shù)為0的Wrap占總Wrap數(shù)的30%左右，有效點(diǎn)數(shù)越多的Wrap所占的百分比越小，如圖3所示。有效點(diǎn)不為0或者32的Wrap內(nèi)會(huì)出現(xiàn)分支條件，且有效點(diǎn)數(shù)為0的Wrap是無效的，其內(nèi)部所有線程均不作任何有效計(jì)算。并且，GPU平臺最優(yōu)的是共享存儲器，而不是全局存儲器，故訪問全局存儲器相對耗時(shí)。如果能避免這些Wrap的無效點(diǎn)，那么訪問全局存儲器的數(shù)目就會(huì)大大降低，從而提高運(yùn)算速度。為此，本文引入了分支條件，對于無效格子點(diǎn)，既不讀取任何數(shù)據(jù)也不執(zhí)行任何計(jì)算。該分支條件對有效點(diǎn)數(shù)為32的Wrap不會(huì)產(chǎn)生任何影響。

　圖5中數(shù)據(jù)表明，經(jīng)過3種顯卡測試，隨著孔隙度的減小，優(yōu)化后均得到了明顯的加速效果?？紫抖仍叫?，Wrap中無效點(diǎn)的數(shù)目就越多，經(jīng)過判斷跳過這些大量的無效點(diǎn)，運(yùn)算時(shí)間大量減少。因?yàn)榭臻g結(jié)構(gòu)優(yōu)化會(huì)受孔隙度的影響，所以不同孔隙度下的運(yùn)行結(jié)果會(huì)有較大差距。
3.2 優(yōu)化結(jié)果
　實(shí)驗(yàn)采用孔隙度為0.26的一組微觀巖石圖像，在9600GT、GT240和GTX560Ti 3種不同性能顯卡上進(jìn)行實(shí)驗(yàn)。圖6為CPU在不同圖片規(guī)模下的運(yùn)算結(jié)果。圖7為加速比，可以看出，經(jīng)過優(yōu)化后的計(jì)算速度有了顯著提升。
　從圖7可以看出，隨著圖像規(guī)模的增大，加速比有明顯的提升。因?yàn)殡S著圖像規(guī)模增大，在格子遷移到邊界的數(shù)目增多，經(jīng)過邊界數(shù)據(jù)處理后粒子直接命中下一個(gè)，分支判斷大量減少，提高了GPU中線程的使用效率，提高了運(yùn)算速度。在圖片規(guī)模為1283、顯卡為GTX560Ti的情況下，優(yōu)化加速比可以達(dá)到12左右。與參考文獻(xiàn)[5]的實(shí)驗(yàn)結(jié)果對比，本文的實(shí)驗(yàn)結(jié)果加速比提高了1.9倍以上。

　本文旨在提高Boltzmann算法在三維巖石微觀圖像上的計(jì)算速度，減少計(jì)算時(shí)間。對邊界數(shù)據(jù)進(jìn)行處理，在邊界遷移時(shí)減少了分支判斷條件，減少了訪問全局存儲器的數(shù)目；針對三維巖石微觀圖像的空間結(jié)構(gòu)特點(diǎn)，減少了Wrap內(nèi)無效點(diǎn)數(shù)目，提升了GPU的線程利用率。三維巖石微觀圖像的流體滲透率的計(jì)算規(guī)模過大，計(jì)算時(shí)間較長的問題得到一定的改善。這擴(kuò)大了Boltzmann算法在計(jì)算流體滲透率中的應(yīng)用，使Boltzmann算法在較大規(guī)模的巖石微觀圖像的應(yīng)用上更加具有實(shí)際意義。預(yù)計(jì)將來GPU計(jì)算將會(huì)代替CPU計(jì)算成為巖石微觀圖像滲流分析的主要手段。
參考文獻(xiàn)
[1] 陳杰，周改英，趙喜亮，等.儲層巖石孔隙結(jié)構(gòu)特征研究方法綜述[J].特種油氣藏，2005，12（4）：11-14.
[2] 趙秀才.數(shù)字巖心及孔隙網(wǎng)絡(luò)模型重構(gòu)方法研究[D].北京：石油工程學(xué)院，2009.
[3] MANWART C， AALTOSALMI U， KOPONEN A， et al. Lattice-Boltzmann and finite-difference simulations for the permeability for three-dimensional porous media[J]. Physical Review E， 2002，66（1）：16702.
[4] GLADROW W. Lattice-gas cellular automata and Lattice Boltzmann models[D]. Springer， 2005.
[5] Zeng Fansong， Teng Qizhi. Application and optimization for Lattice Boltzmann method in calculation of fluid permeability[J]. Computer Engineering and Applications， 2011，47（22）： 156-159.
[6] Ming Pingjian. Numerical simulation of sloshing in rectangular tank with vof based on unstructured grids[J]. Journal of Hydrodynamics， 2010， 22（6）： 856-864.
[7] 仇得元.GPU編程技術(shù)——從GLSL、CUDA到OpenGL[M].北京：機(jī)械工業(yè)出版社，2011.
[8] 張舒，褚艷利，趙開勇，等.GPGPU高性能運(yùn)算之CUDA[M].北京：中國水利水電出版社，2009.
[9] SRIVASTAVA N. Top 10 super computer in the world [EB/OL].http：//www.top10gallery.com/2012/02/top-10-super-computer-in-world-2012.html， 2012-2.
[10] NVIDIA CUDA C Best Practices Guide[Z]. Version 4.0. NVIDIA Corporation， 2011：103-105.
[11] VALDERHAUG T K. The Lattice Boltzmann simulation on multi-GPU systems[D]. Trondheim： Norwegian University of Science and Technology， 2011.
[12] HABICH J. Performance evaluation of numeric compute kernels on NVIDIA GPUs[M]. Master′s Thesis of Friedrich-Alexander-Universit Erlangen-Nurnberg， 2008.
[13] 李博，李曦鵬，張?jiān)疲?耦合Nvidia/AMD兩類GPU的格子玻爾茲曼模擬[J].科學(xué)通報(bào)，2009，54（20）：3177-3184.
[14] 李森，李新亮，王龍，等.基于OpenGL的并行方腔流加速性能分析[J].計(jì)算機(jī)應(yīng)用，2011，28（4）：1401-1403，1421.