c - 从随机分布的粒子到规则网格的最佳并行通信

Question

我正在并行化我的粒子单元代码，我用它来模拟地球内部的 2D 和 3D 变形。代码的几个例程很容易使用 OpenMP 并行化并且扩展性很好。但是，我在处理从粒子到网格单元的插值的代码的关键部分遇到了问题。粒子在每次迭代中四处移动(根据速度场)。许多计算在规则的、不变形的网格上执行是最有效的。因此，每次迭代都涉及从“随机”分布的粒子到网格单元的通信。

问题可以用以下简化的一维代码来说明:

//EXPLANATION OF VARIABLES (all previously allocated and initialized, 1D arrays)
//double *markerval; // Size Nm. Particle values. Are to be interpolated to the grid
//double *grid; // Size Ng=Nm/100 Grid values. 
//uint *markerpos; // Size Nm. Position of particles relative to grid (each particle
// knows what grid cell it belongs to) possible values are 0,1,...Ng-1

//#pragma omp parallel for schedule(static) private(e)
for (e=0; e<Nm; e++) {
    //#pragma omp atomic
    grid[markerpos[e]]+=markerval[e];
}

在最坏的情况下，粒子位置是随机的，但通常情况下，粒子在内存中彼此相邻，在空间中也彼此相邻，因此在网格内存中也是如此。

如何有效地并行化此过程？多个粒子映射到同一个网格单元，因此如果上述循环直接并行化，则竞争条件和缓存交换的可能性非零。使更新原子化可以防止竞争条件，但会使代码比顺序情况慢得多。

我还尝试为每个线程制作一份网格值的私有(private)副本，然后将它们相加。然而，这可能需要在代码中使用太多内存，并且对于这个例子，它不能很好地扩展线程数量(我不确定原因)。

第三种选择可能是从网格映射到粒子，然后循环遍历网格索引而不是粒子索引。但是，我担心这会涉及很多，并且需要对代码进行重大更改，而且我不确定它会有多大帮助，因为它还需要使用计算量也很大的排序例程。

有人遇到过这个或类似的问题吗？

Answer 1

一个选项可以是在线程上手动映射迭代:

#pragma omp parallel shared(Nm,Ng,markerval,markerpos,grid)
{
  int nthreads = omp_get_num_threads();
  int rank     = omp_get_thread_num();
  int factor   = Ng/nthreads;

  for (int e = 0; e < Nm; e++) {
    int pos = markerpos[e];
    if ( (pos/factor)%nthreads == rank )
      grid[pos]+=markerval[e];
  }
}

几点说明:

1. for 循环的迭代不在线程之间共享。相反，每个线程 执行所有迭代。
2. for 循环内的条件决定哪个线程 将更新grid 数组的位置pos。此位置仅属于一个线程，因此不需要atomic 保护。
3. 公式 (pos/factor)%nthreads 只是一个简单的启发式算法。返回 0,...,nthreads-1 范围内值的 pos 的任何函数实际上都可以替换为该表达式，而不会影响最终结果的有效性(所以如果你有更好的镜头，请随意更改它)。请注意，此功能选择不当可能会导致负载平衡问题。