无网格法 · SPH 理论与实现

SPH 邻域搜索、Verlet 缓冲与并行一致性

从全对全搜索的复杂度出发，推导均匀网格与空间哈希的代价，讲清 Verlet 邻居表的重建判据、溢出与漏邻居排查，并给出一个可复现的邻居集合一致性算例。

源码知识库生成器更新于 2026/9/22

SPH邻域搜索空间哈希Verlet并行计算进阶

SPH 邻域搜索、Verlet 缓冲与并行一致性

邻域搜索是 SPH 每一步都要执行的核心操作：对每个粒子找出支持半径内的邻居，供核求和、梯度近似与成对力计算使用。它既是性能瓶颈，也直接决定结果的正确性——漏掉一个邻居，密度与内力就会悄悄出错。本文讲清均匀网格、空间哈希与 Verlet 邻居表的复杂度与更新条件，并给出溢出、漏邻居和并行非确定性的排查方法。

1. 结论与适用场景

结论先行：直接全对全搜索是 $O(N^2)$ ，不可扩展；基于空间划分的搜索在近均匀密度下接近 $O(N)$ ；缓存邻居表必须显式处理溢出与失效。

适用场景按规模与硬件划分。粒子数在数千以内、追求实现简单时，可用暴力两两搜索；上万到百万级的 CPU 计算，用均匀网格或空间哈希加链表；GPU 大规模计算，用固定容量桶或排序后的紧凑邻居数组，并优先保证显式溢出报告。当密度高度非均匀（如射流破碎、喷雾）时，固定网格会出现单元过空或过满，可改用多级网格或按粒子数动态调整单元尺寸。

需要强调的是，邻域搜索的正确性不亚于核函数本身：求和式密度、压力加速度、黏性项都建立在"邻域完整"的假设上。任何静默截断或漏邻居，都会改变物理量却不报错，是最难排查的一类问题。

另一条实用判据是数据布局。把粒子按单元排序、使同一单元的索引连续存放（空间排序），能让 GPU 访存更合并、减少随机读取；代价是每步或每若干步要重排粒子及其所有属性。是否排序，取决于访存带宽与重排成本的权衡。对交互式或小规模模拟，简单桶链表的收益已足够；对大规模 GPU 模拟，空间排序往往是必选项。

2. 核近似与守恒离散

邻域搜索服务于核近似：所有核求和 $\sum_j$ 的 $j$ 都取自支持半径 $r_c=\kappa h$ 内的邻居集合。因此搜索半径必须完整覆盖核的支持域，少一点就会截断核尾部、破坏归一化与一致性；多一点则浪费算力。二者之间存在一个守恒与性能的平衡点。

守恒性还要求邻域搜索"距离对称"：若粒子 $i$ 把 $j$ 当邻居，则 $j$ 也应把 $i$ 当邻居。在成对力实现中，通常只遍历 $i<j$ 的粒子对并各给两个粒子加相反的力，这天然保证两粒子的邻居关系一致。但若用固定容量邻居表，两个方向的容量可能不同，一边溢出而另一边没有，就会出现"力只加了一半"的非物理结果。因此并行与缓存实现必须保证对称性，或在采用非对称列表时做明确的一致性校正。

总之，邻域搜索不是孤立的加速技巧，而是保持核求和与成对力正确性的基础设施，必须与核函数、边界处理一起验证。

边界粒子也应参与邻居列表。壁面通常用若干层固定或镜像粒子表示，它们必须被纳入搜索范围，否则近壁流体粒子一侧缺邻居，会出现贴壁异常。边界粒子自身一般不移动，其邻居表可跨步复用，从而降低总体开销，前提是要处理好它们与流动粒子的双向可见性。

3. 核心公式与推导

暴力两两搜索需要检查的粒子对数为

\frac{N(N-1)}{2}

随粒子数平方增长， $N=10^6$ 时达到 $5\times10^{11}$ 量级，完全不可行。空间划分把搜索限制在邻近单元内。设单元边长 $L_{cell}$ ，要求

