algorithm - CUDA 粒子中的最近邻

Question

编辑 2: 请查看 this crosspost对于 TLDR。

编辑:鉴于粒子被分割成网格单元(比如 16^3 网格)，让每个网格单元运行一个工作组和尽可能多的工作是更好的主意吗？一个工作组中的项目，因为每个网格单元可以有最大数量的粒子？

在那种情况下，我可以将相邻单元格中的所有粒子加载到本地内存中，并通过它们迭代计算某些属性。然后我可以将特定值写入当前网格单元中的每个粒子。

对于所有粒子和每次迭代(大部分时间相同)的邻居，这种方法是否有利于运行内核？

此外，number of particles/number of grid cells 的理想比例是多少？ ？

---

我正在尝试重新实现(和修改)CUDA Particles对于 OpenCL 并使用它来查询每个粒子的最近邻居。我创建了以下结构:

• 缓冲区 P保持所有粒子的 3D 位置 ( float3 )

• 缓冲区 Sp存储 int2粒子 ID 对及其空间哈希。 Sp根据哈希排序。 (哈希只是从 3D 到 1D 的简单线性映射——还没有 Z 索引。)

• 缓冲区 L存储 int2缓冲区中特定空间散列的起始和结束位置对 Sp .示例:L[12] = (int2)(0, 50) .

  • L[12].x是具有空间散列 Sp 的 first 粒子的索引(在 12 中) .
  • L[12].y是具有空间散列 Sp 的last 粒子的索引(在 12 中) .

现在我有了所有这些缓冲区，我想遍历 P 中的所有粒子并为每个粒子迭代其最近的邻居。目前我有一个看起来像这样的内核(伪代码):

__kernel process_particles(float3* P, int2* Sp, int2* L, int* Out) {
  size_t gid             = get_global_id(0);
  float3 curr_particle   = P[gid];
  int    processed_value = 0;

  for(int x=-1; x<=1; x++)
    for(int y=-1; y<=1; y++)
      for(int z=-1; z<=1; z++) {

        float3 neigh_position = curr_particle + (float3)(x,y,z)*GRID_CELL_SIDE;

        // ugly boundary checking
        if ( dot(neigh_position<0,        (float3)(1)) +
             dot(neigh_position>BOUNDARY, (float3)(1))   != 0)
             continue;

        int neigh_hash        = spatial_hash( neigh_position );
        int2 particles_range  = L[ neigh_hash ];

        for(int p=particles_range.x; p<particles_range.y; p++)
          processed_value += heavy_computation( P[ Sp[p].y ] );

      }

  Out[gid] = processed_value;
}

该代码的问题在于它很慢。我怀疑非线性 GPU 内存访问(特别是最内层的 P[Sp[p].y] 循环中的 for)导致了缓慢。

我想做的是使用Z-order curve作为空间散列。这样我就可以只有 1 for在查询邻居时循环遍历连续的内存范围。唯一的问题是我不知道开始和结束 Z-index 值应该是什么。

我想要实现的 chalice :

__kernel process_particles(float3* P, int2* Sp, int2* L, int* Out) {
  size_t gid             = get_global_id(0);
  float3 curr_particle   = P[gid];
  int    processed_value = 0;

  // How to accomplish this??
  // `get_neighbors_range()` returns start and end Z-index values
  // representing the start and end near neighbors cells range
  int2 nearest_neighboring_cells_range = get_neighbors_range(curr_particle);
  int first_particle_id = L[ nearest_neighboring_cells_range.x ].x;
  int last_particle_id  = L[ nearest_neighboring_cells_range.y ].y;

  for(int p=first_particle_id; p<=last_particle_id; p++) {
      processed_value += heavy_computation( P[ Sp[p].y ] );
  }

  Out[gid] = processed_value;
}

Answer 1

您应该仔细研究 Morton Code 算法。 Ericsons 实时碰撞检测很好地解释了这一点。

Ericson - Real time Collision detection

这是另一个很好的解释，包括一些测试:

Morton encoding/decoding through bit interleaving: Implementations

Z-Order 算法仅定义坐标的路径，您可以在其中从 2 或 3D 坐标散列为整数。尽管算法在每次迭代中都变得更深入，但您必须自己设置限制。通常停止索引由哨兵表示。让哨兵停止将告诉您粒子放置在哪个级别。因此，您要定义的最大级别将告诉您每个维度的单元格数。例如，最大级别为 6 时，您有 2^6 = 64。您的系统 (3D) 中将有 64x64x64 个单元格。这也意味着您必须使用基于整数的坐标。如果您使用 float ，则必须像 coord.x = 64*float_x 等那样进行转换。

如果您知道系统中有多少个细胞，您就可以定义限制。您是否尝试使用二进制八叉树？

由于粒子在运动(在该 CUDA 示例中)，您应该尝试并行处理粒子数量而不是细胞数量。

如果您想建立最近邻列表，您必须将粒子映射到细胞。这是通过一个表格完成的，该表格随后按细胞分类到粒子。您仍然应该遍历粒子并访问其邻居。

关于您的代码:

The problem with that code is that it's slow. I suspect the nonlinear GPU memory access (particulary P[Sp[p].y] in the inner-most for loop) to be causing the slowness.

记住 Donald Knuth。您应该测量瓶颈的位置。您可以使用 NVCC Profiler 并寻找瓶颈。不确定 OpenCL 有什么分析器。

    // ugly boundary checking
    if ( dot(neigh_position<0,        (float3)(1)) +
         dot(neigh_position>BOUNDARY, (float3)(1))   != 0)
         continue;

我认为你不应该那样分支，当你调用 heavy_computation 时返回零怎么样。不确定，但也许您在这里有某种分支预测。尝试以某种方式删除它。

仅当您没有对粒子数据的写访问权时，在单元上并行运行才是一个好主意，否则您将不得不使用原子。如果您越过粒子范围，您反而会读取对单元格和邻居的访问权限，但您会并行创建总和，并且不会被迫采用某种竞争条件范式。

Also, what is the ideal ratio of number of particles/number of grid cells?

确实取决于您的算法和您所在域内的粒子堆积，但在您的情况下，我会定义与粒子直径相等的单元大小，并且只使用您获得的单元数。

因此，如果您想使用 Z 顺序并实现您的 chalice ，请尝试使用整数坐标并对其进行哈希处理。

也尝试使用更大量的粒子。您应该考虑大约 65000 个像 CUDA 示例使用的粒子，因为这样并行化最有效；正在运行的处理单元被利用(更少的空闲线程)。