CN102945295A - 一种格子玻尔兹曼方法的并行加速方法及系统 - Google Patents

Abstract

本发明公开了一种格子玻尔兹曼方法的并行加速方法及系统，克服目前使用X86服务器集群来进行LBM处理在时间等消耗较大的缺陷，该方法中：CPU端确定网格所有格点上的宏观参量及其初始值，定义数据结构和存储方式用于存储每个格点各个方向的平衡分布函数和每个格点的宏观参量，根据宏观参量计算出所有格点上各方向的平衡分布函数，设置MIC卡内核的线程执行配置；MIC端根据宏观参量及其初始值、线程执行配置及网格所有格点上各方向的平衡分布函数进行迁移和碰撞的并行计算和并行的边界处理，并根据网格所有格点各方向的平衡分布函数获得收敛状态。本申请的实施例提高了格子Boltzmann方法的处理性能，满足了流体模拟的需求。

Description

一种格子玻尔兹曼方法的并行加速方法及系统

技术领域

本发明涉及计算流体力学领域，尤其涉及一种格子玻尔兹曼方法的并行加速方法及系统。

背景技术

格子玻尔兹曼方法(Lattice Boltzmann Method，LBM)是介于微观分子动力学方法和基于连续介质假设的宏观方法之间的一种介观方法。该方法与传统的流体模拟方法不同，它基于分子动理论，通过跟踪粒子分布函数的输运而后对分布函数求矩来获得宏观平均特性。

格子Boltzmann方法的动理论特性使得它在模拟许多复杂流动的时候更有效，如多孔介质流动、悬浮流、多相流、多组分流等。格子Boltzmann方法具有天生的并行特性，以及边界处理简单、程序易于实施等优点。

采用LBM方法求解物理问题时的基本过程如图1所示。

对于一个特定的物理问题，首先进行如图1中步骤S110至步骤S130的辅助步骤。

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

步骤S120，进行网格划分，确定节点。假设网格划分大小为NX*NY，其中，NX为x方向的大小，NY为y方向的大小。

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

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

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

步骤S150，求解离散后的控制方程，例如，采用迁移碰撞规则求解LBGK方程。

步骤S160，根据边界条件，在相应边界格点上实施边界处理格式。

步骤S170，基于不同格子Boltzmann模型的宏观量的定义法则，计算各格点上的宏观参量。

步骤S180，判断计算是否收敛，若计算收敛，则转步骤S190；否则返回步骤S150继续求解。

步骤S190，输出计算结果。

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

$f_i(\stackrel{\→}{x}+{\stackrel{\→}{e}}_i\mathrm{\δt},t+\mathrm{\δt})=f_i(\stackrel{\→}{x},t)-\frac{1}{\mathrm{\τ}}[f_i(\stackrel{\→}{x},t)-f_i^{\mathrm{eq}}(\stackrel{\→}{x},t)],i=0...N$ 式(1)

这里，N为网格大小，f_i为粒子分布函数，代表时间t在处存在以微观速度

运动的粒子的概率。松弛时间τ代表达到局部平衡的速率，与流体的动粘性系数相关。平衡分布函数

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

式(2)

表示处的离子密度，

表示离子运动速度。

其中，在D2Q9模型中：

式(3)

流体密度和速度则可由粒子分布函数依如下表达式计算获得。

$\underset{i}{\mathrm{\Σ}}{f}_i=\mathrm{\ρ}$ 式(4)

$\underset{i}{\mathrm{\Σ}}{f}_i{e}_i=\mathrm{\ρu}$ 式(5)

其中，ρ表示密度，u表示速度。

离散速度矢量

和粒子分布函数的个数N，取决于所选定的格子-Boltzmann模型，在D2Q9模型中

有9个分量，对应的粒子分布函数的个数也为9个，见图2。

至强融核(Many Integrated Core，MIC)是英特尔(Intel)公司推出的一款众核处理器，跟通用的多核至强处理器相比，MIC众核架构具有更小的内核和硬件线程，众核处理器计算资源密度更高，片上通信开销显著降低，更多的晶体管和能量，能够胜任更为复杂的并行应用。

OpenMP是用于共享内存并行系统上的多线程程序设计的一套指导性注释，在MIC平台上同样支持OpenMP编程模型，这样减少了MIC并行程序的开发周期，对传统的并行编程语言具有很好的支持，因此，可以利用OpenMP并行编程模型在MIC平台上快速实现高性能应用并行软件，快速得到性能的提高。

