CN103425833B - 一种基于熵格子波尔兹曼模型的并行流体计算实现方法 - Google Patents

Abstract

本发明公开了一种基于熵格子波尔兹曼模型的并行流体计算实现方法，提出了基于目前通用计算领域内主流显卡nVIDIA的GPU上熵格子波尔兹曼模型的并行实现方式，通过衡量流体并行计算加速比和比较每秒更新格子数的性能指标，应用nVIDIA显卡的GPU上统一计算设备编程架构（CUDA）实现ELBM模拟时间比CPU上模拟时间缩短三分之一；直接近似求解参数α的方法比迭代方法更加有效，即平均可以减少31.7%的时间。本发明可以充分利用系统的硬件资源，并且从实际操作层面验证了熵格子波尔兹曼模型的并行化计算方式，从而显著提供整个流体计算的效率。

Description

一种基于熵格子波尔兹曼模型的并行流体计算实现方法

技术领域

本发明属于计算机和流体力学交叉领域，特别涉及一种基于熵格子波尔兹曼模型的并行流体计算实现方法。

背景技术

传统意义上的计算流体力学(CFD)是指用数值方法求解Navier-Stokes(N-S）方程或Euler方程,属于典型的计算密集型学科,对计算能力和资源有较高要求。一般来说,对全机进行CFD计算(包括气动弹性、推进系统影响和控制系统作用），采用N-S方程的浮点运算量达到1015～1018flops,然而只有在计算机集群甚至超级计算机上才能达到1014flops,这在单颗CPU上计算也是异常耗时的。高效率低成本地在计算机上通过并行计算完成对机翼的CFD计算,采用N-S方程的运算量大小通常也一直是CFD研究的关键问题。近5年来,基于图形处理器(GPU)的通用计算技术的发展,如nVIDIA的统一计算设备架构(ComputeUnifiedDeviceArchitecture,CUDA),AMD的流体处理器CTM(CloseTotheMetal)等,为解决CFD计算问题提供了一种有利途径。从典型的CPU和GPU的性能对比，可以看出，GPU具有高出CPU10倍以上的潜在计算能力。充分利用GPU的高浮点计算性能，在个人计算机(PC)上快速进行一般规模(如机翼结构)的CFD计算，或者在由多块GPU组成的高性能计算机(HPC)上完成大规模(如全机)的CFD计算任务，将会对CFD研究和应用产生巨大的促进作用，然而相关研究却一直没有出现。

格子Boltzmann方法（LatticeBoltzmannmethod，LBM）是目前计算流体中最新的研究方法之一，它采用微观粒子模型，用离散的模型模拟连续的流体运动，不同于传统的N-S方程求解方法，具有很高的并行性。LBM是基于统计物理，并以极其简单的形式描述粒子的微观行为，但在宏观层次上正确反映流体的运动。由于它计算简单、本质并行和易于边界处理的优点，使得在这十几年里，在许多领域的各种数值问题上取得很大的成功。目前LBM方法的成功，被称为不可压缩流体计算领域的一场革命。LBM对复杂边界和复杂流场，如湍流、多相流、多孔多介质流、化学反应流、燃烧等问题的成功模拟，为我们展示了广阔的应用前景，为计算流体力学的数值模拟开辟了一条革命性的道路。尽管LBM取得了很大的成功，但是在处理低粘性流动问题的时候，LBM存在较为严重的数值不稳定性问题，标准的LBM在计算雷诺数（Reynolds）为2000的时候，顶盖驱动流就已很不稳定，这就限制了高雷诺数Reynolds流体流动的应用。

近年来发展的一些方法能部分解决此问题，其中有多松弛方法（MRT），熵格子波尔兹曼方法（ELBM）和基于Ehrenfest的粗粒度方法。ELBM主要思想是：在给出时空离散的格子系统中， $f_i(x+e_i{\mathrm{\δ}}_t,t+{\mathrm{\δ}}_t)=f_i(x,t)+\mathrm{\γ}[f_i^{\mathrm{eq}}(x,t)-f_i(x,t)],$ 其中，i为空间速度分布模型DdQq中第i个方向(0≤i≤q-1,q为方向数目，d为空间维数)、f_i为单粒子分布函数、为平衡态分布函数、γ为松弛时间参数和e_i为离散格子速度分量。对于格子系统中，最主要的问题就是选取合适的平衡态分布，使其在统计意义下满足系统熵最大，其中离散H熵函数定义为：这里ω_i为Gauss-Hermiter积分权数。对于标准LBM，松弛时间参数γ一般设定为常数；而对于ELBM，松弛时间参数γ必须调节，以使其满足H函数理论，使得H函数满足极值条件。H函数单调性约束由两个过程组成：1）单粒子分布函数偏离平衡态差量此过程中H函数保持为常数；2）耗散步，导致熵函数增加，在ELBM中松弛时间参数定义为γ=αβ，这里参数α,β定义为τ是标准LBM松弛时间。特别注意的是，当系统处于平衡时，α的极限值为2。

然而，在ELBM计算松弛时间参数过程中，其计算量将变得较多且耗时长，本发明针对此问题提出缩短计算时间、提高计算效率的综合解决方案。

发明内容

本发明提供了一种基于熵格子波尔兹曼模型的并行流体计算实现方法，其目的在于，克服现有技术中ELBM计算松弛时间参数过程中，计算量将变得较多且耗时长的问题。

一种基于熵格子波尔兹曼模型的并行流体计算实现方法，包括以下步骤：

步骤1：初始化参数及分配CPU与GPU内存空间；

根据应用场景中的粒子密度ρ、粒子粘度ν和选用的雷诺数Re对参数进行初始化及分配CPU与GPU内存空间；

初始化参数包括空间速度分布模型中的粒子速度u=0.1、特征线度L=1、松弛时间参数α₀=2、粒子密度分布函数f={f_i=0}_i=0,...,q-1，q是空间速度分布模型的方向数以及平衡态分布函数 $f^{\mathrm{eq}}={\left\{f_i^{\mathrm{eq}}\right\}}_{i=0,...,q-1},$ $f_i^{\mathrm{eq}}={\mathrm{\&rho;\&omega;}}_i\underset{j=1}{\overset{d}{\mathrm{\&Pi;}}}(2-\sqrt{1+u_j^2}){\left(\frac{{2u}_j+\sqrt{1+3u_j^2}}{1-u_j}\right)}^{\frac{v_{\mathrm{ij}}}{v}},$ 其中，d为空间维数，u_j=0为第j维的初始速度值，ω_i为Gauss-Hermiter积分权数，0<ω_i<1，从0、+1或-1中取值；

以空间速度分布模型中速度数组、密度分布数组以及平衡态密度分布数组对应的方向维度大小作为方向数，数组类型定义为浮点型，累加所有方向数，以累加的方向数乘以浮点型数组的字节长度得到的字节长度分配CPU与GPU内存空间，即对CPU和GPU分配的内存空间大小一样；

步骤2：将CPU内存中存储的空间速度分布模型中的速度数组、密度分布数组以及平衡态密度分布数组数据转移到GPU的全局内存Globalmemory；

步骤3：分配共享内存Sharedmemory空间；

按照步骤1所述的分配CPU或GPU内存空间大小对共享内存Sharedmemory空间进行分配；

步骤4：利用nVIDIAGPU编程语言CUDA实现三个函数；

1）H-α求解器函数；

该函数定义为设备端函数，修饰符为__device___，将在设备GPU端被调用，且在设备上运行，用于调整和计算松弛时间参数α；

采用泰勒一次展开法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}-\frac{F\left({\mathrm{\&alpha;}}^{*}\right)}{F^{\&prime;}\left({\mathrm{\&alpha;}}^{*}\right)}$

其中，F(α)=H(f+αΔ)-H(f)，α^*为前一次计算得到的α，α初始值为α₀，f_i(x,t)记为f_i，x和t分别表示空间中的粒子及该粒子所在的迭代时间，f_i(x,t)是粒子在迭代时间t下在位置x时的粒子概率值，F′(α^*)是F(α^*)对α^*进行求导；

2）传播碰撞内核函数；

该函数定义为设备端函数，修饰符为__global__，将在主机CPU端被调用，在设备GPU上运行；

在每次迭代时，通过标准的LBM模型中碰撞、传播表达式和每步调整后的松弛时间参数计算此迭代时间步下的粒子速度、粒子密度函数值以及粒子平衡态分布函数值；

3）边界处理内核函数；

该函数定义为设备端函数，修饰符为__global__，将在主机CPU端被调用，在设备GPU上运行；

采用反弹格式，从流体射向边界的微粒，被边界反射，以被碰撞璧前的速度即相同的速率沿原路返回；

所述修饰符是设置函数的所属类型；

步骤5：依据设定的精度要求和设定的迭代次数，在每次迭代过程中依次调用步骤4实现的三个函数，计算出密度函数值以及平衡态密度函数值，直到前后两次计算结果的差值满足精度要求或者到达迭代次数，迭代结束；

步骤6：迭代结束后，将GPU全局内存中存储的步骤5计算得到的密度函数值和平衡密度函数值转移到CPU内存，并且释放GPU内存空间；

步骤7：根据密度函数值以及平衡密度函数值，利用paraview软件得到流函数等值线图，同时释放主机内存空间，完成并行流体计算。

采用泰勒二次展开法实现所述步骤4中的H-α求解器函数；

泰勒二次展示法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}+\frac{{-F}_1\left({\mathrm{\&alpha;}}^{*}\right)+\sqrt{F_1{\left({\mathrm{\&alpha;}}^{*}\right)}^2-4F_2\left({\mathrm{\&alpha;}}^{*}\right)F\left({\mathrm{\&alpha;}}^{*}\right)}}{2F_2\left({\mathrm{\&alpha;}}^{*}\right)}$

