CUDA粒子中最近的邻居

Question

编辑2: 请看一下TLDR的这个交叉路口.

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

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

这种方法是否有利于为所有粒子运行内核以及每次迭代(大多数时间是相同的)邻居？

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

---

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

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

• 缓冲Sp存储int2粒子id对及其空间哈希值.Sp根据哈希值排序.(哈希只是从3D到1D的简单线性映射 - 还没有Z-indexing.)

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

  • L[12].x是具有空间散列Sp的第一个粒子的索引(in )12.
  • L[12].y是具有空间散列Sp的最后一个粒子的索引(in )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顺序曲线作为空间散列.这样for,在查询邻居时,我只能有一个循环遍历连续的内存范围.唯一的问题是我不知道应该是什么开始和停止Z-index值.

我想要实现的圣杯:

__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

您应该仔细研究莫顿代码算法。爱立信的实时碰撞检测很好地解释了这一点。

Ericson - 实时碰撞检测

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

通过位交织进行莫顿编码/解码：实现

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

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

由于粒子是运动的（在 CUDA 示例中），您应该尝试并行化粒子数量而不是单元数量。

如果你想建立最近邻居的列表，你必须将粒子映射到单元格。这是通过一个表格完成的，该表格随后按细胞到粒子进行排序。尽管如此，您仍然应该迭代粒子并访问其邻居。

关于您的代码：

该代码的问题是速度很慢。我怀疑非线性 GPU 内存访问（特别是最里面的 for 循环中的 P[Sp[p].y]）导致速度缓慢。

请记住唐纳德·高德纳 (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？不确定，但也许你在这里有某种分支预测。尝试以某种方式删除它。

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

另外，粒子数/网格单元数的理想比率是多少？

实际上取决于您的算法和域内的颗粒堆积，但在您的情况下，我将定义等于颗粒直径的单元尺寸，并仅使用您获得的单元数量。

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

还可以尝试使用更多的颗粒。您应该考虑使用大约 65000 个粒子（例如 CUDA 示例），因为这样并行化最有效；正在运行的处理单元被利用（更少的空闲线程）。