[发明专利]一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法

权利要求书

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异构并行架构的基础上，能在较短周期内，较容易地大幅度加速大涡模拟的运算。