L_{cell}\geq r_c

则只需检查当前单元及相邻的 $3^D-1$ 个单元。内部粒子的期望邻居数为

N_{nb}\approx \frac{4}{3}\pi\left(\frac{r_c}{\Delta p}\right)^3

在近均匀密度下，各单元粒子数近似恒定，构建与查询的总代价降为 $O(N)$ 。空间哈希把整数单元坐标映射到桶，常用形式为

H = (c_x p_1 \oplus c_y p_2 \oplus c_z p_3)\bmod N_b

其中 $c_x,c_y,c_z$ 为整数单元坐标， $p_1,p_2,p_3$ 为互质大素数， $N_b$ 为桶数。哈希值相同不代表空间单元相同，桶内仍需保存并核对真实单元坐标，否则会把远处的粒子误当邻居。

Verlet 邻居表用略大的搜索半径预先缓存邻居，在若干步内复用：

r_s = r_c + r_{skin}

求力时仍按 $r_c$ 筛选。重建判据为任一粒子相对建表位置的位移超过 $r_{skin}/2$ ：

\max_i\left|\mathbf{x}_i(t)-\mathbf{x}_i(t_{build})\right| > \frac{r_{skin}}{2}

$r_{skin}$ 越大则重建越稀疏、缓存浪费越多；越小则重建越频繁。经验上取 $r_{skin}\approx 0.15\,r_c\sim 0.3\,r_c$ 可获得较好折中。

单元占据数给出候选数量的估计。单元边长取 $L_{cell}\approx r_c$ 时，每个单元的粒子数约为 $(L_{cell}/\Delta p)^D$ ，一次查询要检查 $3^D$ 个单元，因此候选数大致是期望邻居数的若干倍，随后再按 $r_c$ 精确筛选。这个过滤步骤不可省略，否则核求和会把远邻计入。

4. 数值实现与参数（h/粒子间距、CFL、人工黏性）

单元尺寸与支持半径： $L_{cell}$ 略大于 $r_c$ 最省内存；若小于 $r_c$ ，单次查询需遍历更多单元，重复判定成本上升。GPU 上常把单元尺寸设为 $r_c$ 的整数倍并做批量归约。

粒子间距与平滑长度：搜索半径由 $r_c=\kappa h=\kappa\eta\Delta p$ 决定，因此 $\Delta p$ 、 $\eta$ 、 $\kappa$ 任一改变都要重新评估网格与邻居表容量。

CFL 与重建频率：时间步越大，粒子每步位移越大，Verlet 表失效越快、重建越频繁。若把 $\Delta t$ 调大，需同步检查重建次数是否抵消了收益。

人工黏性与截断：高黏性或强耗散会让粒子分布更规则、邻居数更稳定；反之，低黏性高剪切区容易产生空洞与团聚，邻居数波动大，应调高容量上限或改用动态结构。

容量与溢出：固定容量表必须记录每个粒子的实际邻居数并报告溢出，绝不能静默截断；稳健做法是把容量设为期望邻居数的一点五到两倍，并对溢出粒子单独告警。

并行策略：把粒子按单元划分给线程，使同一单元的求和对同一批数据反复读取，可提高缓存命中。GPU 上常用两种方案：按粒子并行，每线程自行遍历邻域，实现简单但访存随机；按单元并行，先归约到共享内存再做配对，访存更合并但同步开销更大。选择取决于粒子密度与硬件。

内存与容量规划：邻居表总容量为粒子数与每粒子容量的乘积，百万级粒子要精确预算显存；宁可把容量略调高并保留溢出告警，也不要静默丢邻居，否则正确性无法保证。

5. 可复现示例（伪代码/小算例）

用一小组随机粒子，把加速搜索与暴力全对全结果逐粒子比较，验证邻居集合完全一致。

