CN102681972A - 一种利用GPU加速格子-Boltzmann的方法 - Google Patents

Abstract

本发明提供了一种利用GPU加速格子-Boltzmann的方法，涉及主机端的CPU和设备端的GPU。该方法包括：主机端根据物理问题，给定计算域、参考长度、来流速度、密度和雷诺数等参数，并根据网格划分设计内核的线程数；设备端通过这些宏观参量(密度、速度、雷诺数和黏性系数等)计算出所有格点上各个方向的平衡态分布函数，以此作为计算的初场，并行求解离散方程和处理边界，并把最终迭代得到的结果回传到主机端。本发明利用GPU设备端快速计算的特点参与对格子-Boltzmann中的迁移和碰撞的计算，通过CPU主机端和GPU设备端的协调运作加速格子-Boltzmann的迭代过程。

Description

一种利用GPU加速格子-Boltzmann的方法

技术领域

本发明涉及计算机高性能计算领域、计算流体力学领域，具体涉及一种利用NVIDIA的GPU加速格子-Boltzmann的方法。

背景技术

格子Boltzmann方法（Lattice Boltzmann Method，LBM）在过去的20年里已经发展成为一种有效的数值模拟方法，它是介于微观分子动力学方法和基于连续介质假设的宏观方法之间的一种介观方法。该方法与传统的流体模拟方法不同，它基于分子动理论，通过跟踪粒子分布函数的输运而后对分布函数求矩来获得宏观平均特性。格子 Boltzmann 方法的动理论特性使得它在模拟许多复杂流动的时候更有效， 如多孔介质流动、悬浮流、多相流、多组分流等。格子Boltzmann 方法具有天生的并行特性，以及边界处理简单、程序易于实施等优点。

采用LBM方法求解物理问题时的基本过程如图1所示。对于一个特定的物理问题，首先进行辅助步骤：

1)基于各种简化假设，进行物理建模，确定出计算区域、初始条件以及边界条件等，并根据物理问题的不同，选择相应的格子Boltzmann模型；

2)进行网格划分，假设网格划分大小为NX*NY；

3)根据不同格子Boltzmann模型，选择控制方程，并对其进行离散。如采用标准格子 Boltzmann 方法对等温不可压缩流动进行模拟，则离散后的控制方程为 LBGK 方程；

前面这3步是在数值模拟之前就进行的。随后进入数值模拟阶段：

4)根据物理问题，给定所有格点上的宏观参量(密度、速度、黏性系数等)，并由此计算出所有格点上各个方向的平衡态分布函数，以此作为计算的初场；

5)求解离散后的控制方程，例如，采用迁移碰撞规则求解 LBGK 方程；

6)根据边界条件，在相应边界格点上实施边界处理格式；

7)基于不同格子 Boltzmann 模型的宏观量的定义法则，计算各格点上的宏观参量；

8)判断计算是否收敛；

9)若计算收敛，则输出计算结果；否则返回第 4 步，继续求解，直到收敛为止。

广泛应用的格子-Boltzmann单松弛时间近似BGK模型基于如下的演化方程：                                                

这里， 为粒子分布函数，代表时间t，在处存在以微观速度运动的粒子的概率。松弛时间

代表达到局部平衡的速率，与流体的动粘性系数相关。平衡分布函数

是Maxwell-Boltzmann方程的低马赫数近似，取决于流体的密度和流动速度。它们之间的关系由下式决定

    其中，在D2Q9模型中：

 

    流体密度和速度则可由粒子分布函数依下式算出：

离散速度矢量

，和粒子分布函数的个数N取决于所选定的格子-Boltzmann模型，在D2Q9模型中有9个分量，对应的粒子分布函数的个数也为9个，见图2。

自从2006年英伟达（NVIDIA）公司推出图形处理器G80（包含了128个流式多处理器，最新的Fermi架构包含了512个多处理器）以来，图形处理器（GPU，Graphic Processing Unit）在某些大规模并行计算的应用上，相对于CPU来说性能提高可达100倍以上。GPU拥有更多的晶体管，用于数据处理而不是像CPU那样去处理数据cache和指令控制，这意味着GPU具有巨大的并行计算能力。

CUDA C作为GPU的并行编程语言。CUDA C编程将CPU称之为主机，将GPU作为一个协处理器称为设备。在CUDA编程中，多个线程同时执行在一个GPU上，并由多个线程组成一个线程块（Block），多个线程块又组织成网格（Grid）；另外，每32个线程组成一个束（warp）。CUDA编程中常用到的优化技术有合理的网格配置，每个流多处理器（Stream Multiprocessor，SM）上有足够多的warp可隐藏访问延迟，并进行全局存储器的合并访问，共享存储器、纹理存储器和常量存储器的使用，寄存器的合理使用等等。