LBM方法需要大量的计算。以方腔的计算为例，假设网格大小为1024*1024，迭代10000次，每个网格点都要在一次迭代中做一次迁移和碰撞计算，在主频为3.00GHz的QuadCore Intel Xeon X5450进行计算需要数小时的时间，而计算更大的网格迭代更多次数时需要数天的时间，这严重影响了LBM方法的性能。

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

发明内容

本发明所要解决的技术问题是克服目前使用X86服务器集群来进行LBM处理在时间、电力以及维护上消耗较大的缺陷。

为了解决上述技术问题，本发明提供了一种格子玻尔兹曼方法的并行加速方法，CPU端执行如下处理：

CPU端确定网格所有格点上的宏观参量及其初始值，定义数据结构和存储方式用于存储每个格点各个方向的平衡分布函数和每个格点的宏观参量，并根据所述宏观参量计算出所有格点上各个方向的平衡分布函数，并设置至强融核(MIC)卡内核的线程执行配置；

至强融核(MIC)端根据所述宏观参量及其初始值、所述线程执行配置以及所述网格所有格点上各个方向的平衡分布函数进行迁移和碰撞的并行计算和并行的边界处理，并根据所述网格所有格点各个方向的平衡分布函数获得收敛状态并传递给所述CPU端。

优选地，所述宏观参量包括密度、速度、参考长度、雷诺数和黏性系数。

优选地，所述CPU端设置MIC卡内核的线程执行配置的步骤，包括：

所述CPU端根据所述MIC卡的核数设置并行执行的线程数。

优选地，所述并行的边界处理包括：

反弹法或者非平衡外推法的边界处理。

优选地，所述MIC端根据所述CPU端确定的网格所有格点上的宏观参量及其初始值，计算出所有格点上各个方向的平衡分布函数。

优选地，所述MIC端的MIC卡内核中的每个线程计算一行网格点的迁移和碰撞处理。

本申请的实施例还提供了一种格子玻尔兹曼方法的并行加速系统，包括CPU端和至强融核(MIC)端，其中：

所述CPU端用于确定网格所有格点上的宏观参量及其初始值，定义数据结构和存储方式用于存储每个格点各个方向的平衡分布函数和每个格点的宏观参量，并根据所述宏观参量计算出所有格点上各个方向的平衡分布函数；以及设置至强融核(MIC)卡内核的线程执行配置；

所述MIC端用于根据所述宏观参量及其初始值、所述线程执行配置以及所述网格所有格点上各个方向的平衡分布函数进行迁移和碰撞的并行计算以及并行的边界处理；并用于根据所述网格所有格点各个方向的平衡分布函数获得收敛状态传递给所述CPU端。

优选地，所述宏观参量比如包括密度、速度、参考长度、雷诺数和黏性系数。

优选地，所述CPU端用于根据MIC端的MIC卡的核数设置并行执行的线程数。

优选地，所述MIC端用于进行反弹法或者非平衡外推法的边界处理。

优选地，所述MIC端用于根据所述CPU端确定的网格所有格点上的宏观参量及其初始值，计算出所有格点上各个方向的平衡分布函数。

优选地，所述MIC端的MIC卡内核中的每个线程计算一行网格点的迁移和碰撞处理。

与现有技术相比，本申请的实施例提高了格子Boltzmann方法的处理性能，满足了流体模拟的需求，而且能降低功耗，减少机房构建成本和管理、运行、维护费用。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

附图用来提供对本发明技术方案的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明的技术方案，并不构成对本发明技术方案的限制。

图1为现有技术中LBM方法模拟求解的基本流程示意图。

图2为现有技术中D2Q9模型示意图。

图3为本发明实施例格子玻尔兹曼方法的并行加速方法的流程示意图。

图4为本发明实施例中的内核中内层循环的向量化的示意图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。本发明实施例以及实施例中的各个特征在不相冲突前提下的相互结合，均在本发明的保护范围之内。

如图3所示，本申请实施例的格子Boltzmann的并行加速方法主要包括如下步骤。

步骤S310，CPU端根据物理问题，根据物理问题进行网格划分，确定网格所有格点上的宏观参量及其初始值，传递给MIC端。其中该宏观参量比如包括密度、速度、参考长度、雷诺数和黏性系数等。

