代码之家 › 专栏 › 技术社区 › 2b-t

OpenMP的扩展问题

scientific-computing physics openmp parallel-processing c++

2b-t · 技术社区 · 8 年前

我为一种特殊类型的三维cfd模拟编写了一个代码,即格子boltzmann方法(与timm kr_四分之一ger et alii的《格子boltzmann方法》一书中提供的代码非常相似)。用openmp多线程处理程序我遇到了一些我不太理解的问题:结果很大程度上依赖于整个域的大小。

其基本原理是,在离散方向上,一个3d域的每个单元被赋予19个分布函数(0-18)的特定值。它们被放置在堆中分配的两个线性数组中(一个填充被放置在一个单独的数组中):某个单元格的18个填充在内存中是连续的,连续x值的值彼此相邻,依此类推(所以行主要排序:填充->x->y->z)。这些分布函数根据单元格中的某些值重新分布,然后流到相邻的单元格。因此我有两个种群f1和f2。该算法从f1获取值,重新分配并将其复制到f2。然后指针交换,算法重新开始。代码在单个核上运行得非常好,但是当我尝试在多个核上并行时,我得到的性能取决于域的总体大小:对于非常小的域(10^3个单元),算法的速度相当慢,每秒有1500万个单元,对于非常小的域(30^3个单元),该算法的速度非常快,每秒超过6000万个单元,对于任何大于此速度的域,性能再次下降到每秒约3000万个单元。在单个核心上执行代码只会导致每秒1500万个单元的性能相同。当然,这些结果在不同的处理器之间有所不同,但在质量上仍然存在相同的问题!

代码的核心归结为这个反复执行的并行循环,指向f1和f2的指针被交换:

#pragma omp parallel for default(none) shared(f0,f1,f2) schedule(static)
    for(unsigned int z = 0; z < NZ; ++z)
    {
        for(unsigned int y = 0; y < NY; ++y)
        {
            for(unsigned int x = 0; x < NX; ++x)
            {
                /// temporary populations
                double ft0  = f0[D3Q19_ScalarIndex(x,y,z)];
                double ft1  = f1[D3Q19_FieldIndex(x,y,z,1)];
                double ft2  = f1[D3Q19_FieldIndex(x,y,z,2)];
                double ft3  = f1[D3Q19_FieldIndex(x,y,z,3)];
                double ft4  = f1[D3Q19_FieldIndex(x,y,z,4)];
                double ft5  = f1[D3Q19_FieldIndex(x,y,z,5)];
                double ft6  = f1[D3Q19_FieldIndex(x,y,z,6)];
                double ft7  = f1[D3Q19_FieldIndex(x,y,z,7)];
                double ft8  = f1[D3Q19_FieldIndex(x,y,z,8)];
                double ft9  = f1[D3Q19_FieldIndex(x,y,z,9)];
                double ft10 = f1[D3Q19_FieldIndex(x,y,z,10)];
                double ft11 = f1[D3Q19_FieldIndex(x,y,z,11)];
                double ft12 = f1[D3Q19_FieldIndex(x,y,z,12)];
                double ft13 = f1[D3Q19_FieldIndex(x,y,z,13)];
                double ft14 = f1[D3Q19_FieldIndex(x,y,z,14)];
                double ft15 = f1[D3Q19_FieldIndex(x,y,z,15)];
                double ft16 = f1[D3Q19_FieldIndex(x,y,z,16)];
                double ft17 = f1[D3Q19_FieldIndex(x,y,z,17)];
                double ft18 = f1[D3Q19_FieldIndex(x,y,z,18)];

                /// microscopic to macroscopic
                double r    = ft0 + ft1 + ft2 + ft3 + ft4 + ft5 + ft6 + ft7 + ft8 + ft9 + ft10 + ft11 + ft12 + ft13 + ft14 + ft15 + ft16 + ft17 + ft18;
                double rinv = 1.0/r;
                double u    = rinv*(ft1 - ft2 + ft7 + ft8  + ft9   + ft10 - ft11 - ft12 - ft13 - ft14);
                double v    = rinv*(ft3 - ft4 + ft7 - ft8  + ft11  - ft12 + ft15 + ft16 - ft17 - ft18);
                double w    = rinv*(ft5 - ft6 + ft9 - ft10 + ft13 -  ft14 + ft15 - ft16 + ft17 - ft18);

                /// collision & streaming
                double trw0 = omega*r*w0;                   //temporary variables
                double trwc = omega*r*wc;
                double trwd = omega*r*wd;
                double uu   = 1.0 - 1.5*(u*u+v*v+w*w);

                double bu = 3.0*u;
                double bv = 3.0*v;
                double bw = 3.0*w;

                unsigned int xp = (x + 1) % NX;             //calculate x,y,z coordinates of neighbouring cells
                unsigned int yp = (y + 1) % NY;
                unsigned int zp = (z + 1) % NZ;
                unsigned int xm = (NX + x - 1) % NX;
                unsigned int ym = (NY + y - 1) % NY;
                unsigned int zm = (NZ + z - 1) % NZ;

                f0[D3Q19_ScalarIndex(x,y,z)]      = bomega*ft0  + trw0*(uu);                        //redistribute distribution functions and stream to neighbouring cells
                double cu = bu;
                f2[D3Q19_FieldIndex(xp,y, z,  1)] = bomega*ft1  + trwc*(uu + cu*(1.0 + 0.5*cu));
                cu = -bu;
                f2[D3Q19_FieldIndex(xm,y, z,  2)] = bomega*ft2  + trwc*(uu + cu*(1.0 + 0.5*cu));
                cu = bv;
                f2[D3Q19_FieldIndex(x, yp,z,  3)] = bomega*ft3  + trwc*(uu + cu*(1.0 + 0.5*cu));
                cu = -bv;
                f2[D3Q19_FieldIndex(x, ym,z,  4)] = bomega*ft4  + trwc*(uu + cu*(1.0 + 0.5*cu));
                cu = bw;
                f2[D3Q19_FieldIndex(x, y, zp, 5)] = bomega*ft5  + trwc*(uu + cu*(1.0 + 0.5*cu));
                cu = -bw;
                f2[D3Q19_FieldIndex(x, y, zm, 6)] = bomega*ft6  + trwc*(uu + cu*(1.0 + 0.5*cu));
                cu = bu+bv;
                f2[D3Q19_FieldIndex(xp,yp,z,  7)] = bomega*ft7  + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = bu-bv;
                f2[D3Q19_FieldIndex(xp,ym,z,  8)] = bomega*ft8  + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = bu+bw;
                f2[D3Q19_FieldIndex(xp,y, zp, 9)] = bomega*ft9  + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = bu-bw;
                f2[D3Q19_FieldIndex(xp,y, zm,10)] = bomega*ft10 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = -bu+bv;
                f2[D3Q19_FieldIndex(xm,yp,z, 11)] = bomega*ft11 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = -bu-bv;
                f2[D3Q19_FieldIndex(xm,ym,z, 12)] = bomega*ft12 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = -bu+bw;
                f2[D3Q19_FieldIndex(xm,y, zp,13)] = bomega*ft13 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = -bu-bw;
                f2[D3Q19_FieldIndex(xm,y, zm,14)] = bomega*ft14 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = bv+bw;
                f2[D3Q19_FieldIndex(x, yp,zp,15)] = bomega*ft15 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = bv-bw;
                f2[D3Q19_FieldIndex(x, yp,zm,16)] = bomega*ft16 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = -bv+bw;
                f2[D3Q19_FieldIndex(x, ym,zp,17)] = bomega*ft17 + trwd*(uu + cu*(1.0 + 0.5*cu));
                cu = -bv-bw;
                f2[D3Q19_FieldIndex(x, ym,zm,18)] = bomega*ft18 + trwd*(uu + cu*(1.0 + 0.5*cu));
            }
        }
    }