LBM方法需要大量的计算，以方腔的计算为例，假设网格大小为1024*1024，迭代10000次，每个网格点都要在一次迭代中做一次迁移和碰撞计算，在主频为3.00GHz的QuadCore Intel Xeon X5450进行计算需要数小时的时间，而我们计算更大的网格迭代更多次数时需要数天的时间，这严重影响了LBM方法的性能。目前，人们往往使用大规模的X86服务器集群来进行LBM处理，其原理是先将计算负载进行划分然后分配到各个格点，然后由各个格点单独进行计算并在每个迭代步后交换数据，最后将结果汇总输出。这种做法由于CPU峰值浮点计算能力较低而且网络传输开销极大而消耗了大量的时间、电力和维护费用。而且，随着人们对流体模拟的周期要求越来越短，精度要求越来越高，PC服务器集群的规模越做越大，在系统构建成本、数据中心机房空间、功耗散热和电力限制、可管理性、编程简易性、扩展性、管理维护费用等方面都面临着巨大的挑战。

可见，为满足流体模拟的需求，需要一种提高LBM的计算性能的方法，并且此方法能降低机房构建成本和管理、运行、维护费用。而GPU可以很好的解决这些问题。

发明内容

本发明的目的在于加速格子Boltzmann方法，提高其处理性能，使CPU与GPU进行协同计算，从而满足流体模拟的需求，并且降低机房构建成本和管理、运行、维护费用。提供一种利用NVIDIA GPU加速格子-Boltzmann的方法。

本发明的目的是按以下方式实现的，将基本参数初始化计算放在CPU端执行，而把耗时并且并行性非常好的平衡分布函数计算、宏观量统计、离散方程求解和边界处理部分利用CUDA技术进行并行化改造，使其在GPU端并行执行， CPU与GPU进行协同计算，最终实现加速格子Boltzmann方法，该方法可以分为以下步骤：

定位格子Boltzmann方法的性能瓶颈；

利用LBM进行流体模拟时，最耗时的计算部分为求解离散方程和边界处理的过程，该过程占据了整个模拟的绝大部分时间，其它部分几乎不耗时，因此，求解离散方程和边界处理的迭代过程是LBM中的性能瓶颈。

1)并行性分析；

根据对LBM算法中求解离散方程和边界处理的串行算法的分析，每个网格点的迁移、碰撞、宏观量统计、平衡态分布函数计算和边界处理的计算都是数据并行的；

2)求解离散方程和边界处理的CUDA并行化改造；

a)求解离散方程可以采用迁移碰撞的过程，宏观量统计、平衡态分布函数计算和碰撞过程中对每个网格的计算之间没有任何依赖性，因此，可以让CUDA中的每个线程负责一个网格划分中的格点的计算，网格划分中的总格点数即为总的线程数；分布函数的迁移只涉及到该格点周围的其他格点，也可以通过单个线程对全局存储器中相关分布函数的读操作来实现；

b)在LBM算法中，对边界要做特殊的处理（非平衡外推，反弹），对于边界上的每个格点之间的计算也没有数据的依赖性，因此，也可以让每个线程负责一个格点的计算；

c)CUDA的线程模型设计：根据GPU存储器的访问特点、GPU的寄存器资源及并行算法来设计CUDA的线程模型，确定总线程数和总线程块数；

d)求解离散方程和边界处理的GPU内核代码编写。

本发明的有益效果是：不但提高了格子Boltzmann方法的处理性能，满足了流体模拟的需求，而且能降低功耗，减少机房构建成本和管理、运行、维护费用，并且这种方法实现简单，需要的开发成本低。

附图说明

图1是 LBM方法模拟求解的基本流程图；

图2是利用GPU加速LBM方法实施例的流程图；

图3 是D2Q9模型；

图4 是迁移过程示意图。

具体实施方式

参照说明书附图对本发明的作以下详细地说明：

为了使本发明的目的、技术方案和优点更加清晰，下面结合附图和实施例，对本发明作以下详细说明。

本发明的目的在于加速格子Boltzmann方法，提高其处理性能，使CPU与GPU进行协同计算，从而满足流体模拟的需求，并且降低机房构建成本和管理、运行、维护费用。本发明中，将需要把初始化计算放在CPU端执行，而把耗时并且并行性非常好的求解离散方程和边界处理部分利用CUDA技术进行并行化改造，使其在GPU端并行执行，CPU与GPU进行协同计算，最终实现加速格子Boltzmann方法，如图3所示。具体步骤及实施过程如下：

1)根据物理问题，在主机端给定计算域上的宏观参量(密度、速度、黏性系数等) ，传递给设备端；

2)定义设备端的数据结构和存储方式，用于存储每个格点各个方向的平衡分布函数和每个格点的速度、密度等宏观参量，由主机端传来的宏观参量计算出所有格点上各个方向的平衡态分布函数，以此作为计算的初场；

3)设计迁移碰撞内核，设计每个block的thread数目为BLOCKSIZE（取值在64-512），线程结构为：Block(BLOCKSIZE, 1)，Grid((NX+BLOCKSIZE-1)/BLOCKSIZE, NY)，并且让内核中的每个线程计算一个格点的迁移和碰撞过程，如图4所示，内核伪代码如下；