F₁(α)=H′(f+αΔ)Δ，α^*为前一次计算得到的α，α初始值为α₀。

所述步骤1中对CPU和GPU内存空间进行分配采用Malloc函数以及cudaMalloc函数。

通过cudaMemcpyHostToDevice函数将CPU内存数据转移到GPU的全局内存Globalmemory。

通过cudaMemcpyDeviceToHost函数将GPU全局内存数据即密度函数值和平衡密度函数值转移到CPU内存。

有益效果

本发明的技术主要构思是在并行计算流体实现过程中，提出了基于目前通用计算领域内主流显卡nVIDIA的GPU上熵格子波尔兹曼模型的并行实现方式，通过衡量流体并行计算加速比和比较每秒更新格子数的性能指标，本发明的明显效果如下：

1)应用nVIDIA显卡的GPU上统一计算设备编程架构（CUDA）实现ELBM模拟时间比CPU上模拟时间缩短三分之一；

2)通过比较本发明中三种熵格子波尔兹曼模型的松弛时间参数求解方式，得出直接近似求解参数α的方法比迭代方法更加有效，即平均可以减少31.7%的时间。

综上所述本发明可以充分利用系统的硬件资源，并且从实际操作层面验证了熵格子波尔兹曼模型的并行化计算方式，从而显著提供整个流体计算的效率。更进一步来说，本发明通过利用熵格子波尔兹曼模型并行流体计算达到高稳定和高精度的要求，克服了ELBM计算松弛时间参数过程中，计算量将变得较多且耗时长的问题。

