物理模拟
SPGird
00 min
2025-1-6
2025-1-6
type
status
date
slug
summary
tags
category
icon
password
SPGrid
引言贡献SPGridArray layoutAddress computation and multi-channel arrays访问模式

引言

计算机图形学研究领域探索了利用流体技术进行空间自适应模拟的可能性，包括自适应笛卡尔网格，自适应四面体网格，以及基于粒子的方法。在相同自由度下，自适应技术能够显著提高视觉效果的丰富度，但其计算成本也显著增加。在统一网格上，由于可优化内存访问模式以提升缓存命中率，同时简化并行性中的域分区，模板计算和稀疏线性代数等常见操作的效率得以大幅提升。
该论文提出了一种名为 SPGrid（Sparse Paged Grid 的简称）的全新数据结构，专为稀疏填充的统一网格设计，支持高效的存储与流处理。SPGrid 在保持统一网格计算效率和编程便利性的同时，其存储空间需求仅略高于实际稀疏网格的数据量。该论文基于此数据结构开发了自适应流体模拟的离散化方法，用稀疏的统一网格金字塔替代了传统依赖树结构与指针的表示方法（见下图）。我们证明了在自适应计算流体力学模拟中，所有主要的计算核心（如对流、离散梯度、散度与拉普拉斯运算）均可在框架中实现，而无需显式的树结构遍历，这得益于高效的统一网格操作内核。此外，还提出了一种轻量化的自适应多重网格预条件 Krylov 求解器，其在收敛性与可伸缩性方面优于基于不完全分解的传统方法。
 
notion image

贡献

我们的核心贡献总结如下：
  1. 提出了一种新颖且轻量化的数据结构（SPGrid），专为处理大型稀疏笛卡尔网格而设计。该结构利用虚拟内存子系统的硬件加速机制，在性能上可与密集统一网格相媲美。同时，SPGrid 支持多通道数据的协同存储，即使不同通道的数据域存在差异。SPGrid 通过 Morton 编码将三维数据映射为线性内存空间，并实现了无需转换为几何坐标即可直接操作线性化（编码后）地址的模板运算。
  1. 将传统基于指针的八叉树样式自适应流体离散化重构为基于 SPGrid 金字塔的形式。此方法将昂贵的计算内核（如对流、梯度、散度与拉普拉斯运算）分解为简单且易于并行化的统一网格内核操作，即使在非梯度八叉树（相邻单元分辨率可能相差超过一级）的情况下依然适用。
  1. 提出了一个自适应多重网格预条件 Krylov 求解器，能够实现与分辨率无关的收敛速率，并在相同的 SPGrid 金字塔中以无矩阵形式实现，兼容其他流体模拟内核。
  1. 开发了一种插值方案，既降低了对流中的耗散效应，又简化了网格拓扑的动态自适应调整。
尽管这些贡献在我们的最终成果中紧密结合，形成了完整的解决方案，但它们也具有模块化与可复用性。例如，窄带水平集应用可以单独使用 SPGrid，而无需采用自适应扩展。同样，在需要外存处理的场景中，可以用 OpenVDB 替代 SPGrid 存储稀疏网格，同时保留金字塔自适应概念。此外，我们的补充材料中还详细说明了如何利用 SPGrid 优化存储和应用不完全 Cholesky 预条件器作为多重网格的替代方案。
 

SPGrid

