GPU 上的颗粒

颗粒内核计算密集，并且原则上可以与流体推进的部分异步处理。 MFIX-Exa 中颗粒方法的核心组件是：

  • 邻居列表构建

  • 颗粒-颗粒碰撞

  • 颗粒-壁面碰撞

在这些操作中，邻居列表的构建需要最为注意。 邻居列表是预先计算的列表，列出给定颗粒在接下来的 n 个时间步长中可以与之交互的所有邻居。 邻居列表通常通过按交互距离对颗粒进行分类来构建，然后仅对邻近的分类中的颗粒执行 N2 距离检查。详细来说，邻居列表算法的 CPU 版本如下：

  • 对于每个级别上的每个块，遍历颗粒，识别其所属的分类。

  • 将颗粒添加到“拥有”它的单元的链表中。

  • 对于每个单元，遍历所有颗粒，然后遍历所有可能的碰撞伙伴在相邻单元中。

  • 如果碰撞伙伴足够接近，将其添加到该颗粒的邻居列表中。

为了将该算法移植到 GPU，我们使用了作为 CUDA 工具包一部分分发的并行算法库 Thrust。Thrust 提供并行排序、搜索和前缀求和算法，这些算法在移植颗粒算法时尤其有用。为了在 GPU 上构建邻居列表，我们遵循由颗粒协同设计中心开发的 Canaba 所使用的基本方法：

  • 使用并行计数排序，通过分类对每个网格上的颗粒进行排序。我们使用 Thrust 的 exclusive_scan 函数实现排序的前缀和阶段，其余部分使用手写内核实现。此步骤实际上不涉及重新排列颗粒数据，而是计算一个置换，可以在不实际重新排序的情况下将颗粒按顺序排列。

  • 一旦颗粒按分类排序，我们可以遍历相邻分类中的颗粒。我们对颗粒进行两次遍历。首先，我们启动一个内核来计算每个颗粒的碰撞伙伴数量。

  • 然后，我们对这些数字求和并为我们的邻居列表分配空间。

  • 最后，我们再次遍历颗粒，将它们放入适当位置的列表中。

请注意，我们构建了一个完整的邻居列表，这意味着如果颗粒 i 出现在颗粒 j 的列表中，那么颗粒 j 也会出现在颗粒 i 的列表中。这简化了使用这些列表时的力计算步骤，因为可以在不使用原子操作的情况下更新给定颗粒的力和力矩。

最终的网格邻居列表数据结构由两个数组组成。首先，我们有一个邻居列表本身，存储为一个大的、一维的颗粒索引数组。然后，我们有一个 offsets 数组，存储对于每个颗粒在邻居列表数组中的位置。这种数据结构的细节已隐藏在一个迭代器中，因此用户代码可以如下所示：

// 现在我们遍历邻居列表并计算力
 AMREX_FOR_1D ( np, i,
 {
     ParticleType& p1 = pstruct[i];
     p1.rdata(PIdx::ax) = 0.0;
     p1.rdata(PIdx::ay) = 0.0;
     p1.rdata(PIdx::az) = 0.0;

     for (const auto& p2 : nbor_data.getNeighbors(i))
     {
         Real dx = p1.pos(0) - p2.pos(0);
         Real dy = p1.pos(1) - p2.pos(1);
         Real dz = p1.pos(2) - p2.pos(2);

         ...
     }

请注意，由于我们使用托管内存来存储颗粒数据和邻居列表，上述代码在为 CPU 或 GPU 编译时都能工作。

上述算法处理单个网格上颗粒的邻居列表构建。 当使用域分解时，还必须制作相邻网格上颗粒的副本，可能需要为与其他进程关联的网格执行必要的 MPI 通信。 计算哪些颗粒需要作为虚拟颗粒传输到哪个网格的例程 fillNeighbors，以及更新已作为虚拟颗粒传输的颗粒最新数据的例程 updateNeighbors，也已使用类似于 AMReX 的 Redistribute 例程的方法迁移到 GPU。

注意，当颗粒密集时，只有一小部分虚拟颗粒是网格内颗粒的邻居。AMReX 提供了一个函数 selectActualNeighbors 用于过滤掉不用于构建邻居列表的虚拟颗粒。因此，后续对 updateNeighbors 的调用可以避免传输这些未使用的虚拟颗粒，这在一些测试中显著降低了通信成本。要使用此优化，请在 MFIX-Exa 的输入中将 particles.reduceGhostParticles 设置为 true。

一旦构建了邻居列表，就可以轻松处理与颗粒和壁面的碰撞。

MFIX-Exa 目前运行在仅 CPU 和混合 CPU/GPU 架构上。