附图说明

图1为基于熵格子波尔兹曼模型的并行流体计算实现方法流程图；

图2为粒子二维速度空间分布模型D2Q9；

图3为求解器三种方法的计算时间比较图；

图4为雷诺数1000的等值流线比较图，其中，图(a)为本发明ELBM在GPU上模拟后的流函数等值线，图(b)为参考文献[1]ELBM在CPU上模拟后的流函数等值线，图(c)为参考文献[1]LBM在CPU上模拟后的流函数等值线。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述。

图1为基于熵格子波尔兹曼模型的并行流体计算实现方法流程图；

熵格子波尔兹曼模型是在标准的格子波尔兹曼模型的基础上选取合适的平衡态分布函数，使其在统计意义下满足系统熵最大，其中对于离散型H熵函数一般定义为：这里ω_i为Gauss-Hermiter积分权数。同时，对于熵格子波尔兹曼方程而言，松弛时间参数γ必须调整，以使其满足H函数理论，使得H函数满足极值条件。H函数单调性约束由两个过程组成：1）单粒子分布函数偏离平衡态差量此过程中H函数保持为常数；2）耗散步，导致熵函数增加，在ELBM中松弛时间参数定义为γ=αβ，这里参数α,β定义为τ是标准LBM松弛时间。

并行流体计算具体实现：