在不同的应用领域中，稀疏和自适应存储方案（如八叉树、哈希网格、RLE 等）通常会根据具体需求加以选择。对于计算流体力学（CFD）工作负载，则有一些特定的设计要求。首先，顺序访问和模板访问是最常见的使用场景；即使是“随机”访问模式，也往往表现出显著的空间相干性（例如半拉格朗日方法中的对流运算）。此外，迭代求解器（如共轭梯度法）以及直接求解方法或其组件（如不完全 Cholesky 分解）通常受到内存带宽的限制，因此要求数据存储方案在内存带宽和内存占用方面经过优化。这一点在多处理器（或 SIMD）平台上尤为关键，因为如果内存利用效率不高，将严重影响系统性能。此外，流体模拟需要在多个数据通道上运行，包括标量（如压力、密度）和向量值（如速度，以及 Krylov 求解器的辅助量等）。这些通道的几何分布可能不一致，但它们往往占据相似的空间范围。
论文提出了一种针对稀疏均匀网格的数据结构，可应用于解决诸如窄带水平集问题之类的挑战。需要注意的是，原则上可以通过直接分配一个密集的均匀网格并仅使用其稀疏的相关子集来满足上述算法的需求，前提是内存占用不是问题。在这种情况下，模板访问可以得到极大优化，因为网格位置的特定邻居始终位于相对于参考数据点地址的固定内存偏移量中。然而，对于我们研究中的典型网格情况（网格占用率可能低至 0.1%），这种直接分配密集网格的方法显然不切实际。然而，这一思路可以类比于虚拟内存的工作方式：CPU 允许进程访问一个巨大的虚拟内存地址空间，只要实际访问的稀疏子集（以 4KB 页面为单位）足够小，就可以被映射到物理内存中。
我们提出的数据结构结合了从二维或三维网格位置到一维内存跨度的局部保留映射，并利用虚拟内存机制避免了未使用区域的物理分配。这种方法有效减少了存储开销，同时保留了高效的模板访问特性，适用于内存敏感的稀疏网格应用场景。
 

Array layout

我们提出的 SPGrid 数据结构为多维（2D/3D）数组提供了一种抽象表示，用于存储稀疏子集的数据内容。类似于 C++ 的静态数组，每个数组坐标——无论是否填充——都在（虚拟）内存中拥有一个固定的存储位置。与基于指针的八叉树、VDB 或哈希网格不同，这些方法需要在首次访问时动态分配存储，并且在几何坐标与实际内存位置之间维护显式映射。而 SPGrid 则通过保留一个巨大的虚拟内存地址范围，实现了稀疏网格的表示。尽管如此，只有一小部分被访问的数组条目会实际占用物理内存空间。
与传统 C++ 数组采用字典顺序的内存布局不同，SPGrid 采用自定义的映射模式，将 2D/3D 几何坐标映射到线性化的内存偏移量。首先，我们将数组划分为尺寸精心选择的小矩形块（例如，在我们的流体模拟中采用 4×4×4 的块）。然后，块内的数据以字典顺序排列，而块之间按照莫顿曲线的顺序进行布局。莫顿排序是一种局部优化的映射方式，能够确保几何上接近的数组条目在内存中也被存储在相邻位置。因此，如果数组的数据分布在空间上具有较高的局部性，对应的内存布局也会表现出高度聚集性。为了支持这种布局，SPGrid 为整个数组预留了一个虚拟内存地址范围，但仅在需要访问内容时才实际分配物理内存。
现代操作系统支持预留虚拟内存而不立即分配物理内存的功能。在 Linux 中，这一功能可以通过 mmap 系统调用实现。我们在实现中使用了如下调用：
上述调用会立即为数组保留一个巨大的虚拟内存范围（受架构限制，每个进程最大可达 128TB），但不会分配任何物理内存。通过指定 MAP_ANONYMOUS 和 MAP_NORESERVE 参数，我们进一步指示操作系统在该范围内永不启用交换机制，从而确保以下约束：虚拟地址范围内的页面要么为物理内存驻留，要么从未被访问过（如违反此条件，将触发内存不足错误）。
在完成虚拟内存范围的预留后，对该范围内内存地址的访问会导致以下两种情况之一：
  1. 页面故障：如果请求的地址所在的虚拟页面尚未映射，则会触发页面故障，操作系统会分配一个物理页面，在页面表中记录映射，并将物理页面内容初始化为零。
  1. 快速访问：如果请求的地址所在的页面已被映射（即之前已被访问过），地址翻译（包括 TLB 异常处理和页表遍历）将完全由硬件完成，无需操作系统干预。
第二种情况的性能与访问通过其他机制（如 malloc 或 new）分配的内存无异。需要注意的是，为使 SPGrid 的方案生效，虚拟内存的预留必须通过 mmap（或类似的 VirtualAlloc 机制）完成，而不能使用 malloc 或默认的 C++ 操作符 new，因为后者会立即将虚拟地址范围映射到物理内存。
最后，我们假定数组的内存分配始终与 4KB 的页面边界对齐。如果使用 mmap，这一要求会自动得到满足。
 

