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

    研究巖石的滲透特性,、孔隙特性以及變化規(guī)律是預測油氣形成、分布規(guī)律與探索提高采收率的主要方法[1],。滲透特性反映流體通過巖石空間的能力,，以往大多通過實驗來定性研究，但近年來,，為從微觀尺度對滲流特性進行定量分析，國內(nèi)外學者在該領域展開了大量研究,。為了使研究的滲流結果具有實際意義,，建立一個能夠真實反映巖心孔隙特征的數(shù)字巖心往往是開展研究的第一步。目前,，國內(nèi)外學者能夠借助數(shù)值重建方法或物理實驗手段建立數(shù)字巖心[2],。利用數(shù)字巖心進行滲透率分析的方法主要有兩種：孔隙網(wǎng)絡模型和格子Boltzman方法LBM（Lattice Boltzmann Method）,。孔隙網(wǎng)絡模型是以孔隙和喉道為基本組成單元的模型,，它對不規(guī)則的孔隙喉道進行規(guī)則化抽象處理,，將其等效為具有任意橫截面的管道，并且采用形狀因子記錄真實的形狀信息,，使幾何計算得以簡化,，從而使得更多的物理公式能直接用于滲流模擬分析?？紫毒W(wǎng)絡模型是根據(jù)研究的問題定義具體的規(guī)則進行流動模擬的[3],。LBM把模擬介質看作由細微的“格子”狀的基本單元組成，通過模擬微觀粒子的運動過程和分布規(guī)律來模擬滲流過程,。LBM是一種成熟的算法,，被廣泛應用于多種領域的流體運動模擬[4-5]。但是傳統(tǒng)算法的計算規(guī)模過大,，運算時間太長,，在實際應用中極為耗時。為了解決模擬時間過長的問題,，國內(nèi)外學者都在研究一些基于并行技術的加速方案,，這些方案多采用集群并行計算機來加速LBM模擬的方法。
　本文主要以D3Q19為模型來模擬流體,，并針對LBM的算法特點對訪問的邊界數(shù)據(jù)進行特殊處理,，減少了訪問邊界數(shù)據(jù)時的分支跳轉；同時,，根據(jù)三維微觀巖石的空間結構特點,，剔除了無效點，減少了運算時間,。在國家自然科學基金項目中,，參考文獻[5]在計算結構和粒子空間上進行了優(yōu)化，本文在此基礎上進一步優(yōu)化,，取得了一定的效果,。
1 LBM和CUDA平臺
1.1 LBM的基本方法
　在LBM中，復雜流體運動被轉化為規(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的高度并行化架構的優(yōu)勢,，通用計算圖形處理器GPGPU技術逐步發(fā)展起來,，并且在科學計算領域得到了應用[7]。在2006年,，NVIDIA公司發(fā)布了GeForce 8系列GPU,，同時發(fā)布了GPU運用于通用計算的統(tǒng)一計算設備架構CUDA（Compute Unified Device Architecture）[8]。從此,，GPU技術的發(fā)展進入了全新的時代,。目前，GPU通用計算技術已廣泛應用于流體滲流模擬,、天文計算,、實時數(shù)字圖像處理與醫(yī)療診斷等眾多科學計算領域。在全球超級計算機前10排行榜中,，有3個采用了CPU加GPU的異構方式[9],。
2 CUDA平臺上LBM并行算法優(yōu)化
　本文主要從兩個方面對CUDA平臺上LBM并行算法進行優(yōu)化：一是針對邊界數(shù)據(jù)進行處理，減少訪問存儲器次數(shù),；二是針對巖石的多孔介質的空間特點,，從實質上提高運算過程中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 空間結構優(yōu)化
　對于三維巖石微觀圖像而言，其孔隙空間具有多孔介質形態(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%的巖石孔隙結構中,，有效點數(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)化后均得到了明顯的加速效果,?？紫抖仍叫。琖rap中無效點的數(shù)目就越多,，經(jīng)過判斷跳過這些大量的無效點,，運算時間大量減少。因為空間結構優(yōu)化會受孔隙度的影響,，所以不同孔隙度下的運行結果會有較大差距,。
3.2 優(yōu)化結果
　實驗采用孔隙度為0.26的一組微觀巖石圖像，在9600GT,、GT240和GTX560Ti 3種不同性能顯卡上進行實驗,。圖6為CPU在不同圖片規(guī)模下的運算結果。圖7為加速比,，可以看出,，經(jīng)過優(yōu)化后的計算速度有了顯著提升。
　從圖7可以看出,，隨著圖像規(guī)模的增大,，加速比有明顯的提升。因為隨著圖像規(guī)模增大,，在格子遷移到邊界的數(shù)目增多,，經(jīng)過邊界數(shù)據(jù)處理后粒子直接命中下一個，分支判斷大量減少,，提高了GPU中線程的使用效率,，提高了運算速度。在圖片規(guī)模為1283,、顯卡為GTX560Ti的情況下,，優(yōu)化加速比可以達到12左右。與參考文獻[5]的實驗結果對比,，本文的實驗結果加速比提高了1.9倍以上,。

　本文旨在提高Boltzmann算法在三維巖石微觀圖像上的計算速度，減少計算時間。對邊界數(shù)據(jù)進行處理,，在邊界遷移時減少了分支判斷條件,，減少了訪問全局存儲器的數(shù)目；針對三維巖石微觀圖像的空間結構特點,，減少了Wrap內(nèi)無效點數(shù)目,，提升了GPU的線程利用率。三維巖石微觀圖像的流體滲透率的計算規(guī)模過大,，計算時間較長的問題得到一定的改善。這擴大了Boltzmann算法在計算流體滲透率中的應用,，使Boltzmann算法在較大規(guī)模的巖石微觀圖像的應用上更加具有實際意義,。預計將來GPU計算將會代替CPU計算成為巖石微觀圖像滲流分析的主要手段。
參考文獻
[1] 陳杰,，周改英,，趙喜亮，等.儲層巖石孔隙結構特征研究方法綜述[J].特種油氣藏,，2005,，12（4）：11-14.
[2] 趙秀才.數(shù)字巖心及孔隙網(wǎng)絡模型重構方法研究[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.