一.模拟场景描述

在满足熵格子波尔兹曼模型的前提下，并行算法的目标在于最小化模拟流体计算的运行时间，其工作重点在通过将密集型计算过程转移到GPU上较多线程执行来优化流体计算的应用程序。与CPU上多线程的熵格子波尔兹曼模型并行不同，GPU上粒子碰撞和传播过程可以在较CPU上百倍甚至千倍的线程上进行密度函数、平衡态分布函数以及参数计算，而且CPU和GPU上的相互通信时间代价是固定的。

首先，模拟场景中定义相关概念如下：对于二维顶盖驱动模型，采用空间速度分布模型D2Q9，如图2所示，定义二维规则方形的x轴方向格子数为nx，y轴方向格子数为ny。三维情形下，则增加z轴方向格子数nz。由于我们在模拟不同格子数下的熵格子波尔兹曼模型的并行计算性能，所以为便于对于不同情形下方便操作，实际模拟中的nx,ny均可配置。对于粒子9个不同的速度方向矢量e_i(0≤i<9)定义如下：

$e_i=\left\{\begin{array}{c}\left(\mathrm{0,0}\right),\\ (\cos \frac{(i-1)\mathrm{\&pi;}}{2},\sin \frac{(i-1)\mathrm{\&pi;}}{2}),\\ (\sqrt{2}\cos (\frac{\mathrm{\&pi;}}{4}+\frac{(i-5)\mathrm{\&pi;}}{2}),\sqrt{2}\sin (\frac{\mathrm{\&pi;}}{4}+\frac{(i-5)\mathrm{\&pi;}}{2}))\end{array}\right.$

同时，Gauss-Hermiter积分权数ω_i(0≤i<9)定义为：

${\mathrm{\&omega;}}_0=\frac{4}{9},{\mathrm{\&omega;}}_1={\mathrm{\&omega;}}_2={\mathrm{\&omega;}}_3={\mathrm{\&omega;}}_4=\frac{1}{9},{\mathrm{\&omega;}}_5={\mathrm{\&omega;}}_6={\mathrm{\&omega;}}_7={\mathrm{\&omega;}}_8=\frac{1}{36}.$

对于二维格子坐标(x,y),其中0≤x≤nx-1,0≤y≤ny-1。对应格子下的粒子速度数组u[i][j][k](0≤i<nx,0≤j<ny,0≤k<9)，对应格子下的粒子密度分布函数和平衡态分布函数数组分别为f[i][j][k],feq[i][j][k](0≤i<nx,0≤j<ny,0≤k<9)。

流体密度和黏度分别为ρ,υ，特征线度为L，顶盖速度为u，雷诺数定义为定义为γ=αβ，这里参数α,β定义为δ_t为间隔时间步长，取值为1，τ是标准LBM松弛时间。

二、算法设计

对于GPU上的标准格子波尔兹曼模型而言，不同格子上的粒子碰撞和传播过程均可以不同的线程进行控制和计算；具体来说，根据对应GPU上的线程ID对应不同的格子颗粒，同时进行计算粒子的碰撞后密度分布函数和平衡态分布函数。而对于熵格子波尔兹曼模型并行计算而言，松弛时间参数也可以通过GPU的线程进行计算和调整。对于调整松弛时间参数的计算方式，众所周知，它将对应求解相应粒子下的非线性方程组。我们分别通过牛顿迭代法、泰勒二次展开法以及泰勒三次展开法可进行求解，根据粒子波尔兹曼方程、数据分析以及数值计算知识，可以得到具体计算公式：

1）牛顿迭代法公式为：

${\mathrm{\&alpha;}}_{n+1}={\mathrm{\&alpha;}}_n-\frac{F\left({\mathrm{\&alpha;}}_n\right)}{F^{\&prime;}\left({\mathrm{\&alpha;}}_n\right)}$

这里F(α)=H(f+α(f^eq-f))-H(f)=H(f+αΔ)-H(f),F′(α)为对参数α求导；

2）采用泰勒一次展开法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}-\frac{F\left({\mathrm{\&alpha;}}^{*}\right)}{F^{\&prime;}\left({\mathrm{\&alpha;}}^{*}\right)}$

F(α)=H(f+αΔ)-H(f)，Δ=f^eq-f，α^*为前一次计算得到的α，α初始值为α₀，即为2；

