CN103324531A - 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 - Google Patents
一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 Download PDFInfo
- Publication number
- CN103324531A CN103324531A CN2013102291613A CN201310229161A CN103324531A CN 103324531 A CN103324531 A CN 103324531A CN 2013102291613 A CN2013102291613 A CN 2013102291613A CN 201310229161 A CN201310229161 A CN 201310229161A CN 103324531 A CN103324531 A CN 103324531A
- Authority
- CN
- China
- Prior art keywords
- mic
- thread
- cpu
- num
- grid
- 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
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法,包括采用LBM方法的大涡模拟方法、CPU端、MIC众核协处理器端以及CPU+MIC协同计算模式,其中:CPU端负责将要进行大涡模拟的网格数据分割,向MIC卡传递大涡模拟所需要的值,CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作,而且在整个网格的计算任务中,CPU也会以Openmp多线程模式,依次通过迁移碰撞,边界处理的多次迭代获取得速度、密度和流函数的宏观参量;MIC众核协处理器负责网格点的迁移和碰撞过程,对边界进行处理,根据分布函数并行求得速度、密度和流函数等宏观参量,在MIC卡上也采用openmp多线程的方式来运算。
Description
技术领域
本发明涉及计算机应用技术领域,具体地说是一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法。
背景技术
MIC(Many Integrated Core)是Intel公司推出的众核处理器,跟通用的多核至强处理器相比,MIC众核架构具有更小的内核和硬件线程,众核处理器计算资源密度更高,片上通信开销显著降低,更多的晶体管和能量,能够胜任更为复杂的并行应用。Intel MIC产品基于X86架构,基于重核的众核处理器,包含50个以上的核心,以及512bit的向量位宽,双精性能超过1TFlops。
MIC拥有极其灵活的编程方式,MIC卡可以作为一个协处理器存在,也可以被看作是一个独立的节点。基本的MIC编程模型是将MIC看作一个协处理器,CPU根据程序的指令,将一部分代码运行在MIC端。此时存在两类设备,即CPU端和MIC众核协处理器端。
大涡模拟,英文简称LES(Large eddy simulation),是近几十年才发展起来的一个流体力学中重要的数值模拟研究方法。它区别于直接数值模拟(DNS)和雷诺平均(RANS)方法。其基本思想是通过精确求解某个尺度以上所有湍流尺度的运动,从而能够捕捉到RANS方法所无能为力的许多非稳态、非平衡过程中出现的大尺度效应和拟序结构,同时又克服了直接数值模拟需要求解所有湍流尺度而带来的巨大计算开销的问题,因而被认为是最具有潜力的湍流数值模拟发展方向。由于计算耗费依然很大,目前大涡模拟还无法在工程上广泛应用,但是大涡模拟技术对于研究许多流动机理问题提供了更为可靠的手段,可为流动控制提供理论基础,并可为工程上广泛应用的RANS方法改进提供指导。LES算法过程如图1所示。
LBM是计算流体力学领域内一种不同于传统数值方法的建模和计算方法,是对Boltzmann方程的一种特殊的离散格式的求解。求解过程是时间推进式的,并且求解过程具有良好的区域性,所以特别适合并行求解。
一般对Lattice Boltzmann方程求解,可以分解为两部分:
1)碰撞项:
2)迁移项:
图2所示的LBM的D2Q9模型中,平衡分布函数具体定义下:
数得到,具体计算式为:
(5)
边界条件的处理方法在LBM 的应用中对数值精度和计算的稳定性都有着极大的影响。在本算法中,采用非平衡外推边界条件方法来处理固壁边界条件。其基本思想是将固壁边界格点上的未知分布函数分解为平衡态和非平衡态两部分,然后利用一阶精度的外推格式得到非平衡部分。该格式可以表示为:
发明内容
本发明的目的是提供一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法。
本发明所要解决的技术问题是提供一种基于格子Boltzmann理论的大涡模拟的高效方法,能够基于MIC众核架构并且使用CPU+MIC协同计算的快速大涡模拟。
包括采用LBM方法的大涡模拟方法、CPU端、MIC众核协处理器端以及CPU+MIC协同计算模式,其中:
CPU端负责将要进行大涡模拟的网格数据分割,向MIC卡传递大涡模拟所需要的值,CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作,而且在整个网格的计算任务中,CPU也会以Openmp多线程模式,依次通过迁移碰撞,边界处理的多次迭代获取得速度、密度和流函数的宏观参量;
MIC众核协处理器负责网格点的迁移和碰撞过程,对边界进行处理,根据分布函数并行求得速度、密度和流函数等宏观参量。在MIC卡上也采用openmp多线程的方式来运算;
进一步地,CPU端将要进行大涡模拟的网格数据分割,向MIC卡传递大涡模拟所需要的值,以openmp多线程模式,并行执行迁移碰撞,边界处理来获取速度、密度和流函数的宏观参量,具体包括:
LES_MIC算法中,根据网格的划分,在对划分边界要做一些特殊的处理,在迁移过程中,每个点的计算需要其周围9个方向的分布函数值,因此,需要在分布函数中多存储1或2行用于下一步时层的计算,并且在每次迭代之后在节点与节点之间、MIC卡之间进行边界数据的交换;下表表示各MIC卡节点上数组大小;
进一步地,CPU负责CPU+MIC协同计算模式的框架搭建以及任务调度,具体包括:
单节点服务器采用由双路6核CPU和2块KNF MIC卡组成的桌面服务器,在CPU+MICs协同计算中,把双路CPU和MIC卡都作为计算设备,这样每个单节点就相当于有3个计算设备,每个设备通过一个OpenMP线程进行控制;
本方法是数据级并行,所以采用静态数据划分的方式,每次各个设备读取设定好的网格数据以及边界处理需要的数据,然后分别进行数据处理,相邻设备需要交换数据,迭代多次,直到所有设备计算完成所有网格数据,由CPU输出结果;
进一步地,MIC众核协处理器负责网格点的迁移和碰撞过程,对边界进行处理,根据分布函数并行求得速度、密度和流函数的宏观参量,在MIC卡上也采用openmp多线程的方式来运算;
根据对LBM算法中求解离散方程和边界处理的串行算法的热点分析和并行性分析,每个网格点的迁移、碰撞、宏观量统计、平衡态分布函数计算和边界处理的计算都是数据并行的;
进一步地,求解离散方程可以采用迁移碰撞的过程,宏观量统计、平衡态分布函数计算和碰撞过程中对每个网格的计算之间没有任何依赖性,因此,让MIC中的每个线程负责一个网格划分中的一行网格点的计算,每行网格点的计算利用MIC上的向量化技术进一步加速;分布函数的迁移只涉及到该格点周围的其他格点,或通过单个线程对全局存储器中相关分布函数的读操作来实现;
进一步地,在LBM算法中,对边界要做特殊的处理,包括:非平衡外推,反弹,对于边界上的每个格点之间的计算也没有数据的依赖性,因此,利用OpenMP多线程负责边界格点的计算;
再进一步地,OpenMP的线程模型设计:根据MIC核心数目设置内核的线程数。
包括如下步骤:
步骤一:CPU+MIC协同计算框架搭建
在单节点上,共有M+1个设备(一个CPU+M个MIC卡),采用OpenMP的fork-join模式搭建单节点上的框架,当程序开始执行的时候只有一个主线程存在,需要进行并行计算的时候,主线程派生出附加线程,即启用M+1个OpenMP线程, 0~M-1号进程控制所有MIC设备,M号线程控制CPU设备,按照数据的分配设计,每个设备进行输入数据的分发读写,
如在本专利中,CPU+2MICs平台上, 主线程控制输入数据的动态分发,0号线程控制MIC0设备,1号线程控制MIC1设备,2号线程控制CPU设备,如图3所示。
在CPU和MIC卡上,各个设备的数据大小的划分如表2所示,对于分布函数,其大小要比网格大小大两行,同时方便代码的书写,我们对最上面的设备和最下面的设备也都做了+2行的数据定义(有一行在代码中不需要使用)。
单节点上CPU+MICs协同计算的外围框架的伪代码如下所示:
定义一些变量
//在迭代范围之内
int DEVICE_NUM //设备数
for(int i=0;i<steps;i++)
{
if(i%2==0) //奇数步和偶数步输入和输出交换,所以有一个迭代的判断
{
omp_set_nested(true);
#pragma omp parallel for private(…), num_threads(DEVICE_NUM)
for(int thread=0;thread<DEVICE_NUM;thread++)
{
if(thread==0) //mic 0 computing
{
#pragma offload target(mic:0) \
in(fs0_in0_up,…:length(nx) alloc_if(0) free_if(0))\
out(fn1_out0_up,:length(nx) alloc_if(0) free_if(0))\
nocopy(fr0_mic0,…:length((hh+1)*nx) alloc_if(0) free_if(0)) \
nocopy(fr1_mic0,…:length((hh+1)*nx) alloc_if(0) free_if(0))
{
……
LBCollProp(DEVICE_NUM, thread_num_mic,thread,… );
LBBC_LR(DEVICE_NUM, thread_num_mic,thread, …);
LBBC_DOWN(DEVICE_NUM, thread_num_mic,thread , …);
……
}//mic0 end
}
else if(thread==DEVICE_NUM-1) //cpu computing
{
……
LBCollProp(DEVICE_NUM, thread_num_omp,thread,…);
LBBC_LR(DEVICE_NUM, thread_num_omp,thread,…);
LBBC_UP(DEVICE_NUM, thread_num_omp,thread,…);
}
else //other mic computing
{
#pragma offload target(mic:1) \
in(thread_num_mic,thread,nx,…)\
in(fs1_in1_up,…:length(nx) alloc_if(0) free_if(0))\
out(fn0_out1_up,…:length(nx) alloc_if(0) free_if(0))\
nocopy(fr0_mic1,…:length((hh+2)*nx) alloc_if(0) free_if(0)) \
…
{//mic1 compute
……
LBCollProp(DEVICE_NUM, thread_num_mic,thread,… );
LBBC_LR(DEVICE_NUM, thread_num_mic,thread, …);
LBBC_DOWN(DEVICE_NUM, thread_num_mic,thread, nx, hh, fr0_mic0, fe0_m …);
……
}
else //奇数步
{
//函数内容和偶数步是一样的,只是输入输出和偶数步互换
}
步骤二:CPU/MIC内核实现
(3) 设计迁移碰撞内核,设计线程数为T=4*M,M为MIC卡的核数,并且让内核中的每个线程计算一行网格点的迁移和碰撞过程,以及利用#pragma ivdep实现内核中内层循环的向量化,如图4所示,内核伪代码如下;
1: #pragma omp parallel for private (i, j, k,…) num_threads(T)//T为线程数
2: for (i=1;i<NY-1;i++)
3: #pragma ivdep //向量化
4: for(j=1;j<NX-1;j++)
5: {
6: k=i*NX+j; //k表示网格的标号
7: fr = fr0[k];//0代表上一时层的分布函数
8: fe = fe0[k-1];
9: fn = fn0[k-NX];
10: fw = fw0[k+1];
11: fs = fs0[k+NX];
12: fne = fne0[k-NX-1];
13: fnw = fnw0[k-NX+1];
14: fsw = fsw0[k+NX+1];
15: fse = fse0[k+NX-1];
16: /*碰撞过程*/
17: 根据迁移后的分布函数fr-fse 求宏观量
18: 根据宏观量求各个方向的平衡分布函数
f1,f2,f3,f4,f5,f6,f7,f8;
19: 根据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];
20: }
(4)在MIC端对边界进行处理,边界处理可以采用反弹法、非平衡外推法等方法,在对边界的处理时同样设计T个线程处理边界节点的计算;
步骤三:数据传递方式设计
节点间的数据传递
LES算法根据格点所在区域将格点划分到不同的设备上,因此当每个格点更新自己的分布函数并迁移时,各个划分区域中边界格点的分布函数要传递给邻近的节点。如图5所示。图5中位于H、L虚线之间的黑实线将计算域一分为二,但为了在计算中方便该实线的迁移的计算,需要邻近网格,M(i+1,j)的计算域增加了该实线下边邻近网格L,M(i,j)的计算域增加了该实线上边邻近网格H。M(i+1,j)的计算域中H的分布函数需要传递给M(i,j),同时也需要接收从M(i,j)传递过来的L的分布函数;对M(i,j)也类似处理。
具体优化1:LES_MIC优化,实施过程如下:
1)向量化
对于MIC内核,我们设计外层for并行,内层for采用向量化的方案,对于各个内核函数,我们可以采取自动向量化的方案进行优化;
2)减少offload次数
在迭代过程中,尽量减少offload次数,以及CPU和MIC之间的I/O次数;
3)减少节点与节点的数据传递以及MIC卡之间的数据传递
在每次迭代之后,相邻节点之间,以及相邻的MIC卡之间要进行边界数据的传递,而每个网格有9个方向的分布函数,然而在内核计算中,并不需要边界上的全部9个方向上的值,
只需要3个方向的值即可,如图6所示。对于节点i,只需要接收节H中的fsw,fs,fse的值,同样,对于节点i+1,只需要接收节L中的fnw,fn,fne的值。
本项目中LES应用案例测试用的软硬件环境如表3所示。
LES性能采用格子点更新速率(LUPS, Lattice Unit Per Second)统计,常为MLUPS(每秒更新的百万网格数),计算方法:
P=NX*NY*S/T
其中NX、NY为网格宽和高,S为流场迭代次数,T为计算时间,P为格子点更新速率。
表3 LES实验环境
平台 | Inspur NF5280M3,2节点 |
CPU | Intel Xeon CPU E5-2680 2.7GHz,双路8核 |
Memory | DDR3 1333MHz 128GB |
MIC | KNC,61核,1.1GHz,GDDR5 8GB memory,5.5GT/s |
OS | Red Hat Enterprise Linux Server release 6.1,64bit |
驱动 | KNC_beta_oem-2.1.3653-8-rhel-6.1 |
编译器 | l_ccompxe_intel64_2013.0.079 |
测试用例 | 雷诺数:10000;迭代次数:10000 |
基于单节点CPU+MIC协同计算LES并行算法的实验结果如图7所示,
由本发明给出的上述测试结果可以看出,在MIC众核架构平台上,基于格子Boltzmann算法和CPU+MIC异构并行架构的基础上,可以在较短周期内,较容易地大幅度加速大涡模拟的运算。
本发明的有益效果是: 本发明不仅提高了格子Boltzmann方法的处理性能,满足了流体模拟的需求,而且能充分利用CPU以及MIC协处理器的降低功耗,减少机房构建成本和管理、运行、维护费用,并且这种方法实现简单,需要的开发成本低。
附图说明
图1是采用迭代重构图像方法实现图像重构的过程图;
图2是CPU+MIC协同计算结构图;
图3是CPU+MIC协同计算程序框架图;
图4是算法代码流程图;
图5是同一节点上MIC卡之间的数据传递图;
图6是数据传递边界部分迁移过程图;
图7是单节点CPU+MIC协同计算实验结果图;
图8 是单节点CPU+MIC协同计算相对串行加速比示意图;
图9 是单节点CPU+MIC协同计算相对单节点OpenMP多线程加速比示意图。
具体实施方式
以下结合附图和优选实施例对本发明的技术方案进行详细地阐述。
本发明提供了一种采用LBM方法的大涡模拟方法,涉及CPU端,还涉及MIC众核协处理器端以及CPU+MIC协同计算模式,其中:
CPU端负责将要进行大涡模拟的网格数据分割,向MIC卡传递大涡模拟所需要的值,CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作,而且在整个网格的计算任务中,CPU也会以Openmp多线程模式,依次通过迁移碰撞,边界处理的多次迭代获取得速度、密度和流函数等宏观参量。
MIC众核协处理器负责网格点的迁移和碰撞过程,对边界进行处理,根据分布函数并行求得速度、密度和流函数等宏观参量。在MIC卡上也采用openmp多线程的方式来运算。
进一步地,CPU端将要进行大涡模拟的网格数据分割,向MIC卡传递大涡模拟所需要的值,以openmp多线程模式,并行执行迁移碰撞,边界处理来获取速度、密度和流函数等宏观参量,具体包括:
LES_MIC算法中,根据网格的划分,在对划分边界要做一些特殊的处理,在迁移过程中,每个点的计算需要其周围9个方向的分布函数值,因此,我们需要在分布函数中多存储1或2行用于下一步时层的计算,并且在每次迭代之后在节点与节点之间、MIC卡之间进行边界数据的交换。
表1 表示各MIC卡节点上数组大小
进一步地,CPU负责CPU+MIC协同计算模式的框架搭建以及任务调度,具体包括:
单节点服务器采用由双路6核CPU和2块KNF MIC卡组成的桌面服务器。在CPU+MICs协同计算中,我们把双路CPU和MIC卡都作为计算设备,这样每个单节点就相当于有3个计算设备,每个设备通过一个OpenMP线程进行控制。CPU+MIC协同计算框架如图3所示。
本专利是数据级并行,所以采用静态数据划分的方式,每次各个设备读取设定好的网格数据以及边界处理需要的数据,然后分别进行数据处理,相邻设备需要交换数据,迭代多次,直到所有设备计算完成所有网格数据,由CPU输出结果。
进一步地,MIC众核协处理器负责网格点的迁移和碰撞过程,对边界进行处理,根据分布函数并行求得速度、密度和流函数等宏观参量。在MIC卡上也采用openmp多线程的方式来运算。
根据对LBM算法中求解离散方程和边界处理的串行算法的热点分析和并行性分析,每个网格点的迁移、碰撞、宏观量统计、平衡态分布函数计算和边界处理的计算都是数据并行的。
进一步地,求解离散方程可以采用迁移碰撞的过程,宏观量统计、平衡态分布函数计算和碰撞过程中对每个网格的计算之间没有任何依赖性,因此,可以让MIC中的每个线程负责一个网格划分中的一行网格点的计算,每行网格点的计算利用MIC上的向量化技术进一步加速;分布函数的迁移只涉及到该格点周围的其他格点,也可以通过单个线程对全局存储器中相关分布函数的读操作来实现。
进一步地,在LBM算法中,对边界要做特殊的处理(非平衡外推,反弹),对于边界上的每个格点之间的计算也没有数据的依赖性,因此,可以利用OpenMP多线程负责边界格点的计算;
再进一步地,OpenMP的线程模型设计:根据MIC核心数目设置内核的线程数。
实施例
本发明首先分析用格子Boltzmann方法进行大涡模拟的性能瓶颈以及将串行方法移植到其他高性能平台上的难度,找到耗时的热点代码,测试其在大涡模拟的整个过程中所占的时间比例以及分析提高整个方法性能的难度以及开发轴端。
测试结果表明在采用格子Boltzmann方法进行大涡模拟时,大部分的时间消耗在求解离散方程和边界处理的过程,求解离散方程是通过迁移和碰撞的过程实现的,因此,求解离散方程和边界处理的迭代过程是格子Boltzmann方法中的性能瓶颈。通过发明内容中的分析,MIC中的每个线程负责一个网格划分中的一行网格点的计算,每行网格点的计算利用MIC上的向量化技术进一步加速;分布函数的迁移只涉及到该格点周围的其他格点,也可以通过单个线程对全局存储器中相关分布函数的读操作来实现。
整个算法代码实现的流程如图4所示。首先初始化密度,速度,平衡分布函数,对分布函数赋值,碰撞过程求出f* ,迁移过程,求出下一时刻的分布函数,由网格点上的分布函数统计密度和速度,利用求出的密度值和速度值计算平衡分布函数。
本发明针对上述分析,提供了一种基于格子Boltzmann方法且使用CPU+MIC协同计算模式的大涡模拟并行方法的实施例,包括如下步骤:
步骤一:CPU+MIC协同计算框架搭建
在单节点上,共有M+1个设备,包括一个CPU和M个MIC卡,采用OpenMP的fork-join模式搭建单节点上的框架,当程序开始执行的时候只有一个主线程存在,需要进行并行计算的时候,主线程派生出附加线程,即启用M+1个OpenMP线程, 0~M-1号进程控制所有MIC设备,M号线程控制CPU设备,按照数据的分配设计,每个设备进行输入数据的分发读写,
如在本专利中,CPU+2MICs平台上, 主线程控制输入数据的动态分发,0号线程控制MIC0设备,1号线程控制MIC1设备,2号线程控制CPU设备,如图3所示。
在CPU和MIC卡上,各个设备的数据大小的划分如表2所示,对于分布函数,其大小要比网格大小大两行,同时方便代码的书写,我们对最上面的设备和最下面的设备也都做了+2行的数据定义(有一行在代码中不需要使用)。
单节点上CPU+MICs协同计算的外围框架的伪代码如下所示:
定义一些变量
//在迭代范围之内
int DEVICE_NUM //设备数
for(int i=0;i<steps;i++)
{
if(i%2==0) //奇数步和偶数步输入和输出交换,所以有一个迭代的判断
{
omp_set_nested(true);
#pragma omp parallel for private(…), num_threads(DEVICE_NUM)
for(int thread=0;thread<DEVICE_NUM;thread++)
{
if(thread==0) //mic 0 computing
{
#pragma offload target(mic:0) \
in(fs0_in0_up,…:length(nx) alloc_if(0) free_if(0))\
out(fn1_out0_up,:length(nx) alloc_if(0) free_if(0))\
nocopy(fr0_mic0,…:length((hh+1)*nx) alloc_if(0) free_if(0)) \
nocopy(fr1_mic0,…:length((hh+1)*nx) alloc_if(0) free_if(0))
{
……
LBCollProp(DEVICE_NUM, thread_num_mic,thread,… );
LBBC_LR(DEVICE_NUM, thread_num_mic,thread, …);
LBBC_DOWN(DEVICE_NUM, thread_num_mic,thread , …);
……
}//mic0 end
}
else if(thread==DEVICE_NUM-1) //cpu computing
{
……
LBCollProp(DEVICE_NUM, thread_num_omp,thread,…);
LBBC_LR(DEVICE_NUM, thread_num_omp,thread,…);
LBBC_UP(DEVICE_NUM, thread_num_omp,thread,…);
}
else //other mic computing
{
#pragma offload target(mic:1) \
in(thread_num_mic,thread,nx,…)\
in(fs1_in1_up,…:length(nx) alloc_if(0) free_if(0))\
out(fn0_out1_up,…:length(nx) alloc_if(0) free_if(0))\
nocopy(fr0_mic1,…:length((hh+2)*nx) alloc_if(0) free_if(0)) \
…
{//mic1 compute
……
LBCollProp(DEVICE_NUM, thread_num_mic,thread,… );
LBBC_LR(DEVICE_NUM, thread_num_mic,thread, …);
LBBC_DOWN(DEVICE_NUM, thread_num_mic,thread, nx, hh, fr0_mic0, fe0_m …);
……
}
else //奇数步
{
//函数内容和偶数步是一样的,只是输入输出和偶数步互换
}
步骤二:CPU/MIC内核实现
(3) 设计迁移碰撞内核,设计线程数为T=4*M,M为MIC卡的核数,并且让内核中的每个线程计算一行网格点的迁移和碰撞过程,以及利用#pragma ivdep实现内核中内层循环的向量化,如图4所示,内核伪代码如下;
1: #pragma omp parallel for private (i, j, k,…) num_threads(T)//T为线程数
2: for (i=1;i<NY-1;i++)
3: #pragma ivdep //向量化
4: for(j=1;j<NX-1;j++)
5: {
6: k=i*NX+j; //k表示网格的标号
7: fr = fr0[k];//0代表上一时层的分布函数
8: fe = fe0[k-1];
9: fn = fn0[k-NX];
10: fw = fw0[k+1];
11: fs = fs0[k+NX];
12: fne = fne0[k-NX-1];
13: fnw = fnw0[k-NX+1];
14: fsw = fsw0[k+NX+1];
15: fse = fse0[k+NX-1];
16: /*碰撞过程*/
17:根据迁移后的分布函数fr-fse 求宏观量
18:根据宏观量求各个方向的平衡分布函数f1,f2,f3,f4,f5,f6,f7,f8;
19: 根据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];
20: }
在MIC端对边界进行处理,边界处理可以采用反弹法、非平衡外推法等方法,在对边界的处理时同样设计T个线程处理边界节点的计算;
步骤三:数据传递方式设计
节点间的数据传递
LES算法根据格点所在区域将格点划分到不同的设备上,因此当每个格点更新自己的分布函数并迁移时,各个划分区域中边界格点的分布函数要传递给邻近的节点。如图5所示。图5中位于H、L虚线之间的黑实线将计算域一分为二,但为了在计算中方便该实线的迁移的计算,需要邻近网格,M(i+1,j)的计算域增加了该实线下边邻近网格L,M(i,j)的计算域增加了该实线上边邻近网格H。M(i+1,j)的计算域中H的分布函数需要传递给M(i,j),同时也需要接收从M(i,j)传递过来的L的分布函数;对M(i,j)也类似处理。
具体优化1:LES_MIC优化,实施过程如下:
1)向量化
对于MIC内核,我们设计外层for并行,内层for采用向量化的方案,对于各个内核函数,我们可以采取自动向量化的方案进行优化;
2)减少offload次数
在迭代过程中,尽量减少offload次数,以及CPU和MIC之间的I/O次数;
3)减少节点与节点的数据传递以及MIC卡之间的数据传递
在每次迭代之后,相邻节点之间,以及相邻的MIC卡之间要进行边界数据的传递,而每个网格有9个方向的分布函数,然而在内核计算中,并不需要边界上的全部9个方向上的值,只需要3个方向的值即可,如图6所示。对于节点i,只需要接收节H中的fsw,fs,fse的值,同样,对于节点i+1,只需要接收节L中的fnw,fn,fne的值。
本项目中LES应用案例测试用的软硬件环境如表3所示。
LES性能采用格子点更新速率(LUPS, Lattice Unit Per Second)统计,常为MLUPS(每秒更新的百万网格数),计算方法:
P=NX*NY*S/T
其中NX、NY为网格宽和高,S为流场迭代次数,T为计算时间,P为格子点更新速率。表3是LES实验环境
平台 | Inspur NF5280M3,2节点 |
CPU | Intel Xeon CPU E5-2680 2.7GHz,双路8核 |
Memory | DDR3 1333MHz 128GB |
MIC | KNC,61核,1.1GHz,GDDR5 8GB memory,5.5GT/s |
OS | Red Hat Enterprise Linux Server release 6.1,64bit |
驱动 | KNC_beta_oem-2.1.3653-8-rhel-6.1 |
编译器 | l_ccompxe_intel64_2013.0.079 |
测试用例 | 雷诺数:10000;迭代次数:10000 |
基于单节点CPU+MIC协同计算LES并行算法的实验结果如图7所示,
由本发明给出的上述测试结果可以看出,在MIC众核架构平台上,基于格子Boltzmann算法和CPU+MIC异构并行架构的基础上,可以在较短周期内,较容易地大幅度加速大涡模拟的运算。
以上说明仅为本发明较佳的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应该以权力要求书的保护范围为准。
除说明书所述的技术特征外,均为本专业技术人员的已知技术。
Claims (1)
1.一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法, 其特征在于包括采用LBM方法的大涡模拟方法、CPU端、MIC众核协处理器端以及CPU+MIC协同计算模式,其中:
CPU端负责将要进行大涡模拟的网格数据分割,向MIC卡传递大涡模拟所需要的值,CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作,而且在整个网格的计算任务中,CPU也会以Openmp多线程模式,依次通过迁移碰撞,边界处理的多次迭代获取得速度、密度和流函数的宏观参量;
MIC众核协处理器负责网格点的迁移和碰撞过程,对边界进行处理,根据分布函数并行求得速度、密度和流函数等宏观参量,在MIC卡上也采用openmp多线程的方式来运算;
CPU端将要进行大涡模拟的网格数据分割,向MIC卡传递大涡模拟所需要的值,以openmp多线程模式,并行执行迁移碰撞,边界处理来获取速度、密度和流函数的宏观参量,具体包括:
LES_MIC算法中,根据网格的划分,在对划分边界要做一些特殊的处理,在迁移过程中,每个点的计算需要其周围9个方向的分布函数值,因此,需要在分布函数中多存储1或2行用于下一步时层的计算,并且在每次迭代之后在节点与节点之间、MIC卡之间进行边界数据的交换;
CPU负责CPU+MIC协同计算模式的框架搭建以及任务调度,具体包括:
单节点服务器采用由双路6核CPU和2块KNF MIC卡组成的桌面服务器,在CPU+MICs协同计算中,把双路CPU和MIC卡都作为计算设备,这样每个单节点就相当于有3个计算设备,每个设备通过一个OpenMP线程进行控制;
本方法是数据级并行,所以采用静态数据划分的方式,每次各个设备读取设定好的网格数据以及边界处理需要的数据,然后分别进行数据处理,相邻设备需要交换数据,迭代多次,直到所有设备计算完成所有网格数据,由CPU输出结果;
MIC众核协处理器负责网格点的迁移和碰撞过程,对边界进行处理,根据分布函数并行求得速度、密度和流函数的宏观参量,在MIC卡上也采用openmp多线程的方式来运算;
根据对LBM算法中求解离散方程和边界处理的串行算法的热点分析和并行性分析,每个网格点的迁移、碰撞、宏观量统计、平衡态分布函数计算和边界处理的计算都是数据并行的;
求解离散方程可以采用迁移碰撞的过程,宏观量统计、平衡态分布函数计算和碰撞过程中对每个网格的计算之间没有任何依赖性,因此,让MIC中的每个线程负责一个网格划分中的一行网格点的计算,每行网格点的计算利用MIC上的向量化技术进一步加速;分布函数的迁移只涉及到该格点周围的其他格点,或通过单个线程对全局存储器中相关分布函数的读操作来实现;
在LBM算法中,对边界要做特殊的处理,包括:非平衡外推,反弹,对于边界上的每个格点之间的计算也没有数据的依赖性,因此,利用OpenMP多线程负责边界格点的计算;
OpenMP的线程模型设计:根据MIC核心数目设置内核的线程数;
具体步骤如下:
步骤一:CPU+MIC协同计算框架搭建
在单节点上,共有M+1个设备,一个CPU+M个MIC卡,采用OpenMP的fork-join模式搭建单节点上的框架,当程序开始执行的时候只有一个主线程存在,需要进行并行计算的时候,主线程派生出附加线程,即启用M+1个OpenMP线程, 0~M-1号进程控制所有MIC设备,M号线程控制CPU设备,按照数据的分配设计,每个设备进行输入数据的分发读写;
在CPU+2MICs平台上, 主线程控制输入数据的动态分发,0号线程控制MIC0设备,1号线程控制MIC1设备,2号线程控制CPU设备;
在CPU和MIC卡上,对于分布函数,其大小要比网格大小大两行,同时方便代码的书写,要对最上面的设备和最下面的设备要做+2行的数据定义,其中有一行在代码中不需要使用;
单节点上CPU+MICs协同计算的外围框架的伪代码如下所示:
定义一些变量
//在迭代范围之内
int DEVICE_NUM //设备数
for(int i=0;i<steps;i++)
{
if(i%2==0) //奇数步和偶数步输入和输出交换,所以有一个迭代的判断
{
omp_set_nested(true);
#pragma omp parallel for private(…), num_threads(DEVICE_NUM)
for(int thread=0;thread<DEVICE_NUM;thread++)
{
if(thread==0) //mic 0 computing
{
#pragma offload target(mic:0) \
in(fs0_in0_up,…:length(nx) alloc_if(0) free_if(0))\
out(fn1_out0_up,:length(nx) alloc_if(0) free_if(0))\
nocopy(fr0_mic0,…:length((hh+1)*nx) alloc_if(0) free_if(0)) \
nocopy(fr1_mic0,…:length((hh+1)*nx) alloc_if(0) free_if(0))
{
……
LBCollProp(DEVICE_NUM, thread_num_mic,thread,… );
LBBC_LR(DEVICE_NUM, thread_num_mic,thread, …);
LBBC_DOWN(DEVICE_NUM, thread_num_mic,thread , …);
……
}//mic0 end
}
else if(thread==DEVICE_NUM-1) //cpu computing
{
……
LBCollProp(DEVICE_NUM, thread_num_omp,thread,…);
LBBC_LR(DEVICE_NUM, thread_num_omp,thread,…);
LBBC_UP(DEVICE_NUM, thread_num_omp,thread,…);
}
else //other mic computing
{
#pragma offload target(mic:1) \
in(thread_num_mic,thread,nx,…)\
in(fs1_in1_up,…:length(nx) alloc_if(0) free_if(0))\
out(fn0_out1_up,…:length(nx) alloc_if(0) free_if(0))\
nocopy(fr0_mic1,…:length((hh+2)*nx) alloc_if(0) free_if(0)) \
…
{//mic1 compute
……
LBCollProp(DEVICE_NUM, thread_num_mic,thread,… );
LBBC_LR(DEVICE_NUM, thread_num_mic,thread, …);
LBBC_DOWN(DEVICE_NUM, thread_num_mic,thread, nx, hh, fr0_mic0, fe0_m …);
……
}
else //奇数步
{
//函数内容和偶数步是一样的,只是输入输出和偶数步互换
}
步骤二:CPU/MIC内核实现
设计迁移碰撞内核,设计线程数为T=4*M,M为MIC卡的核数,并且让内核中的每个线程计算一行网格点的迁移和碰撞过程,以及利用#pragma ivdep实现内核中内层循环的向量化,内核伪代码如下;
1: #pragma omp parallel for private (i, j, k,…) num_threads(T)//T为线程数
2: for (i=1;i<NY-1;i++)
3: #pragma ivdep //向量化
4: for(j=1;j<NX-1;j++)
5: {
6: k=i*NX+j; //k表示网格的标号
7: fr = fr0[k];//0代表上一时层的分布函数
8: fe = fe0[k-1];
9: fn = fn0[k-NX];
10: fw = fw0[k+1];
11: fs = fs0[k+NX];
12: fne = fne0[k-NX-1];
13: fnw = fnw0[k-NX+1];
14: fsw = fsw0[k+NX+1];
15: fse = fse0[k+NX-1];
16: /*碰撞过程*/
17:根据迁移后的分布函数fr-fse 求宏观量
18: 根据宏观量求各个方向的平衡分布函数f1,f2,f3,f4,f5,f6,f7,f8;
19: 根据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];
20: }
在MIC端对边界进行处理,边界处理可以采用反弹法、非平衡外推法等方法,在对边界的处理时同样设计T个线程处理边界节点的计算;
步骤三:数据传递方式设计
节点间的数据传递:LES算法根据格点所在区域将格点划分到不同的设备上,因此当每个格点更新自己的分布函数并迁移时,各个划分区域中边界格点的分布函数要传递给邻近的节点;
LES_MIC优化过程如下:
1)向量化
对于MIC内核,设计外层for并行,内层for采用向量化的方案,对于各个内核函数,采取自动向量化的方案进行优化;
2)减少offload次数
在迭代过程中,尽量减少offload次数,以及CPU和MIC之间的I/O次数;
3)减少节点与节点的数据传递以及MIC卡之间的数据传递,在每次迭代之后,相邻节点之间,以及相邻的MIC卡之间要进行边界数据的传递,而每个网格有9个方向的分布函数,然而在内核计算中,并不需要边界上的全部9个方向上的值,只需要3个方向的值即可,如图6所示,对于节点i,只需要接收节H中的fsw,fs,fse的值,同样,对于节点i+1,只需要接收节L中的fnw,fn,fne的值;
LES性能采用格子点更新速率LUPS, Lattice Unit Per Second统计,常为MLUPS,每秒更新的百万网格数,计算方法:
P=NX*NY*S/T
其中NX、NY为网格宽和高,S为流场迭代次数,T为计算时间,P为格子点更新速率;
通过测试结果看出,在MIC众核架构平台上,基于格子Boltzmann算法和CPU+MIC异构并行架构的基础上,能在较短周期内,较容易地大幅度加速大涡模拟的运算。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN2013102291613A CN103324531A (zh) | 2013-06-09 | 2013-06-09 | 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN2013102291613A CN103324531A (zh) | 2013-06-09 | 2013-06-09 | 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN103324531A true CN103324531A (zh) | 2013-09-25 |
Family
ID=49193293
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN2013102291613A Pending CN103324531A (zh) | 2013-06-09 | 2013-06-09 | 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103324531A (zh) |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103778098A (zh) * | 2014-02-17 | 2014-05-07 | 浪潮(北京)电子信息产业有限公司 | 基于格子Boltzmann理论实现协同计算大涡模拟系统及方法 |
CN104331320A (zh) * | 2014-10-30 | 2015-02-04 | 浪潮电子信息产业股份有限公司 | 一种利用mic加速三维纵横波分离的弹性波方程模拟的方法 |
CN105787227A (zh) * | 2016-05-11 | 2016-07-20 | 中国科学院近代物理研究所 | 结构材料辐照损伤的多gpu分子动力学模拟方法 |
CN105893151A (zh) * | 2016-04-01 | 2016-08-24 | 浪潮电子信息产业股份有限公司 | 一种基于cpu+mic异构平台的高维数据流的处理方法 |
CN106383961A (zh) * | 2016-09-29 | 2017-02-08 | 中国南方电网有限责任公司电网技术研究中心 | Cpu+mic异构平台下的大涡模拟算法优化处理方法 |
CN106487784A (zh) * | 2016-09-28 | 2017-03-08 | 东软集团股份有限公司 | 一种会话迁移的方法、装置及防火墙 |
CN107102895A (zh) * | 2016-02-19 | 2017-08-29 | 中国石油化工股份有限公司 | 一种并行网格处理器自适应分配方法和系统 |
CN107515987A (zh) * | 2017-08-25 | 2017-12-26 | 中国地质大学(北京) | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
CN107636637A (zh) * | 2015-04-17 | 2018-01-26 | 微软技术许可有限责任公司 | 用于使用软处理器执行软件线程的系统和方法 |
CN108595277A (zh) * | 2018-04-08 | 2018-09-28 | 西安交通大学 | 一种基于OpenMP/MPI混合编程的CFD仿真程序的通信优化方法 |
CN109408867A (zh) * | 2018-09-12 | 2019-03-01 | 西安交通大学 | 一种基于mic协处理器的显式r-k时间推进加速方法 |
CN111105341A (zh) * | 2019-12-16 | 2020-05-05 | 上海大学 | 一种低功耗高运算性能求解计算流体动力学的框架方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101354729A (zh) * | 2007-07-23 | 2009-01-28 | 南车四方机车车辆股份有限公司 | 高速列车头部纵向对称面型线的低气动噪声优化方法 |
CN103064819A (zh) * | 2012-10-25 | 2013-04-24 | 浪潮电子信息产业股份有限公司 | 一种利用MIC快速实现格子Boltzmann并行加速的方法 |
US20130136788A1 (en) * | 2007-05-07 | 2013-05-30 | Insmed Incorporated | Method for treating pulmonary disorders with liposomal amikacin formulations |
-
2013
- 2013-06-09 CN CN2013102291613A patent/CN103324531A/zh active Pending
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130136788A1 (en) * | 2007-05-07 | 2013-05-30 | Insmed Incorporated | Method for treating pulmonary disorders with liposomal amikacin formulations |
CN101354729A (zh) * | 2007-07-23 | 2009-01-28 | 南车四方机车车辆股份有限公司 | 高速列车头部纵向对称面型线的低气动噪声优化方法 |
CN103064819A (zh) * | 2012-10-25 | 2013-04-24 | 浪潮电子信息产业股份有限公司 | 一种利用MIC快速实现格子Boltzmann并行加速的方法 |
Non-Patent Citations (1)
Title |
---|
王恩东 等: "《MIC高性能计算编程指南》", 30 November 2012, 中国水利水电出版社 * |
Cited By (18)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103778098A (zh) * | 2014-02-17 | 2014-05-07 | 浪潮(北京)电子信息产业有限公司 | 基于格子Boltzmann理论实现协同计算大涡模拟系统及方法 |
CN104331320A (zh) * | 2014-10-30 | 2015-02-04 | 浪潮电子信息产业股份有限公司 | 一种利用mic加速三维纵横波分离的弹性波方程模拟的方法 |
CN107636637A (zh) * | 2015-04-17 | 2018-01-26 | 微软技术许可有限责任公司 | 用于使用软处理器执行软件线程的系统和方法 |
CN107102895A (zh) * | 2016-02-19 | 2017-08-29 | 中国石油化工股份有限公司 | 一种并行网格处理器自适应分配方法和系统 |
CN105893151B (zh) * | 2016-04-01 | 2019-03-08 | 浪潮电子信息产业股份有限公司 | 一种基于cpu+mic异构平台的高维数据流的处理方法 |
CN105893151A (zh) * | 2016-04-01 | 2016-08-24 | 浪潮电子信息产业股份有限公司 | 一种基于cpu+mic异构平台的高维数据流的处理方法 |
CN105787227B (zh) * | 2016-05-11 | 2018-10-09 | 中国科学院近代物理研究所 | 结构材料辐照损伤的多gpu分子动力学模拟方法 |
CN105787227A (zh) * | 2016-05-11 | 2016-07-20 | 中国科学院近代物理研究所 | 结构材料辐照损伤的多gpu分子动力学模拟方法 |
CN106487784A (zh) * | 2016-09-28 | 2017-03-08 | 东软集团股份有限公司 | 一种会话迁移的方法、装置及防火墙 |
CN106487784B (zh) * | 2016-09-28 | 2019-06-25 | 东软集团股份有限公司 | 一种会话迁移的方法、装置及防火墙 |
CN106383961B (zh) * | 2016-09-29 | 2019-07-19 | 中国南方电网有限责任公司电网技术研究中心 | Cpu+mic异构平台下的大涡模拟算法优化处理方法 |
CN106383961A (zh) * | 2016-09-29 | 2017-02-08 | 中国南方电网有限责任公司电网技术研究中心 | Cpu+mic异构平台下的大涡模拟算法优化处理方法 |
CN107515987A (zh) * | 2017-08-25 | 2017-12-26 | 中国地质大学(北京) | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
CN108595277A (zh) * | 2018-04-08 | 2018-09-28 | 西安交通大学 | 一种基于OpenMP/MPI混合编程的CFD仿真程序的通信优化方法 |
CN108595277B (zh) * | 2018-04-08 | 2021-01-19 | 西安交通大学 | 一种基于OpenMP/MPI混合编程的CFD仿真程序的通信优化方法 |
CN109408867A (zh) * | 2018-09-12 | 2019-03-01 | 西安交通大学 | 一种基于mic协处理器的显式r-k时间推进加速方法 |
CN109408867B (zh) * | 2018-09-12 | 2021-04-20 | 西安交通大学 | 一种基于mic协处理器的显式r-k时间推进加速方法 |
CN111105341A (zh) * | 2019-12-16 | 2020-05-05 | 上海大学 | 一种低功耗高运算性能求解计算流体动力学的框架方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103324531A (zh) | 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法 | |
Pérez-Hurtado et al. | A membrane parallel rapidly-exploring random tree algorithm for robotic motion planning | |
Rinaldi et al. | A Lattice-Boltzmann solver for 3D fluid simulation on GPU | |
Xiong et al. | GPU-accelerated adaptive particle splitting and merging in SPH | |
Rehbach et al. | Comparison of parallel surrogate-assisted optimization approaches | |
CN102681972A (zh) | 一种利用GPU加速格子-Boltzmann的方法 | |
CN103064819A (zh) | 一种利用MIC快速实现格子Boltzmann并行加速的方法 | |
Gunow et al. | Simplemoc-a performance abstraction for 3d moc | |
Correia et al. | Abelian–Higgs cosmic string evolution with multiple GPUs | |
CN103778098A (zh) | 基于格子Boltzmann理论实现协同计算大涡模拟系统及方法 | |
Stone et al. | GPGPU parallel algorithms for structured-grid CFD codes | |
CN108460195A (zh) | 海啸数值计算模型基于gpu并行的快速执行方法 | |
Koskela et al. | A novel multi-level integrated roofline model approach for performance characterization | |
McClure et al. | Petascale application of a coupled CPU-GPU algorithm for simulation and analysis of multiphase flow solutions in porous medium systems | |
Malinowski et al. | Multi-agent large-scale parallel crowd simulation | |
Obrecht et al. | The TheLMA project: Multi-GPU implementation of the lattice Boltzmann method | |
Al-Hashimi et al. | Evaluating power and energy efficiency of bitonic mergesort on graphics processing unit | |
Yamazaki et al. | New scheduling strategies and hybrid programming for a parallel right-looking sparse LU factorization algorithm on multicore cluster systems | |
Charlton et al. | Fast simulation of crowd collision avoidance | |
Mallinson et al. | Experiences at scale with pgas versions of a hydrodynamics application | |
Coghlan et al. | Argonne applications for the IBM blue gene/Q, Mira | |
Chapuis et al. | Predicting performance of smoothed particle hydrodynamics codes at large scales | |
Holmen et al. | Exploring use of the reserved core | |
Kulikov et al. | Numerical modeling of jellyfish galaxy at intel xeon phi supercomputers | |
Pera | Design and performance evaluation of a Linux HPC cluster |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
WD01 | Invention patent application deemed withdrawn after publication | ||
WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20130925 |