1:  k  = gridDim.y * blockIdx.y + blockIdx.x*blockDim.x + threadIdx.x;//k表示网格格点的下标

2:  /*迁移过程，对当前格子点周围相关格子点的分布函数进行汇总式的读操作*/

3:  fr = fr0[k];//0代表上一时层的分布函数

4:  fe = fe0[k-1];

5:  fn = fn0[k-NX];

6:  fw = fw0[k+1];

7:  fs = fs0[k+NX];

8:  fne = fne0[k-NX-1];

9:  fnw = fnw0[k-NX+1];

10: fsw = fsw0[k+NX+1];

11: fse = fse0[k+NX-1];

12: /*碰撞过程*/

13：根据迁移后的分布函数fr-fse 求宏观量

14: 根据宏观量求各个方向的平衡分布函数f1,f2,f3,f4,f5,f6,f7,f8；

15: 根据f1,f2,f3,f4,f5,f6,f7,f8以及迁移后的分布函数fr, fe, fn, fw, fs, fne, fnw, fsw, fsw, fse求碰撞后的分布函数fr1[k]，fe1[k]，fn1[k]，fw1[k]，fs1[k]，fne1[k]，fnw1[k]，fsw1[k]，fse1[k]；

4)在设备端对边界进行处理，边界处理可以采用反弹法、非平衡外推法等方法，在对边界的处理时同样设计每个线程处理一个节点的计算；

5)判断是否迭代完成，完成则输出，否则继续迭代；

6)设备端根据分布函数并行求得速度、密度和流函数等宏观参量并把结果传递给主机端；主机端进行结果输出；

7)性能测试

a)测试环境及测试数据

测试环境包括硬件环境、软件环境、运行软件，其中运行软件包括运行在CPU上的CPU版LBM算法和运行在GPU上的LBM算法；测试数据选取了顶盖驱动方腔流，输入包括网格大小和其他一些输入参数，具体测试环境和测试数据各项参数如下表所示；

b)性能结果

为了保证测试性能结果的稳定性，我们对上述作业进行了10次测试，数据类型为double，CPU版LBM算法在单CPU上运行10次的平均时间是19763秒，而GPU版LBM算法在单GPU上运行上述同样作业10次的平均时间是598秒，GPU版运行的性能是CPU版的19763/598 = 33倍。

由本发明的技术方案可见，本发明通过测试迁移和碰撞以及边界处理部分为LBM算法中的性能瓶颈，并且这部分的数据完全具有独立性，完全适合于GPU上采用CUDA进行并行计算，而把不耗时的初始化参数和结果输出仍然放在CPU端执行，CPU与GPU进行协同计算。通过测试整体性能提高33倍，现在的一个计算GPU计算节点相当于原来的33个以上CPU计算节点集群的计算性能，这样不但满足了流体模拟的需求，而且大大降低了功耗，减少了机房构建成本和管理、运行、维护费用，而且这种方法实现简单，需要的开发成本低。

除说明书所述的技术特征外，均为本专业技术人员的已知技术。

Claims (1)

1.一种利用GPU加速格子-Boltzmann的方法，其特征在于，将基本参数初始化计算放在CPU端执行，而把平衡分布函数计算、宏观量统计、离散方程求解和边界处理部分利用CUDA技术进行并行化改造，使其在GPU端并行执行， CPU与GPU进行协同计算，最终实现加速格子Boltzmann方法，包括以下步骤：

1)定位格子Boltzmann方法的性能瓶颈；

利用LBM进行流体模拟时，最耗时的计算部分为求解离散方程和边界处理的过程，该过程占据了整个模拟的绝大部分时间，其它部分几乎不耗时，因此，求解离散方程和边界处理的迭代过程是LBM中的性能瓶颈；

2)并行性分析；

根据对LBM算法中求解离散方程和边界处理的串行算法的分析，每个网格点的迁移、碰撞、宏观量统计、平衡态分布函数计算和边界处理的计算都是数据并行的；

3)求解离散方程和边界处理的CUDA并行化改造；

a)求解离散方程可以采用迁移碰撞的过程，宏观量统计、平衡态分布函数计算和碰撞过程中对每个网格的计算之间没有任何依赖性，因此，可以让CUDA中的每个线程负责一个网格划分中的格点的计算，网格划分中的总格点数即为总的线程数；分布函数的迁移只涉及到该格点周围的其他格点，也可以通过单个线程对全局存储器中相关分布函数的读操作来实现；

b)在LBM算法中，对边界要做特殊的处理，包括非平衡外推，反弹，对于边界上的每个格点之间的计算也没有数据的依赖性，因此，也可以让每个线程负责一个格点的计算；

c)CUDA的线程模型设计：根据GPU存储器的访问特点、GPU的寄存器资源及并行算法来设计CUDA的线程模型，确定总线程数和总线程块数；

d)求解离散方程和边界处理的GPU内核代码编写。

Images