Address computation and multi-channel arrays

 
notion image
我们要求 SPGrid 支持存储多个数据通道（对应于不同的模拟变量）。SPGrid 并未为每个数据通道单独分配内存，而是选择在块粒度上将多个通道的数据交错存储在内存中。通过选择合适的几何块大小（geometric block size），确保块中所有通道的数据可以精确适配一个 4KB 内存页。例如，上图展示了一个 4KB 内存页，其存储了四个单精度浮点数据通道 (u, v, w, p)。每个通道占用 1KB，对于单精度浮点数（每个数据 4 字节），每通道可以存储 256 个条目。因此，几何块的大小可以设置为 8×8×4（三维空间中）。
为了实现快速访问，我们为每个数据通道预先计算了一个基指针，该指针指向该通道中第一个条目的存储位置（即对应零笛卡尔坐标的条目）。任意通道条目(i,j,k)的字节偏移量可通过基指针计算如下图所示：
notion image
偏移量由 64 位组成，其中：
  • 低 12 位：用于表示在单个 4KB 页面内的字节偏移；
    • 最低两位：表示单个数据四字节
    • 接下来的8位：在单个4KB page（1024字节）里面的字节偏移
    • 剩余最高2位：有四个数据通道
  • 高 52 位：用于表示虚拟页面索引。
 
  1. 交错模式的确定性： 数据交错模式在编译时已知，因为它仅取决于通道的数量和大小。每个坐标(i,j,k)的位顺序在交错过程中保持不变，因此可以通过存储一个 64 位掩码（每个坐标一个）来编码每个坐标的位分布情况。
  1. 对齐与跳跃： 在位模式中插入的零不仅用于根据数据大小对齐，还用于在跳过其他通道后移动到下一个块。
  1. 独立于数组的最大尺寸： 位交错模式不需要关注数组的最大尺寸。在组装 64 位偏移量时，我们构造扩展掩码以尽可能多地利用坐标位，即使实际数组尺寸较小。
位交错计算的复杂度高于字典序访问。在静态 C++ 数组中，字典序偏移量的计算仅需 2 次加法和 2 次乘法，而位交错转换的成本则更高。然而，在 Haswell 架构的 CPU 上，由于硬件支持（如 pdep 指令），这一成本显著降低。对于三维坐标的转换，所需时钟周期约为 16。而在 Ivy Bridge 架构中，优化后的实现需要 85 个时钟周期。
需要注意的是，在完全随机访问、缺乏空间局部性的情况下，转换成本通常并非性能瓶颈，因为缓存/TLB 未命中的代价更大。然而，在流式访问（如顺序访问或模板访问）中，这种成本是不可接受的。为此，我们开发了替代的高效方法以进一步降低开销。

访问模式

在计算流体动力学（CFD）工作负载中，完全随机访问非常罕见。这类应用的大部分计算负载集中在流式内核中，即对数据进行顺序或模板访问。以下描述了我们用于高效执行这些内核的方法，使其性能接近密集数组。
顺序访问
notion image
顺序访问指的是一种逐个流式处理一个或多个数据通道内容的模式（例如按照上图所示的遍历顺序），并仅在具有相同数组索引的数据条目之间执行操作。例如，按常数缩放数据通道、将一个通道复制到另一个通道或 BLAS 的 Saxpy 子程序都是顺序访问的例子。归约操作（如点积、最大值/最小值计算）也符合这种模式，因为它们遵循相同的访问模式。
在每个几何块中，执行顺序内核是相对简单的，因为每个通道的数据在内存中是按顺序排列的，并且对齐的，这允许编译器自动发出 SIMD 指令。此外，TLB 缺失得以减少，因为一个块的所有通道都位于同一内存页中。
唯一的问题是如何迭代稀疏的已填充块集合。为此，我们使用了两个加速结构：
  1. 位图记录：在初始化稀疏数组的拓扑结构时，我们使用位图记录所有被访问过的几何块。每 4KB 内存只需要 1 位，因此位图的大小仅占虚拟内存空间的极小部分。例如，一个存储 16 个浮点通道、网格分辨率为 2K×2K×2K 的 SPGrid（对应 512GB 的“虚拟”大小和 5GB 的实际数据负载，假设 1% 的稀疏性），其位图大小仅为 16MB，占用数据大小的 0.3%。
  1. 线性化块偏移数组：
      2. 为了高效遍历稀疏块集合，我们构建了一个线性化 64 位偏移量数组，记录每个填充块第一个条目的位置。这个偏移量数组会在每次改变 SPGrid 稀疏模式的操作后批量更新。这些偏移量并不特定于通道，因为它们指向每个块的第一个条目，无论通道如何，其较低的 12 位始终为零。这种方法甚至适用于具有不同数据大小的通道（例如浮点和双精度）。
