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

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

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

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

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