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

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

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

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

　本文旨在提高Boltzmann算法在三維巖石微觀圖像上的計算速度,，減少計算時間。對邊界數(shù)據(jù)進行處理,，在邊界遷移時減少了分支判斷條件,，減少了訪問全局存儲器的數(shù)目,；針對三維巖石微觀圖像的空間結(jié)構(gòu)特點，減少了Wrap內(nèi)無效點數(shù)目,，提升了GPU的線程利用率,。三維巖石微觀圖像的流體滲透率的計算規(guī)模過大，計算時間較長的問題得到一定的改善,。這擴大了Boltzmann算法在計算流體滲透率中的應用,，使Boltzmann算法在較大規(guī)模的巖石微觀圖像的應用上更加具有實際意義。預計將來GPU計算將會代替CPU計算成為巖石微觀圖像滲流分析的主要手段,。
參考文獻
[1] 陳杰,，周改英，趙喜亮,，等.儲層巖石孔隙結(jié)構(gòu)特征研究方法綜述[J].特種油氣藏,，2005，12（4）：11-14.
[2] 趙秀才.數(shù)字巖心及孔隙網(wǎng)絡模型重構(gòu)方法研究[D].北京：石油工程學院,，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編程技術——從GLSL,、CUDA到OpenGL[M].北京：機械工業(yè)出版社,，2011.
[8] 張舒，褚艷利,，趙開勇,，等.GPGPU高性能運算之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] 李博,，李曦鵬，張云,，等.耦合Nvidia/AMD兩類GPU的格子玻爾茲曼模擬[J].科學通報,，2009,，54（20）：3177-3184.
[14] 李森，李新亮,，王龍,，等.基于OpenGL的并行方腔流加速性能分析[J].計算機應用，2011,，28（4）：1401-1403,，1421.