本申请的实施例，根据物理问题的要求将流场计算域进行网格划分，比如网格大小为NX*NY，NX为x方向，NY为y方向，网格上的节点数为N＝NX*NY。

步骤S320，CPU端定义数据结构和存储方式，用于存储每个格点各个方向的平衡分布函数和每个格点的宏观参量，并根据宏观参量计算出所有格点上各个方向的平衡分布函数，以此作为计算的初场。

步骤S330，CPU端设置内核的线程执行配置。CPU端根据MIC卡的核数设置并行执行迁移碰撞等迭代计算的线程数(比如MIC卡的核数为M，则迁移碰撞计算的线程数T可以等于4*M)，使得内核中的每个线程计算一行网格点的迁移和碰撞过程。

步骤S340，MIC端根据CPU端确定的宏观参量及其初始值、CPU端设置的线程执行配置以及CPU端确定的网格所有格点上各个方向的平衡分布函数，进行迁移、碰撞并行计算，利用#pragma ivdep实现如图4所示的内核中内层循环的向量化。

内核伪代码如下：

步骤S350，MIC端进行并行的边界处理，边界处理可以采用反弹法或者非平衡外推法等方法，边界处理时同样设计T个线程处理边界节点的分布函数计算。

步骤S360，判断是否迭代完成，完成则转步骤S370，否则返回步骤S340继续迭代。

步骤S370，MIC端根据分布函数并行求得速度、密度和流函数等宏观参量并把流场的收敛状态传递给CPU端。

步骤S380，CPU端输出流场的收敛状态。

本申请的实施例中，MIC端根据CPU端确定的宏观参量及其初始值，以及CPU端设置的线程执行配置使用相应的多线程并行地根据所有格点上各个方向的平衡分布函数，通过迁移、碰撞和边界处理的迭代处理，得到流场的收敛状态。其中，MIC端根据所有格点上各个方向的平衡分布函数，通过迁移、碰撞和边界处理，得到下一个时层的分布函数，然后通过迭代处理得到流场的收敛状态。其中，MIC端根据宏观参量及其初始值，使用一个求分布函数的内核，计算出所有格点上各个方向的平衡分布函数。

MIC端根据线程执行配置使用相应的多线程并行地根据流场当前时层的分布函数，通过迁移、碰撞和边界处理得到该流场下一时层的分布函数，具体过程可以是MIC端使用T个线程并行地对流体网格的N个格点根据初始的分布函数F_i ⁽⁰⁾或上一步计算得到的分布函数F_i ^(k)执行迁移、碰撞和边界处理，计算获取所述网格格点的下一时层的分布函数F_i ^(k+1)，所述i取0-b共b+1个值，分别表示格点上的b+1个方向的分布函数，所述k为大于等于1的整数。

本申请的实施例中，CPU端控制迭代次数，MIC端获得的最终结果也回传到CPU端。其中，CPU端控制迭代ITR次，即内核迭代调用ITR次，所述ITR为流体模拟中进行的迭代次数。

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

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

求解离散方程可以采用迁移碰撞的过程，宏观量统计、平衡分布函数计算和碰撞过程中对每个网格的计算之间没有任何依赖性。因此，可以让MIC中的每个线程负责一个网格划分中的一行网格点的计算，每行网格点的计算利用MIC上的向量化技术进一步加速。分布函数的迁移只涉及到该格点周围的其他格点，也可以通过单个线程对全局存储器中相关分布函数的读操作来实现。

在LBM算法中，对边界要做特殊的处理(非平衡外推法或者反弹法)，对于边界上的每个格点之间的计算也没有数据的依赖性，因此，可以利用OpenMP多线程负责边界格点的计算。OpenMP的线程模型根据MIC核心数目设置内核的线程数。

本申请实施例的格子玻尔兹曼方法的并行加速系统，包括CPU端和至强融核(MIC)端。

CPU端用于确定网格所有格点上的宏观参量及其初始值，定义数据结构和存储方式用于存储每个格点各个方向的平衡分布函数和每个格点的宏观参量，并根据宏观参量计算出所有格点上各个方向的平衡分布函数；以及设置至强融核(MIC)卡内核的线程执行配置；

MIC端用于根据宏观参量及其初始值、线程执行配置以及网格所有格点上各个方向的平衡分布函数，进行迁移和碰撞的并行计算；实现内核中内层循环的向量化，进行并行的边界处理；并根据CPU端确定的网格所有格点各个方向的平衡分布函数并行求得宏观参量；获得流场的收敛状态传递给CPU端。