3）泰勒二次展示法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}+\frac{{-F}_1\left({\mathrm{\&alpha;}}^{*}\right)+\sqrt{F_1{\left({\mathrm{\&alpha;}}^{*}\right)}^2-{4F}_2\left({\mathrm{\&alpha;}}^{*}\right)F\left({\mathrm{\&alpha;}}^{*}\right)}}{{2F}_2\left({\mathrm{\&alpha;}}^{*}\right)}$

F₁(α)=H′(f+αΔ)Δ，α^*为前一次计算得到的α，α初始值为α₀，即为2。

一种基于熵格子波尔兹曼模型的并行流体计算实现方法，包括以下步骤：

步骤1：初始化参数及分配CPU与GPU内存空间；

根据应用场景中的粒子密度ρ、粒子粘度ν和选用的雷诺数Re对参数进行初始化及分配CPU与GPU内存空间；

初始化参数包括空间速度分布模型中的粒子速度u=0.1、特征线度L=1、松弛时间参数α₀=2、粒子密度分布函数f={f_i=0}_i=0,...,q-1以及平衡态分布函数 $f^{\mathrm{eq}}={\left\{f_i^{\mathrm{eq}}\right\}}_{i=0,...,q-1},$ $f_i^{\mathrm{eq}}={\mathrm{\&rho;\&omega;}}_i\underset{j=1}{\overset{d}{\mathrm{\&Pi;}}}(2-\sqrt{1+u_j^2}){\left(\frac{{2u}_j+\sqrt{1+{3u}_j^2}}{1-u_j}\right)}^{\frac{v_{\mathrm{ij}}}{v}},$ 其中，q是空间速度分布模型的方向数，d为空间维数，u_j=0为第j维的初始速度值，ω_i为Gauss-Hermiter积分权数，0<ω_i<1，从0、+1或-1中取值；

以空间速度分布模型中速度数组、密度分布数组以及平衡态密度分布数组对应的方向维度大小作为方向数，数组类型定义为浮点型，累加所有方向数，以累加的方向数乘以浮点型数组的字节长度得到的字节长度分配CPU与GPU内存空间，即对CPU和GPU分配的内存空间大小一样；

步骤2：将CPU内存中存储的空间速度分布模型中速度数组、密度分布数组以及平衡态密度分布数组数据转移到GPU的全局内存Globalmemory；

步骤3：分配共享内存Sharedmemory空间；

按照步骤1所述的分配CPU或GPU内存空间大小对共享内存Sharedmemory空间进行分配；

步骤4：利用nVIDIAGPU编程语言CUDA实现三个函数；

1）H-α求解器函数；

该函数定义为设备端函数，修饰符为__device___，将在设备GPU端被调用，且在设备上运行，用于调整和计算松弛时间参数α；

采用泰勒一次展开法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}-\frac{F\left({\mathrm{\&alpha;}}^{*}\right)}{F^{\&prime;}\left({\mathrm{\&alpha;}}^{*}\right)}$

其中，F(α)=H(f+αΔ)-H(f)，α^*为前一次计算得到的α，α初始值为α₀，f_i(x,t)记为f_i，x和t分别表示空间中的粒子及该粒子所在的迭代时间，f_i(x,t)是粒子在迭代时间t下在位置x时的粒子概率值，F′(α^*)是F(α^*)对α^*进行求导；

2）传播碰撞内核函数；

该函数定义为设备端函数，修饰符为__global__，将在主机CPU端被调用，在设备GPU上运行；

在每次迭代时，通过标准的LBM模型中碰撞、传播表达式和每步调整后的松弛时间参数计算此迭代时间步下的粒子速度、粒子密度函数值以及粒子平衡态分布函数值；

3）边界处理内核函数；

该函数定义为设备端函数，修饰符为__global__，将在主机CPU端被调用，在设备GPU上运行；

采用反弹格式，从流体射向边界的微粒，被边界反射，以被碰撞璧前的速度即相同的速率沿原路返回；

所述修饰符是设置函数的所属类型；

步骤5：依据设定的精度要求和设定的迭代次数，在每次迭代过程中依次调用步骤4实现的三个函数，计算出密度函数值以及平衡态密度函数值，直到前后两次计算结果的差值满足精度要求或者到达迭代次数，迭代结束；

