摘 要：针对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.