CN107515987A - 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 - Google Patents
基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 Download PDFInfo
- Publication number
- CN107515987A CN107515987A CN201710741065.5A CN201710741065A CN107515987A CN 107515987 A CN107515987 A CN 107515987A CN 201710741065 A CN201710741065 A CN 201710741065A CN 107515987 A CN107515987 A CN 107515987A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- gpu
- distribution function
- lattice
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，属于地下水数值模拟领域，该方法运行于带有若干个GPU的计算机或服务器，该方法包括：将地下水区域离散化为若干个网格，给定初始宏观物理量以及初始分布函数，并将网格平均分为若干块传入若干GPU；每个GPU调用多个线程同时进行初始化；迭代开始，每个GPU调用多个线程同时根据多松弛格子玻尔兹曼方程进行块内格点碰撞操作；各个GPU进行格点之间的迁移操作；各个GPU调用多个线程得到非边界格点的宏观物理量；格点的边界处理，得到边界格点的分布函数和宏观物理量；判断稳定性。本发明占地面积小、造价低、耗能低，加速明显。
Description
技术领域
本发明涉及地下水数值模拟领域，特别是指一种基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法。
背景技术
地下水数值模拟能够方便、快捷地在计算机上模拟复杂多变的地下水系统，是开展地下水系统研究必不可少的一种工具。格子玻尔兹曼方法进行地下水模拟既可以描述流体的内部细节，又无需构造大规模线性方程组，且具有算法步骤简单、易于网格划分、天然的并行性以及易处理的边界条件等优点。
多松弛格子玻尔兹曼模型是一种常用的计算模型，与单松弛模型相比，它允许多个方向具有不同的松弛时间，在收敛性、数值稳定性以及应用范围等方面具有更大的优势。但通过多松弛格子玻尔兹曼模型模拟地下水流动过程需要进行实际区域的网格划分，网格划分得越细，需要求解的网格数量也就越多，计算时间也随之增加，由于CPU计算核心有限，欲达到较好的加速效果，需要配制多个CPU，占地面积大，昂贵，能耗较高。
发明内容
本发明提供一种占地面积小、造价低、耗能低并且加速明显的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法。
为解决上述技术问题，本发明提供技术方案如下：
本发明提供一种基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，所述基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法运行于带有若干个GPU的计算机或服务器，所述基于多松弛格子玻尔兹曼模型的地下水流动模拟加速系统的加速方法，包括：
步骤1：将地下水区域离散化为若干个网格，每个网格格点给定其初始宏观物理量以及初始分布函数，并将所有网格平均分为若干块传入若干个已开辟好空间的GPU；
步骤2：每个GPU调用多个线程同时进行各个块内网格格点的初始宏观物理量及其对应的初始分布函数的初始化；
步骤3：迭代开始，每个GPU调用多个线程同时根据多松弛格子玻尔兹曼方程进行块内网格格点碰撞操作；
步骤4：各个GPU将块内网格格点的分布函数组合到0号GPU上，由0号GPU作为主线程进行格点之间的迁移操作，将之后再次获得的分布函数分发给其他GPU；
步骤5：各个GPU调用多个线程将块内各个网格格点的分布函数值求和，即得到网格非边界格点的宏观物理量；
步骤6：各个GPU分别进行块内网格格点的边界处理，得到网格边界格点的分布函数和宏观物理量；
步骤7：将步骤5得到的网格非边界格点的宏观物理量和步骤6得到的网格边界格点的宏观物理量均传回CPU，与开始迭代时的对应的网格格点的宏观物理量进行对比判断稳定性，若误差大于设定阈值，则继续迭代，执行步骤3，否则，迭代结束，释放所有GPU的内存。
进一步的，所述步骤1中，网格格点的初始分布函数f的数量由不同的DnQb模型来决定，n表示维度，b表示不同方向分布函数f的数量。
进一步的，所述步骤2包括：
给定地下水流动偏微分方程：
中渗透系数Kxx、Kyy和Kzz，单位体积流量W，孔隙介质的贮水率Ss的宏观物理量的初始值，其中，公式(1)中：Kxx、Kyy和Kzz分别为x、y和z方向上的渗透系数分量；h为所求宏观物理量，表示水头；W为单位体积流量，用来表示流入或流出流体区域的水量；Ss为孔隙介质的贮水率；t为时间；
根据公式(2)、(3)初始化每个网格格点的水头并根据该水头得出分布函数值，一般初始分布函数设为其平衡态分布函数，
fi＝fi (eq)(ρ0，u0，T0) (2)
其中，公式(2)中，ρ0、u0和T0是初始时刻的密度、速度、温度的宏观物理量，在地下水流动模拟中，宏观物理量为各个网格格点的水头h。
进一步的，所述步骤3中，多松弛格子玻尔兹曼方程表示为
fi(x+ciδt，t+δt)-fi(x，t)＝-Λij[fi(x，t)-fi (eq)(x，t)] (4)
其中Λ是一个b×b的矩阵，表示碰撞因子，经过一段时间t的变化，公式(4)改写为
其中m和m(eq)分别为分布函数f和平衡态分布函数f(eq)的矩空间，即m＝Mf，m(eq)＝Mf(eq)，S为松弛对角矩阵，S＝MΛM-1。
进一步的，所述步骤4中，所述迁移操作根据DnQb模型中分布函数的方向进行，即
fi(x+ciδt，t+δt)＝fi′(x，t) (6)
其中ci表示离散速度，对应DnQb模型中的方向坐标。
进一步的，所述步骤5中，根据DnQb模型，每个网格格点具有b个不同方向的分布函数，所以分布函数的和即为该网格非边界格点的宏观物理量，即
∑ifi＝ρ (7)
其中ρ表示任一宏观物理量，ρ表示为水头h。
进一步的，所述步骤6中，边界条件包括定水头边界、变水头边界以及无渗流边界，对于定水头边界和变水头边界，采用非平衡态外推格式来处理，即
fi(xb，t)＝fi (eq)(xb，t)+fi (ne)(xb，t) (8)
其中fi (ne)表示非平衡态部分，对于边界处有一定速度的边界来说，边界格点xb处的宏观uω速度和其相邻格点xf的宏观量密度ρf都是已知的，代入公式(8)中，得到速度边界的非平衡态外推格式，即
对于无渗流边界，将边界处格点的分布函数设为临近点对应位置网格边界格点的分布函数，即
fi(xb，t)＝fi(xf，t) (10)。
进一步的，所述步骤7中，通过当前宏观物理量与上一时刻宏观物理量误差的最大值，判断系统是否达到稳定，当误差小于设定阈值时，即
则认为系统稳定，否则继续迭代，直到特定时刻或者系统达到稳定状态。
进一步的，所述步骤3中，利用部分CUBLAS库以及CUDA的kernel函数结合OpenMP技术，实现其在多个GPU上的线程级并行。
进一步的，基于所述CUDA，通过设定线程块block和块内线程thread的数量，设定并行使用的GPU总线程数量，实现线程级并行；
线程块或线程之间在一定的内存中读取数据，内存包括线程本地内存和寄存、线程块内的共享内存以及块间网格内的全局内存、常量内存和纹理内存。
本发明具有以下有益效果：
本发明通过多松弛格子玻尔兹曼模型模拟了地下水流动过程并将其在具有多块GPU卡作为加速部件的服务器实现了并行，在占地更小、成本更低、能耗更小的情况下获得更大的加速效果，在保证精度的基础上具有较高的运行效率。
附图说明
图1为本发明的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法的流程图；
图2为本发明的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法的部分DnQb模型示意图；
图3为本发明的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法的实施例1的模拟区域示意图；
图4为本发明的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法的CUDA内存模式与GPU存储模式示意图；
图5为本发明的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法的实施例1的实验结果图；
图6为本发明的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法的实施例1的加速比对比图。
具体实施方式
为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。
本发明提供一种基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，如图1～6所示，基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法运行于带有若干个GPU的计算机或服务器，基于多松弛格子玻尔兹曼模型的地下水流动模拟加速系统的加速方法，包括：
步骤1：将地下水区域离散化为若干个网格，每个网格格点给定其初始宏观物理量以及初始分布函数，并将所有网格平均分为若干块传入若干个已开辟好空间的GPU，块内包含多个格点，一个块为单位在每个GPU进行计算；
步骤2：每个GPU调用多个线程同时进行各个块内网格格点的初始宏观物理量及其对应的初始分布函数的初始化；
步骤3：迭代开始，每个GPU调用多个线程同时根据多松弛格子玻尔兹曼方程进行块内网格格点碰撞操作；
步骤4：各个GPU将块内网格格点的分布函数组合到0号GPU上，由0号GPU作为主线程进行格点之间的迁移操作，将之后再次获得的分布函数分发给其他GPU；
步骤5：各个GPU调用多个线程将块内各个网格格点的分布函数值求和，即得到网格非边界格点的宏观物理量；
步骤6：各个GPU分别进行块内网格格点的边界处理，如果块内格点是边界格点，按照上下格点、左右格点个前后格点不同的处理格式进行处理，处理后得到网格边界格点的分布函数和宏观物理量；
步骤7：将步骤5得到的网格非边界格点的宏观物理量和步骤6得到的网格边界格点的宏观物理量均传回CPU，与开始迭代时的对应的网格格点的宏观物理量进行对比判断系统稳定性，若误差大于设定阈值，则继续迭代，执行步骤3，否则，迭代结束，释放所有GPU的内存。
本发明的有益效果在于：
本发明通过多松弛格子玻尔兹曼模型模拟了地下水流动过程并将其在具有多块GPU卡作为加速部件的服务器实现了并行，在占地更小、成本更低、能耗更小的情况下获得更大的加速效果，在保证精度的基础上具有较高的运行效率。
进一步的，步骤1中，网格格点的初始分布函数f的数量由不同的DnQb模型来决定，n表示维度，b表示不同方向分布函数f的数量。
进一步的，步骤2包括：
给定地下水流动偏微分方程：
中渗透系数Kxx、Kyy和Kzz，单位体积流量W，孔隙介质的贮水率Ss等宏观物理量的初始值，其中，公式(1)中：Kxx、Kyy和Kzz分别为x、y和z方向上的渗透系数分量；h为所求宏观物理量，表示水头；W为单位体积流量，用来表示流入或流出流体区域的水量；Ss为孔隙介质的贮水率；t为时间；
根据公式(2)、(3)初始化每个网格格点的水头并根据该水头得出分布函数值，一般初始分布函数设为其平衡态分布函数，
其中，公式(2)中，ρ0、u0和T0是初始时刻的密度、速度、温度的宏观物理量，在地下水流动模拟中，宏观物理量为各个网格格点的水头h，公式(3)中，ρ同公式(2)中的ρ0，表示该格点的宏观物理量，可以是密度、温度等，在地下水中表示水头h，ωi表示粒子向各个方向的权重，u表示该格点的宏观速度，ci表示格点向各个方向的离散速度，cs表示声速。
优选的，步骤3中，多松弛格子玻尔兹曼方程表示为
fi(x+ciδt，t+δt)-fi(x，t)＝-Λij[fi(x，t)-fi (eq)(x，t)] (4)
其中δt表示时间步长，Λ是一个b×b的矩阵，表示碰撞因子，经过一系列的变化，公式(4)改写为
其中m和m(eq)分别为分布函数f和平衡态分布函数f(eq)的矩空间，即m＝Mf，m(eq)＝Mf(eq)，S为松弛对角矩阵，S＝MΛM-1。
优选的，步骤4中，所述迁移操作根据DnQb模型中分布函数的方向进行，即
fi(x+ciδt，t+δt)＝fi′(x，t) (6)
其中ci表示离散速度，对应DnQb模型中的方向坐标，如图2所示，常见的DnQb模型包括D1Q2、D1Q3、D2Q4、D2Q5、D2Q9以及D3Q19等。
本发明中，步骤5中，根据DnQb模型，每个网格格点具有b个不同方向的分布函数，所以分布函数的和即为该网格非边界格点的宏观物理量，即
∑ifi＝ρ (7)
其中ρ表示任一宏观物理量，ρ表示为水头h。
作为本发明的一种改进，步骤6中，边界条件包括定水头边界、变水头边界以及无渗流边界，对于定水头边界和变水头边界，采用非平衡态外推格式来处理，即
fi(xb，t)＝fi (eq)(xb，t)+fi (ne)(xb，t) (8)
其中fi (ne)表示非平衡态部分，对于边界处有一定速度的边界来说，边界格点xb处的宏观uω速度和其相邻格点xf的宏观量密度ρf都是已知的，代入公式(8)中，得到速度边界的非平衡态外推格式，即
fi(xb，t)＝fi (eq)(ρ(xb，t)，uω)+[fi(xf，t)-fi (eq)(xf，t)] (9)
对于无渗流边界，将边界处格点的分布函数设为临近点对应位置网格边界格点的分布函数，即
fi(xb，t)＝fi(xf，t) (10)。
本发明的步骤7中，通过当前宏观物理量与上一时刻宏观物理量误差的最大值，判断系统是否达到稳定，当误差小于设定阈值时，即
则认为系统稳定，否则继续迭代，直到特定时刻或者系统达到稳定状态。
本发明的步骤3中，利用部分CUBLAS库以及CUDA的kernel函数结合OpenMP技术，实现其在多个GPU上的线程级并行，CUBLAS库包括矩阵、向量等操作子函数，例如cublasDcopy()、cublasDaxpy()以及cublasDscal()等。
作为本发明的又一改进，基于CUDA架构，通过设定线程块block和块内线程thread的数量，设定并行使用的GPU总线程数量，实现线程级并行；
线程块或线程之间在一定的内存中读取数据，内存包括线程本地内存和寄存、线程块内的共享内存以及块间网格内的全局内存、常量内存和纹理内存；
通过OpenMP基于共享内存的并行架构机制，使多个GPU同步执行，扩充了总线程的数量，从而实现了更多格点的并行，进一步加快了计算速度。
实施例1：
下面通过一个实施例对本发明进行进一步描述，如图3所示，在地下水流动过程中，左右边界的水头始终为30和10，其他边界无渗流情况，对此区域的具体模拟流动的加速方法如图1所示。
步骤1：将Lx×Ly×Lz的三维地下水区域离散化为Nx×Ny×Nz个网格，分别表示行数、列数和层数，其中Nx和Ny比例为7:6，分别为224和192，Nz设置为4，表示竖直方向一共有四层，根据D3Q19模型，每个网格格点具有格子的水头以及19个方向的分布函数，所有网格平均分为若干块传入若干个已开辟好空间的GPU中；
步骤2：每个GPU调用多个线程同时进行各个块内格点宏观物理量以及19个分布函数的初始化：
对于求解的地下水流动偏微分方程：
其中：Kxx、Kyy和Kzz分别为x、y和z方向上的渗透系数分量；h为所求宏观物理量，表示水头；W为单位体积流量，用来表示流入(汇)或流出(源)流体区域的水量；Ss为孔隙介质的贮水率；t为时间。我们初始化K＝Kxx＝Kyy＝Kzz＝0.5，W＝0，Ss＝0.00001。这时，地下水流动偏微分方程简化为
其中，扩散系数然后对各个格点的水头进行初始化，除右边界处格点的水头设为10以外，其他所有格点的水头设为30；
根据公式(2)、公式(3)，可以通过水头值计算出每个格点的19个方向的分布函数。
fi＝fi (eq)(ρ0，u0，T0) (2)
其中，由于格点在此实施例中宏观速度为0，原平衡态分布函数公式右端括号里的内容简化为1，故公式(3)简化为公式(3-1)
f(eq)＝ωiρ (3-1)
对于D3Q19模型，离散速度和权重分别表示为
步骤3：迭代开始，每个GPU调用多个线程同时进行块内格点根据多松弛格子玻尔兹曼方程进行碰撞操作：
多松弛格子玻尔兹曼方程可表示为
fi(x+ciδt，t+δt)-fi(x，t)＝-Λij[fi(x，t)-fi (eq)(x，t)] (4)
其中Λ是一个b×b的矩阵，表示碰撞因子，经过一系列变化，公式(4)可以改写为
其中，m和m(eq)分别分布函数f和平衡态分布函数f(eq)的矩空间，即m＝Mf，m(eq)＝Mf(eq)，其中M为转换矩阵，对于D3Q19模型，矩阵M可表示为
S为松弛矩阵，S＝MΛM-1，对于D3Q19模型，它是一个19*19的对角矩阵，可表示为
S＝diag(0，se，ss，0，sq，0，sq，0，sq，sv，sπ，sv，sπ，sv，sv，sv，st，st，st)
其中，黏性系数为
步骤4：各个GPU将块内格点分布函数组合到0号GPU上，有0号GPU作为主线程进行格点之间的迁移操作，之后再将获得的分布函数分发给其他几个GPU：
迁移时根据DnQb模型中分布函数的方向进行，即
fi(x+ciδt，t+δt)＝f′i(x，t) (6)
其中ci表示离散速度，对于D3Q19模型，i表示0到18。
步骤5：各个GPU调用多个线程将块内各个格点的分布函数值求和，即得到该格点的宏观物理量。迁移后每个格点具有b个不同方向的分布函数，所有分布函数的和即为该格点的宏观物理量，即
∑ifi＝ρ (7)
其中ρ表示任一宏观物理量，根据地下水流动偏微分方程，ρ可以表示为水头h。对于D3Q19模型，当前格点宏观量水头的值即为该格点19个方向的分布函数的和。
步骤6：各个GPU分别进行块内格点的边界处理，如果块内格点是边界格点，按照上下格点、左右格点和前后格点不同的处理格式进行处理，处理后才可以得到边界格点未知的分布函数和宏观物理量：
本实例中我们设定左右边界值始终为30和10，其余边界包括上下边界和前后边界设为无渗流情况，如图3所示。左右边界的处理使用非平衡态外推格式，即
其他边界使用无渗流边界条件，即
fi(xb，t)＝fi(xf，t) (10)
步骤7：将各个GPU的块内格点宏观物理量传回CPU，与开始迭代时各个格点宏观物理量进行对比判断系统稳定性，如误差大于一个规定的范围则需要继续进行下一步迭代，否则迭代结束，释放所有GPU的内存。
通过当前宏观物理量与上一时刻宏观物理量误差的最大值，判断系统是否达到稳定，当误差小于我们要求的某一特定值时，即
我们认为系统稳定，否则继续迭代，直到特定时刻或者系统达到稳定状态。本实例中，error的数值我们设为1e-5，本发明通过多松弛格子玻尔兹曼模型模拟了地下水流动过程，结果如图5所示。
基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法运行于带有多块GPU卡作为加速部件的服务器，具体的，该发明利用部分CUBLAS库以及CUDA的kernel函数结合OpenMP技术，实现其在多个GPU上的线程级并行，CUBLAS库包括矩阵、向量等操作子函数，例如cublasDcopy()实现向量的复制，cublasDaxpy()实现向量的乘积和，cublasDscal()实现向量积。
通过CUDA实现函数的线程级并行时，需要设定线程块block和块内线程thread的数量，继而设定并行使用的GPU总线程数量，实现线程级并行。我们将线程设为16，线程块设为二维(Nx*Nz/numOfThreads,Ny)，在CUDA调用的设备端kernel函数中，通过Nx*Nz*blockIdx.y+threadIdx.x+blockIdx.x*blockDim.x获得绝对线程序号，从而控制每一个线程执行相应的操作，CUDA的执行模型与GPU的存储模型如图4所示。
线程块或线程之间在一定的内存中读取数据，内存包括线程自己的本地内存和寄存，线程块内的共享内存，以及块间网格内的全局内存、常量内存和纹理内存。由于同一线程块中的线程可以对该线程块的共享内存中的数据进行读取，且读写共享内存的速度很快，所以我们将一些常用的常量(例如离散速度、权重等)通过__shared__定义的方式存储在每一个块的共享内存中，这样在多次调用这些常向量的时候，可以节省很多的时间；对于迁移操作而言，所有格点的分布函数值我们存储在全局内存global memory中，这样即使格点划分进不同的线程块，再由各个线程块中的每个线程同时执行，也保证了每个线程可以从全局内存中找到不同块内临近格点的分布函数值。
进一步的，通过OpenMP基于共享内存的并行架构机制，使多个GPU同步执行，扩充了总线程的数量，从而实现了更多格点的并行，进一步加快了计算速度。
由于每个GPU只进行部分网格的计算，所以对于Nx×Ny×Nz个网格来说，每个GPU实际上只计算了Nx×Ny×Nz/num_gpu个网格，本文实验使用的服务器中GPU数量为6个，为了划分数据的均衡，我们分别进行了2个GPU、4个GPU以及6个GPU的实验。
迁移操作时，每个GPU中格点的分布函数通过cudaMemcpyPeer()传入0号GPU，并在0号GPU上进行组合，然后通过#pragma omp single语句让0号GPU独自完成迁移操作，之后再将计算之后的格点分布函数再通过cudaMemcpyPeer()传回其他各个GPU，继续进行宏观量的计算和边界处理。
本实施例1在具有多块GPU卡作为加速部件的服务器上实现并行获得的加速比如图6所示，采用多个GPU要比一个GPU加速比大的多，加速明显。
综上，本发明具有以下有益效果：
本发明通过多松弛格子玻尔兹曼模型模拟了地下水流动过程并将其在具有多块GPU卡作为加速部件的服务器实现了并行，在占地更小、成本更低、能耗更小的情况下获得更大的加速效果，在保证精度的基础上具有较高的运行效率。
以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (10)
1.一种基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述方法运行于带有若干个GPU的计算机或服务器，所述方法，包括：
步骤1：将地下水区域离散化为若干个网格，每个网格格点给定其初始宏观物理量以及初始分布函数，并将所有网格平均分为若干块传入若干个已开辟好空间的GPU；
步骤2：每个GPU调用多个线程同时进行各个块内网格格点的初始宏观物理量及其对应的初始分布函数的初始化；
步骤3：迭代开始，每个GPU调用多个线程同时根据多松弛格子玻尔兹曼方程进行块内网格格点碰撞操作；
步骤4：各个GPU将块内网格格点的分布函数组合到0号GPU上，由0号GPU作为主线程进行格点之间的迁移操作，将之后再次获得的分布函数分发给其他GPU；
步骤5：各个GPU调用多个线程将块内各个网格格点的分布函数值求和，即得到网格非边界格点的宏观物理量；
步骤6：各个GPU分别进行块内网格格点的边界处理，得到网格边界格点的分布函数和宏观物理量；
步骤7：将步骤5得到的网格非边界格点的宏观物理量和步骤6得到的网格边界格点的宏观物理量均传回CPU，与开始迭代时的对应的网格格点的宏观物理量进行对比判断稳定性，若误差大于设定阈值，则继续迭代，执行步骤3，否则，迭代结束，释放所有GPU的内存。
2.根据权利要求1所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤1中，网格格点的初始分布函数f的数量由不同的DnQb模型来决定，n表示维度，b表示不同方向分布函数f的数量。
3.根据权利要求1所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤2包括：
给定地下水流动偏微分方程：
<mrow>
<mfrac>
<mo>&part;</mo>
<mrow>
<mo>&part;</mo>
<mi>x</mi>
</mrow>
</mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>K</mi>
<mrow>
<mi>x</mi>
<mi>x</mi>
</mrow>
</msub>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>x</mi>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mo>&part;</mo>
<mrow>
<mo>&part;</mo>
<mi>y</mi>
</mrow>
</mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>K</mi>
<mrow>
<mi>y</mi>
<mi>y</mi>
</mrow>
</msub>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>y</mi>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mo>&part;</mo>
<mrow>
<mo>&part;</mo>
<mi>z</mi>
</mrow>
</mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>K</mi>
<mrow>
<mi>z</mi>
<mi>z</mi>
</mrow>
</msub>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>z</mi>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>W</mi>
<mo>=</mo>
<msub>
<mi>S</mi>
<mi>s</mi>
</msub>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
中渗透系数Kxx、Kyy和Kzz，单位体积流量W，孔隙介质的贮水率Ss的宏观物理量的初始值，其中，公式(1)中：Kxx、Kyy和Kzz分别为x、y和z方向上的渗透系数分量；h为所求宏观物理量，表示水头；W为单位体积流量，用来表示流入或流出流体区域的水量；Ss为孔隙介质的贮水率；t为时间；
根据公式(2)、(3)初始化每个网格格点的水头并根据该水头得出分布函数值，一般初始分布函数设为其平衡态分布函数，
fi＝fi (eq)(ρ0，u0，T0) (2)
<mrow>
<msup>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>e</mi>
<mi>q</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>=</mo>
<msub>
<mi>&omega;</mi>
<mi>i</mi>
</msub>
<mi>&rho;</mi>
<mo>&lsqb;</mo>
<mn>1</mn>
<mo>+</mo>
<mfrac>
<mrow>
<msub>
<mi>c</mi>
<mi>i</mi>
</msub>
<mo>&CenterDot;</mo>
<mi>u</mi>
</mrow>
<msubsup>
<mi>c</mi>
<mi>s</mi>
<mn>2</mn>
</msubsup>
</mfrac>
<mo>+</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>c</mi>
<mi>i</mi>
</msub>
<mo>&CenterDot;</mo>
<mi>u</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>c</mi>
<mi>s</mi>
<mn>4</mn>
</msubsup>
</mrow>
</mfrac>
<mo>-</mo>
<mfrac>
<mrow>
<mi>u</mi>
<mo>&CenterDot;</mo>
<mi>u</mi>
</mrow>
<mrow>
<mn>2</mn>
<msubsup>
<mi>c</mi>
<mi>s</mi>
<mn>2</mn>
</msubsup>
</mrow>
</mfrac>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，公式(2)中，ρ0、u0和T0是初始时刻的密度、速度、温度的宏观物理量，在地下水流动模拟中，宏观物理量为各个网格格点的水头h。
4.根据权利要求1所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤3中，多松弛格子玻尔兹曼方程表示为
fi(x+ciδt，t+δt)-fi(x，t)＝-Λij[fi(x，t)-fi(eq)(x，t)] (4)
其中Λ是一个b×b的矩阵，表示碰撞因子，经过一段时间t的变化，公式(4)改写为
<mrow>
<msub>
<mi>m</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>+</mo>
<msub>
<mi>c</mi>
<mi>i</mi>
</msub>
<msub>
<mi>&delta;</mi>
<mi>t</mi>
</msub>
<mo>,</mo>
<mi>t</mi>
<mo>+</mo>
<msub>
<mi>&delta;</mi>
<mi>t</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>m</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>-</mo>
<mi>S</mi>
<mo>&lsqb;</mo>
<msub>
<mi>m</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msubsup>
<mi>m</mi>
<mi>i</mi>
<mrow>
<mo>(</mo>
<mi>e</mi>
<mi>q</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
</mrow>
其中m和m(eq)分别为分布函数f和平衡态分布函数f(eq)的矩空间，即m＝Mf，m(eq)＝Mf(eq)，S为松弛对角矩阵，S＝MΛM-1。
5.根据权利要求2所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤4中，所述迁移操作根据DnQb模型中分布函数的方向进行，即
fi(x+ciδt，t+δt)＝fi′(x，t) (6)
其中ci表示离散速度，对应DnQb模型中的方向坐标。
6.根据权利要求5所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤5中，根据DnQb模型，每个网格格点具有b个不同方向的分布函数，所以分布函数的和即为该网格非边界格点的宏观物理量，即
∑ifi＝ρ (7)
其中ρ表示任一宏观物理量，ρ表示为水头h。
7.根据权利要求1所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤6中，边界条件包括定水头边界、变水头边界以及无渗流边界，对于定水头边界和变水头边界，采用非平衡态外推格式来处理，即
fi(xb，t)＝fi (eq)(xb，t)+fi (ne)(xb，t) (8)
其中fi (ne)表示非平衡态部分，对于边界处有一定速度的边界来说，边界格点xb处的宏观uω速度和其相邻格点xf的宏观量密度ρf都是已知的，代入公式(8)中，得到速度边界的非平衡态外推格式，即
fi(xb，t)＝fi (eq)(ρ(xb，t)，uω)+[fi(xf，t)-fi (eq)(xf，t)] (9)
对于无渗流边界，将边界处格点的分布函数设为临近点对应位置网格边界格点的分布函数，即
fi(xb，t)＝fi(xf，t) (10)。
8.根据权利要求1所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤7中，通过当前宏观物理量与上一时刻宏观物理量误差的最大值，判断系统是否达到稳定，当误差小于设定阈值时，即
<mrow>
<mi>max</mi>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mrow>
<msubsup>
<mi>&rho;</mi>
<mi>i</mi>
<mn>0</mn>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&rho;</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mn>0</mn>
</msubsup>
</mrow>
<mo>|</mo>
<mo>,</mo>
<mo>|</mo>
<mrow>
<msubsup>
<mi>&rho;</mi>
<mi>i</mi>
<mn>1</mn>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&rho;</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mn>1</mn>
</msubsup>
</mrow>
<mo>|</mo>
<mo>,</mo>
<mn>...</mn>
<mo>|</mo>
<mrow>
<msubsup>
<mi>&rho;</mi>
<mi>i</mi>
<mi>n</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&rho;</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</msubsup>
</mrow>
<mo>|</mo>
<mo>)</mo>
</mrow>
<mo><</mo>
<mi>e</mi>
<mi>r</mi>
<mi>r</mi>
<mi>o</mi>
<mi>r</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>11</mn>
<mo>)</mo>
</mrow>
</mrow>
则认为系统稳定，否则继续迭代，直到特定时刻或者系统达到稳定状态。
9.根据权利要求1所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，所述步骤3中，利用部分CUBLAS库以及CUDA的kernel函数结合OpenMP技术，实现其在多个GPU上的线程级并行。
10.根据权利要求9所述的基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法，其特征在于，基于所述CUDA，通过设定线程块block和块内线程thread的数量，设定并行使用的GPU总线程数量，实现线程级并行；
线程块或线程之间在一定的内存中读取数据，内存包括线程本地内存和寄存、线程块内的共享内存以及块间网格内的全局内存、常量内存和纹理内存。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710741065.5A CN107515987A (zh) | 2017-08-25 | 2017-08-25 | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710741065.5A CN107515987A (zh) | 2017-08-25 | 2017-08-25 | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107515987A true CN107515987A (zh) | 2017-12-26 |
Family
ID=60724048
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710741065.5A Pending CN107515987A (zh) | 2017-08-25 | 2017-08-25 | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107515987A (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108460195A (zh) * | 2018-02-08 | 2018-08-28 | 国家海洋环境预报中心 | 海啸数值计算模型基于gpu并行的快速执行方法 |
| CN109388876A (zh) * | 2018-09-28 | 2019-02-26 | 中国地质大学（北京） | 一种地下水溶质运移数值模拟并行加速方法 |
| CN109857543A (zh) * | 2018-12-21 | 2019-06-07 | 中国地质大学（北京） | 一种基于多节点多gpu计算的流线模拟加速方法 |
| CN111680456A (zh) * | 2020-04-28 | 2020-09-18 | 中国科学院深圳先进技术研究院 | 一种流体力学模拟的方法、装置及存储介质 |
| CN112836872A (zh) * | 2021-01-29 | 2021-05-25 | 西安理工大学 | 一种基于多gpu的污染物对流扩散方程高性能数值求解方法 |
| CN113791013A (zh) * | 2021-08-31 | 2021-12-14 | 长江岩土工程有限公司 | 病险水库大坝土质防渗体钻孔连续注水试验方法 |
| CN114155135A (zh) * | 2021-12-02 | 2022-03-08 | 西北核技术研究所 | 基于gpu集群的相对论性波尔兹曼方程计算方法、存储介质及设备 |
| CN115345082A (zh) * | 2022-07-06 | 2022-11-15 | 中山大学 | 一种面向冲击系统的二维九速离散玻尔兹曼方法及装置 |
| CN115495919A (zh) * | 2022-09-30 | 2022-12-20 | 中国船舶科学研究中心 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
| CN116306382A (zh) * | 2023-05-17 | 2023-06-23 | 湖南大学 | 一种非饱和水泥浆体中水分与气体分布的数值模拟方法 |
| CN116842691A (zh) * | 2023-05-24 | 2023-10-03 | 中国水利水电科学研究院 | 一种智能提高地下水数值模拟收敛性的松弛方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102945295A (zh) * | 2012-10-15 | 2013-02-27 | 浪潮(北京)电子信息产业有限公司 | 一种格子玻尔兹曼方法的并行加速方法及系统 |
| CN103324531A (zh) * | 2013-06-09 | 2013-09-25 | 浪潮电子信息产业股份有限公司 | 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 |
| CN103345580A (zh) * | 2013-07-02 | 2013-10-09 | 上海大学 | 基于格子Boltzmann方法的并行CFD方法 |
-
2017
- 2017-08-25 CN CN201710741065.5A patent/CN107515987A/zh active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102945295A (zh) * | 2012-10-15 | 2013-02-27 | 浪潮(北京)电子信息产业有限公司 | 一种格子玻尔兹曼方法的并行加速方法及系统 |
| CN103324531A (zh) * | 2013-06-09 | 2013-09-25 | 浪潮电子信息产业股份有限公司 | 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 |
| CN103345580A (zh) * | 2013-07-02 | 2013-10-09 | 上海大学 | 基于格子Boltzmann方法的并行CFD方法 |
Non-Patent Citations (2)
| Title |
|---|
| LIANKAI YAO 等: "Parallelism of MRT Lattice Boltzmann Method based on Multi-GPUs", 《CONFERENCE: THE 7TH INTERNATIONAL CONFERENCE ON COMPUTER ENGINEERING AND NETWORKS》 * |
| 郭照立 等: "《格子Boltzmann方法的原理及应用》", 31 January 2009 * |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108460195A (zh) * | 2018-02-08 | 2018-08-28 | 国家海洋环境预报中心 | 海啸数值计算模型基于gpu并行的快速执行方法 |
| CN109388876A (zh) * | 2018-09-28 | 2019-02-26 | 中国地质大学（北京） | 一种地下水溶质运移数值模拟并行加速方法 |
| CN109857543A (zh) * | 2018-12-21 | 2019-06-07 | 中国地质大学（北京） | 一种基于多节点多gpu计算的流线模拟加速方法 |
| CN111680456A (zh) * | 2020-04-28 | 2020-09-18 | 中国科学院深圳先进技术研究院 | 一种流体力学模拟的方法、装置及存储介质 |
| CN112836872B (zh) * | 2021-01-29 | 2023-08-18 | 西安理工大学 | 一种基于多gpu的污染物对流扩散方程高性能数值求解方法 |
| CN112836872A (zh) * | 2021-01-29 | 2021-05-25 | 西安理工大学 | 一种基于多gpu的污染物对流扩散方程高性能数值求解方法 |
| CN113791013A (zh) * | 2021-08-31 | 2021-12-14 | 长江岩土工程有限公司 | 病险水库大坝土质防渗体钻孔连续注水试验方法 |
| CN114155135A (zh) * | 2021-12-02 | 2022-03-08 | 西北核技术研究所 | 基于gpu集群的相对论性波尔兹曼方程计算方法、存储介质及设备 |
| CN115345082A (zh) * | 2022-07-06 | 2022-11-15 | 中山大学 | 一种面向冲击系统的二维九速离散玻尔兹曼方法及装置 |
| CN115495919A (zh) * | 2022-09-30 | 2022-12-20 | 中国船舶科学研究中心 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
| CN115495919B (zh) * | 2022-09-30 | 2023-05-26 | 中国船舶科学研究中心 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
| CN116306382A (zh) * | 2023-05-17 | 2023-06-23 | 湖南大学 | 一种非饱和水泥浆体中水分与气体分布的数值模拟方法 |
| CN116306382B (zh) * | 2023-05-17 | 2023-08-08 | 湖南大学 | 一种非饱和水泥浆体中水分与气体分布的数值模拟方法 |
| CN116842691A (zh) * | 2023-05-24 | 2023-10-03 | 中国水利水电科学研究院 | 一种智能提高地下水数值模拟收敛性的松弛方法 |
| CN116842691B (zh) * | 2023-05-24 | 2024-03-08 | 中国水利水电科学研究院 | 一种智能提高地下水数值模拟收敛性的松弛方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107515987A (zh) | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 | |
| Müller et al. | Real time dynamic fracture with volumetric approximate convex decompositions | |
| Stam | Real-time fluid dynamics for games | |
| Odstrcil et al. | Initial coupling of coronal and heliospheric numerical magnetohydrodynamic codes | |
| CN104136942B (zh) | 用于大规模并行储层仿真的千兆单元的线性求解方法和装置 | |
| Ullrich et al. | High-order finite-volume methods for the shallow-water equations on the sphere | |
| Loubere et al. | ReALE: a reconnection-based arbitrary-Lagrangian–Eulerian method | |
| DeVries et al. | A first course in computational physics | |
| CN107590853A (zh) | 一种城市建筑群震害高真实度展示方法 | |
| US8144155B2 (en) | Example-based motion detail enrichment in real-time | |
| Li et al. | GPU-based numerical simulation of multi-phase flow in porous media using multiple-relaxation-time lattice Boltzmann method | |
| Mallison et al. | Practical gridding algorithms for discrete fracture modeling workflows | |
| Castillo et al. | Linear systems arising for second-order mimetic divergence and gradient discretizations | |
| Thürey et al. | Interactive free surface fluids with the lattice Boltzmann method | |
| Rybakin | Modeling of III-D problems of gas dynamics on multiprocessing computers and GPU | |
| CN106096132A (zh) | 一种基于微分域的不同材料衣物褶皱的仿真方法 | |
| Loubère et al. | A totally Eulerian Finite Volume solver for multi-material fluid flows: Enhanced Natural Interface Positioning (ENIP) | |
| Jiang et al. | MHD modeling of solar coronal magnetic evolution driven by photospheric flow | |
| Zhou et al. | Implicit integration for particle‐based simulation of elasto‐plastic solids | |
| Li et al. | Interactive cutting and tearing in projective dynamics with progressive cholesky updates | |
| Wilderotter | An adaptive numerical method for the Richards equation with root growth | |
| CN109829232A (zh) | 基于随机森林算法的分层布料模拟方法 | |
| Akhmed-Zaki et al. | Large-scale simulation of oil recovery by surfactant-polymer flooding | |
| CN108427605B (zh) | 基于质点追踪算法实现流线模拟的加速方法 | |
| CN109726496A (zh) | 一种基于iisph提高不可压缩水体模拟效率的实现方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20171226 |
|
| RJ01 | Rejection of invention patent application after publication |