步骤6：迭代结束后，将GPU全局内存中存储的步骤5计算得到的密度函数值和平衡密度函数值转移到CPU内存，并且释放GPU内存空间；

步骤7：根据密度函数值以及平衡密度函数值，利用paraview软件得到流函数等值线图，同时释放主机内存空间，完成并行流体计算。

采用泰勒二次展开法实现所述步骤4中的H-α求解器函数；

泰勒二次展示法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}+\frac{{-F}_1\left({\mathrm{\&alpha;}}^{*}\right)+\sqrt{F_1{\left({\mathrm{\&alpha;}}^{*}\right)}^2-{4F}_2\left({\mathrm{\&alpha;}}^{*}\right)F\left({\mathrm{\&alpha;}}^{*}\right)}}{{2F}_2\left({\mathrm{\&alpha;}}^{*}\right)}$

F₁(α)=H′(f+αΔ)Δ，α^*为前一次计算得到的α，α初始值为2。

所述步骤1中对CPU和GPU内存空间进行分配采用Malloc函数以及cudaMalloc函数。熵格子波尔兹曼模型的并行计算程序运行过程中，主要是通过主机端调用设备端的kernel进行粒子相关物理量的计算，其中kernel包括碰撞和传播、边界处理以及参数求解器。而对于整个并行计算的性能指标，我们通过加速比以及每秒钟更新格子数进行测评，它们的计算公式如下：

并行算法加速比=串行算法运算时间/并行算法运算时间；

每秒钟更新格子数=各维度格子数乘积*迭代语句数*10^-6/并行计算时间(MULPS)；

我们根据三种不同的并行方式求解松弛时间参数，最终得出的牛顿迭代法N.R.、泰勒一次展开法FOE以及泰勒二次展开法SOE的并行算法加速比分别为3.08、3.16和3.18。同时，对于三种方式求解松弛时间参数的时间占总体计算的时间百分比如图3所示。

对于另一个性能指标每秒钟更新格子数，通过不同格子数不同线程的情形下进行实际计算，可以得出对应于具体的nVIDIA显卡而言，格子数和线程数的选取将很大影响此性能指标。特别的，对于nVIDIA的G210显卡而言，当格子数为256*256，线程数为256时，我们得出的每秒钟更新格子数为57.6MLUPS。这表明程序的运行性能将得到很大的提高，从而对显卡的硬件资源也得到充分的利用。

实施例1：

对于流体密度和黏度分别为ρ,υ，特征距离为d，顶盖速度为u，我们调整物理量可计算雷诺数Re为1000的情况下，通过上述具体实施方式模拟流体在顶盖驱动下的运动情形。其流函数等值线见图4所示，与此同时，给出了与相关文献[1]同雷诺数下的流体模拟的比较。从图和表1中可以看出本发明所述方法在主涡及次涡的位置清晰可见的情况下，ELBM模拟速度基于NVIDIAG210GPU上明显优于在CPU上运行，平均加速比为3，其具体对应数据参见表1。

表1

其中上述提到的文献包括如下：

[1]TosiF,UbertiniS,SucciS,KarlinI.Optimizationstrategiesfortheentropiclatticeboltzmannmethod.JournalofScientificComputing2007;30(3):369–387。

Claims (5)

1.一种基于熵格子波尔兹曼模型的并行流体计算实现方法，其特征在于，包括以下步骤：

步骤1：初始化参数及分配CPU与GPU内存空间；

根据应用场景中的粒子密度ρ、粒子粘度ν和选用的雷诺数Re对参数进行初始化及分配CPU与GPU内存空间；

初始化参数包括空间速度分布模型中的粒子速度u＝0.1、特征线度L＝1、松弛时间参数α₀＝2、粒子密度分布函数f＝{f_i＝0}_{i＝0,...,q-1}，q是空间速度分布模型的方向数以及平衡态分布函数 $f^{eq}={\left\{f_i^{eq}\right\}}_{i=0,\mathrm{...},q-1},f_i^{eq}={\mathrm{\&rho;\&omega;}}_i\underset{j=1}{\overset{d}{\&Pi;}}(2-\sqrt{1+u_j^2}){\left(\frac{2u_j+\sqrt{1+3u_j^2}}{1-u_j}\right)}^{\frac{v_{ij}}{v}},$ 其中，d为空间维数，u_j＝0为第j维的初始速度值，ω_i为Gauss-Hermiter积分权数，0＜ω_i＜1，从0、+1或-1中取值；