如果有人能给我一些提示,告诉我如何找到这种特殊行为的原因,甚至知道是什么导致了这个问题,那就太棒了。如果需要,我可以提供一个完整版本的简化代码! 提前多谢!

1 回复 | 直到 8 年前

dlasalle 8 年前

在共享内存系统(单台机器上的线程代码)上实现扩展是相当棘手的,并且通常需要大量的调整。代码中可能发生的情况是,每个线程的域的一部分适合缓存中“非常小”的问题大小,但随着NX和NY中问题大小的增加,每个线程的数据将停止适合缓存。

为了避免这样的问题,最好将域分解为固定大小的块,这些块的大小不会随域而改变,而是在数量上改变。

const unsigned int numBlocksZ = std::ceil(static_cast<double>(NZ) / BLOCK_SIZE);
const unsigned int numBlocksY = std::ceil(static_cast<double>(NY) / BLOCK_SIZE);
const unsigned int numBlocksX = std::ceil(static_cast<double>(NX) / BLOCK_SIZE);

#pragma omp parallel for default(none) shared(f0,f1,f2) schedule(static,1)
for(unsigned int block = 0; block < numBlocks; ++block)
{
  unsigned int startZ = BLOCK_SIZE* (block / (numBlocksX*numBlocksY));
  unsigned int endZ = std::min(startZ + BLOCK_SIZE, NZ);
  for(unsigned int z = startZ; z < endZ; ++z) {
    unsigned int startY = BLOCK_SIZE*(((block % (numBlocksX*numBlocksY)) / numBlocksX);
    unsigned int endY = std::min(startY + BLOCK_SIZE, NY);
    for(unsigned int y = startY; y < endY; ++y)
    {
      unsigned int startX = BLOCK_SIZE(block % numBlocksX);
      unsigned int endX = std::min(startX + BLOCK_SIZE, NX);
      for(unsigned int x = startX; x < endX; ++x)
      {
        ...
      }
    }  
  }

像上面这样的方法也应该增加 cache locality 通过使用3d blocking(假设这是一个3d模板操作),进一步提高性能。你需要调整block_的大小,以找到在给定系统上能给你最好性能的东西(我会从较小的开始,增加2的功率,例如4,8,16…)。