上述的宏观参量比如包括密度、速度、参考长度、雷诺数和黏性系数。

本申请的实施例中，CPU端用于根据MIC端的MIC卡的核数设置并行执行的线程数。

本申请的实施例中，MIC端可以用于进行反弹法或者非平衡外推法的边界处理。

本申请的实施例中，MIC端用于根据CPU端确定的网格所有格点上的宏观参量及其初始值，计算出所有格点上各个方向的平衡分布函数。

本申请的实施例中，MIC端的MIC卡内核中的每个线程计算一行网格点的迁移和碰撞处理。

本申请的实施例通过CPU端与MIC端进行协同计算，从而满足了流体模拟的需求，降低了机房构建成本和管理、运行、维护费用。本申请的实施例将基本参数初始化计算放在CPU端执行，而把耗时并且并行性非常好的平衡分布函数计算、宏观量统计、离散方程求解和边界处理部分利用OpenMP技术进行并行化设计，使其在MIC端并行执行，CPU与MIC进行协同计算，最终实现加速格子Boltzmann方法。

虽然本发明所揭露的实施方式如上，但所述的内容仅为便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属领域内的技术人员，在不脱离本发明所揭露的精神和范围的前提下，可以在实施的形式及细节上进行任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (12)

1.一种格子玻尔兹曼方法的并行加速方法，包括：

CPU端确定网格所有格点上的宏观参量及其初始值，定义数据结构和存储方式用于存储每个格点各个方向的平衡分布函数和每个格点的宏观参量，并根据所述宏观参量计算出所有格点上各个方向的平衡分布函数，并设置至强融核(MIC)卡内核的线程执行配置；

至强融核(MIC)端根据所述宏观参量及其初始值、所述线程执行配置以及所述网格所有格点上各个方向的平衡分布函数进行迁移和碰撞的并行计算和并行的边界处理，并根据所述网格所有格点各个方向的平衡分布函数获得收敛状态并传递给所述CPU端。

2.根据权利要求1所述的方法，其中：

所述宏观参量包括密度、速度、参考长度、雷诺数和黏性系数。

3.根据权利要求1所述的方法，其中，所述CPU端设置MIC卡内核的线程执行配置的步骤，包括：

所述CPU端根据所述MIC卡的核数设置并行执行的线程数。

4.根据权利要求1所述的方法，其中，所述并行的边界处理包括：

反弹法或者非平衡外推法的边界处理。

5.根据权利要求1所述的方法，其中：

所述MIC端根据所述CPU端确定的网格所有格点上的宏观参量及其初始值，计算出所有格点上各个方向的平衡分布函数。

6.根据权利要求1所述的方法，其中：

所述MIC端的MIC卡内核中的每个线程计算一行网格点的迁移和碰撞处理。

7.一种格子玻尔兹曼方法的并行加速系统，包括CPU端和至强融核(MIC)端，其中：

所述CPU端用于确定网格所有格点上的宏观参量及其初始值，定义数据结构和存储方式用于存储每个格点各个方向的平衡分布函数和每个格点的宏观参量，并根据所述宏观参量计算出所有格点上各个方向的平衡分布函数；以及设置至强融核(MIC)卡内核的线程执行配置；

所述MIC端用于根据所述宏观参量及其初始值、所述线程执行配置以及所述网格所有格点上各个方向的平衡分布函数进行迁移和碰撞的并行计算以及并行的边界处理；并用于根据所述网格所有格点各个方向的平衡分布函数获得收敛状态传递给所述CPU端。

8.根据权利要求7所述的系统，其中：

所述宏观参量比如包括密度、速度、参考长度、雷诺数和黏性系数。

9.根据权利要求7所述的系统，其中：

所述CPU端用于根据MIC端的MIC卡的核数设置并行执行的线程数。

10.根据权利要求7所述的系统，其中：

所述MIC端用于进行反弹法或者非平衡外推法的边界处理。

11.根据权利要求7所述的系统，其中：

所述MIC端用于根据所述CPU端确定的网格所有格点上的宏观参量及其初始值，计算出所有格点上各个方向的平衡分布函数。

12.根据权利要求7所述的系统，其中：

所述MIC端的MIC卡内核中的每个线程计算一行网格点的迁移和碰撞处理。

Images