以空间速度分布模型中速度数组、密度分布数组以及平衡态密度分布数组对应的方向维度大小作为方向数，数组类型定义为浮点型，累加所有方向数，以累加的方向数乘以浮点型数组的字节长度得到的字节长度分配CPU与GPU内存空间，即对CPU和GPU分配的内存空间大小一样；

步骤2：将CPU内存中存储的空间速度分布模型中的速度数组、密度分布数组以及平衡态密度分布数组数据转移到GPU的全局内存Globalmemory；

步骤3：分配共享内存Sharedmemory空间；

按照步骤1所述的分配CPU或GPU内存空间大小对共享内存Sharedmemory空间进行分配；

步骤4：利用nVIDIAGPU编程语言CUDA实现三个函数；

1)H-α求解器函数；

该函数定义为设备端函数，修饰符为__device___，将在设备GPU端被调用，且在设备上运行，用于调整和计算松弛时间参数α；

采用泰勒一次展开法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}-\frac{F\left({\mathrm{\&alpha;}}^{*}\right)}{F^{\&prime;}\left({\mathrm{\&alpha;}}^{*}\right)}$

其中，F(α)＝H(f+αΔ)-H(f)，Δ＝f^eq-f＝{f_i ^eq-f_i}_{i＝0,...,q-1}，α^*为前一次计算得到的α，α初始值为α₀，f_i(x,t)记为f_i，x和t分别表示空间中的粒子及该粒子所在的迭代时间，f_i(x,t)是粒子在迭代时间t下在位置x时的粒子概率值，F′(α^*)是F(α^*)对α^*进行求导；

2)传播碰撞内核函数；

该函数定义为设备端函数，修饰符为__global__，将在主机CPU端被调用，在设备GPU上运行；

3)边界处理内核函数；

该函数定义为设备端函数，修饰符为__global__，将在主机CPU端被调用，在设备GPU上运行；

所述修饰符是设置函数的所属类型；

步骤5：依据设定的精度要求和设定的迭代次数，在每次迭代过程中依次调用步骤4实现的三个函数，计算出密度函数值f_i以及平衡态密度函数值f_i ^eq，直到前后两次计算结果的差值满足精度要求或者到达迭代次数，迭代结束；

步骤6：迭代结束后，将GPU全局内存中存储的步骤5计算得到的密度函数值和平衡密度函数值转移到CPU内存，并且释放GPU内存空间；

步骤7：根据密度函数值以及平衡密度函数值，利用paraview软件得到流函数等值线图，同时释放主机内存空间，完成并行流体计算。

2.根据权利要求1所述的基于熵格子波尔兹曼模型的并行流体计算实现方法，其特征在于，采用泰勒二次展开法实现所述步骤4中的H-α求解器函数；

泰勒二次展示法计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}^{*}+\frac{-F_1\left({\mathrm{\&alpha;}}^{*}\right)+\sqrt{F_1{\left({\mathrm{\&alpha;}}^{*}\right)}^2-4F_2\left({\mathrm{\&alpha;}}^{*}\right)F\left({\mathrm{\&alpha;}}^{*}\right)}}{2F_2\left({\mathrm{\&alpha;}}^{*}\right)}$

F₁(α)＝H′(f+αΔ)Δ，α^*为前一次计算得到的α，α初始值为α₀。

3.根据权利要求1所述的基于熵格子波尔兹曼模型的并行流体计算实现方法，其特征在于，所述步骤1中对CPU和GPU内存空间进行分配采用Malloc函数以及cudaMalloc函数。

4.根据权利要求1所述的基于熵格子波尔兹曼模型的并行流体计算实现方法，其特征在于，通过cudaMemcpyHostToDevice函数将CPU内存数据转移到GPU的全局内存Globalmemory。

5.根据权利要求1所述的基于熵格子波尔兹曼模型的并行流体计算实现方法，其特征在于，通过cudaMemcpyDeviceToHost函数将GPU全局内存数据即密度函数值和平衡密度函数值转移到CPU内存。