输入: 粒子位置 x[N] (3D), 支持半径 rc, 单元边长 Lcell >= rc
# 暴力基准
neighbors_brute[i] = {j : |x[i]-x[j]| <= rc, j != i}
# 均匀网格搜索
grid = {}                        # 单元坐标 -> 粒子索引列表
for i in range(N):
    c = floor(x[i] / Lcell)
    grid[c].append(i)
neighbors_grid = []
for i in range(N):
    c = floor(x[i] / Lcell)
    cand = []
    for d in 相邻27单元(含自身):
        cand += grid.get(c + d, [])
    keep = {j for j in cand if j != i and |x[i]-x[j]| <= rc}
    neighbors_grid.append(keep)
# 对称性检查
for i in range(N):
    for j in neighbors_grid[i]:
        assert i in neighbors_grid[j], "邻域不对称"
# 一致性检查（排序成集合后比较）
assert 排序(neighbors_grid) == 排序(neighbors_brute)
print("邻居数与最大最小比较通过")

验证步骤：其一，与暴力结果完全一致；其二，验证周期边界、负坐标、域角点，以及恰好位于 $r_c$ 上的粒子（用 $\leq$ 判定）；其三，记录最小、平均、最大邻居数与零邻居粒子数；其四，让粒子跨越单元边界与 Verlet 阈值，确认既无漏邻居也无重复计数。对同一算例，还应把上述检查重复到多个时间步，而非只看初始时刻。

6. 常见坑与排查

静默溢出：GPU 固定容量表满时截断邻居，密度与内力被改变却无报错。对策是显式报告溢出，并把容量设为期望邻居数的一点五到两倍。
哈希碰撞未核对单元坐标：不同单元哈希到同一桶，直接使用会让远处粒子被误当邻居。
Verlet 表失效判据只看局部：只跟踪部分粒子或忽略边界运动会导致漏邻居，应以全场最大位移为准。
周期边界未做坐标折叠：跨边界的粒子对搜索不到，界面处出现空洞。
单元边长小于 $r_c$ ：需检查更多单元，且容易重复计数同一粒子对。
并行归约顺序不一致：浮点加法顺序改变导致末位差异；需要严格复现时固定粒子排序与归约顺序，否则应比较物理统计量。
零邻居粒子未告警：孤立粒子在核求和时会得到异常小的密度，应显式统计并排查。
每步都重建 Verlet 表：若重建判据过于保守，缓存便失去意义，性能退化到每步全搜。
按粒子并行又直接写共享数组：会产生竞写、丢失更新，需要原子操作或分块归约。
把哈希表当成精确空间结构：哈希只提供候选，必须回算距离与单元坐标；否则碰撞会静默改变结果。

7. 检查清单与参考

[ ] 加速搜索与暴力结果逐粒子一致（小规模数据集）。
[ ] 周期边界、负坐标、域角点与恰好位于 $r_c$ 的粒子均通过。
[ ] 邻居关系对称： $j\in\mathcal{N}_i$ 当且仅当 $i\in\mathcal{N}_j$ 。
[ ] 记录最小/平均/最大邻居数与零邻居粒子数，容量不溢出。
[ ] Verlet 表重建判据覆盖全场最大位移， $r_{skin}$ 取值合理。
[ ] 不同线程数或 GPU 归约顺序下，物理统计量在容差内一致。
[ ] 粒子跨越单元与缓存阈值时无漏邻居、无重复。

参考资料：

Monaghan J.J., Smoothed Particle Hydrodynamics, Annual Review of Astronomy and Astrophysics, 1992.
Domínguez J.M. et al., Neighbour lists in smoothed particle hydrodynamics, International Journal for Numerical Methods in Fluids, 2011.
Ihmsen M. et al., Parallelized SPH for interactive applications, Computers and Graphics, 2011.
Liu M.B., Liu G.R., Smoothed Particle Hydrodynamics (SPH): An Overview and Recent Developments, Archives of Computational Methods in Engineering, 2010.