在某些情况下，可以替代性地存储每个已占用块的几何坐标；这种情况下，我们只需为块的根条目支付一次转换成本，而剩余块条目可以通过顺序访问轻松摊销掉该成本。不过，我们倾向于避免这种方法，尽可能在“线性化偏移”空间中操作，以提高大多数内核的效率。？
模板访问
从性能角度来看，最重要的访问模式之一是对稀疏索引集的顺序迭代，并在每个位置应用几何模板（例如 7 点拉普拉斯算子）。如前所述，顺序迭代模板应用中心的位置可以通过块偏移数组来实现。主要的挑战是计算 几何邻居 的内存位置，这些邻居需要在模板评估中被拉入。虽然可以实时地将几何邻居坐标转换为内存偏移量，但这对于三维的 7 点拉普拉斯算子来说，每个邻居需要 6 次转换；即使在 Haswell 处理器上，这也会限制我们达到大约 5% 的处理吞吐量，而同样的操作在字典序排列的密集数组上则能达到更高的吞吐量（在那里，模板邻居的偏移量会转化为内存中的常量步长）。
我们提出了一种方法，可以在线性化的偏移空间中本地计算模板邻居的位置，而无需显式地转换为几何坐标。例如，考虑一个模板，它包括一个与模板参考点偏移量为 (−1,2,1)(-1, 2, 1) 的几何邻居。我们通过将整数坐标转换为打包的 64 位偏移量（如图 5 所示）来计算这个相对坐标的扁平化版本。负值采用二进制补码表示，并按常规方式进行位多路复用。我们的关键观察是，直接将这个打包的偏移量加到一个基准（同样是打包的）扁平化索引上，是一种非常高效的操作。更具体地说，给定一个位置 和一个偏移量 的打包翻译，可以直接评估加法结果 的打包翻译，而无需转换为几何（解包的）坐标。
作为第一步，我们描述了一个更简单的操作，我们将其称为 MaskedAdd，该操作根据位掩码从其两个 64 位输入 II 和 JJ 中提取相同的位集合，对提取的值进行正常加法运算，然后根据相同的掩码将加法结果的位分布到 64 位的最终值中：
本质上，MaskedAdd 在每个掩码为零的位位置上将零插入到 II 中，并将一插入到 JJ 中相同的位位置。这个在掩码指定的位位置之间插入的 1 位串，作为传播进位到下一个有效位位置的载体。三个 MaskedAdd 调用可以结合起来，添加两个打包的 3D 坐标，使用三个相应的位掩码：
我们将几何邻居偏移量预先翻译为相应的打包偏移量。然后，通过执行 PackedAdd 函数，我们可以非常高效地获取邻居并应用模板；在 Haswell 处理器上，这个过程大约需要 5 个时钟周期，而在 Ivy Bridge 系统上则需要 7 个时钟周期（如果 II 或 JJ 是编译时常量，则可以进一步节省开销，这在实际应用中经常发生）。值得注意的是，生成的 PackedAdd 汇编代码使用了 21 条指令。然而，这些指令在这些 CPU 的多个功能单元上分布得非常好，使得超标量执行成为可能。
 
上一篇
Shallow Water Equations
下一篇
C++ 编译器导致的优化错误
Catalog