GPU Gems 1 · Chapter 38. Fast Fluid Dynamics Simulation on the GPU

GPU Gems 1 第 38 章的主线是一个很实用的转换：把连续的流体方程离散到二维网格，再把每个网格 cell 交给 GPU 的 fragment pipeline。原章建立在早期浮点纹理、render-to-texture 和双缓冲之上；今天可以把它理解为 GPU texture 或 buffer 上的一组字段，以及由多个 shader pass 组成的 solver，而不是把历史 API 当成唯一实现。

1. 稳定流体：先确定一个可执行的时间步

stable fluids 是 Mark Harris 等人在本章中采用的思路，源自 Stam 的 Stable Fluids。它的取舍很明确：通过半拉格朗日平流获得无条件稳定的时间推进，再用投影维持不可压缩条件；代价是数值耗散、有限网格分辨率和需要多轮迭代的压力求解。

我们不直接在 GPU 上保存一个抽象的“流体对象”，而是保存网格上的字段：速度向量决定染料怎么走，密度代表烟或墨水，压力帮助速度满足约束，障碍物 mask 描述墙的位置。每个时间步按约定顺序消费上一阶段的稳定结果。

最小的 solver 可以写成下面的概念流程。它不是某个 API 的完整实现，却把每个 pass 的输入输出边界说清楚了：

addSource(density, source)
addSource(velocity, force)
diffuse(velocity, viscosity, dt)
advect(velocity, velocity)
advect(density, density, velocity)
project(velocity, pressure, divergence)

源项和边界条件会影响每一个阶段。实际工程中也可能先投影速度再平流密度，或为不同字段使用不同的扩散设置；重要的是每一步都从已完成的输入读取，并且明确哪一个字段在何时成为下一阶段的状态。

2. 把连续场装进纹理

semi-Lagrangian advection 的第一步不是写 shader，而是决定数据布局。二维网格的每个 cell 对应一个 texel，也对应一次 fragment invocation。速度可以把两个分量放在浮点纹理的 R、G 通道；密度、压力和散度可以使用独立纹理，或在带宽与精度允许时打包到不同通道。

渲染一个覆盖计算域的 quad，就能让 fragment 坐标遍历所有 cell。对于每个输出 texel，shader 读取同一坐标的输入字段和邻居字段，应用离散方程，再写入另一张渲染目标。绝不能把当前输入纹理同时绑定为采样源和写入目标：即使某些驱动看起来“能跑”，结果也依赖未定义的读写顺序。

这里的 <Term def="在两个等价的 GPU 表面之间交替读写、避免一个 pass 同时读写同一目标的存储模式">ping-pong surface</Term> 是最小的同步协议：本轮从 A 读、向 B 写；pass 完成后交换 A 和 B。压力的 Jacobi 迭代、速度扩散和平流都可以复用这个协议。字段是否需要单独 surface，要根据格式、带宽和同时读取的通道数做决定，不能为了少几个对象而牺牲清晰的数据契约。

3. 源项、平流与数值耗散

加入源项是最容易理解的一步：用户在鼠标位置注入密度，或把一个局部 force 加到速度。把画笔栅格化到 source texture 后，addSource 只需逐 texel 相加，并对边界区域应用规则。源项应当乘以时间步长或显式强度，避免帧率变化改变注入量。

平流则采用回溯：对当前 cell 位置 x，从速度场读取 u(x)，计算 x_back = x - dt * u(x)，再到上一张字段纹理采样 field(x_back)。双线性采样让结果比最近邻更平滑；clamp 规则保证回溯位置不会访问网格外的无效内存。下面的伪代码刻意省略 API 细节，只突出依赖方向：

float2 backtrace(float2 cell, Texture velocity, float dt) {
  float2 u = sample(velocity, cell);
  return clamp(cell - dt * u, domainMin, domainMax);
}
 
float advectCell(float2 cell, Texture velocity, Texture field, float dt) {
  float2 previousCell = backtrace(cell, velocity, dt);
  return bilinearSample(field, previousCell);
}

半拉格朗日方法的稳定性来自“沿流线回到过去取值”，但每次插值都会抹平高频细节，因此烟雾会逐渐变淡。增大网格分辨率、减少不必要的重采样、加入受控的 vorticity confinement 或提高源强，都可能改善视觉效果；这些都不能替代对时间步和边界的检查。

4. 扩散：用 Jacobi 迭代换取可并行的近似解

扩散和压力求解都可以归结为邻居线性方程。GPU 不适合在一个 fragment 中顺序解完整的全局线性系统，因此使用 <Term def="反复用邻居的旧值更新当前网格值、每轮在两张表面之间交换的线性系统近似求解方法">Jacobi iteration</Term>：每一轮每个 cell 独立读取邻居旧值，计算一个新值，再写入下一张 surface。

以四邻居 stencil 为例，一轮的结构是：

for each cell in parallel:
    neighbors = read(input, left, right, up, down)
    next[cell] = (source[cell] + neighbors.sum * alpha) / beta
swap(input, next)
repeat a fixed number of iterations

它很适合 fragment pass，因为所有输出只依赖上一轮的稳定纹理；它也很适合教学和调参，因为迭代次数直接暴露了精度与成本的关系。轮数太少，压力残差或扩散误差明显；轮数太多，带宽和 pass 调度可能压低帧率。固定轮数通常比等待一个不可预测的全局收敛条件更容易满足交互预算，但应在开发阶段记录残差。

扩散速度时要处理障碍物邻居：墙内的 texel 不应被当成普通流体 cell。密度扩散也可能使用不同的边界规则。把“所有字段一套邻居公式”当成捷径，往往会在墙边产生泄漏。

5. 让速度满足不可压缩约束

源项、扩散和平流之后，速度场通常不是不可压缩的。先用 <Term def="速度场的局部膨胀率，通常通过相邻速度分量的差分计算">divergence</Term> 估计每个 cell 的散度：正值表示局部向外膨胀，负值表示汇聚。然后用 Jacobi 迭代近似求压力场，最后执行 <Term def="先计算散度、再求压力并减去压力梯度，使速度场回到近似无散状态的修正阶段">projection</Term>。

概念上可以记成三步：

1. divergence = ∇ · velocity：把速度场违反约束的程度写到纹理。
2. pressure = solve(divergence)：用多轮邻居更新近似解压力 Poisson 方程。
3. velocity -= gradient(pressure)：从速度中减掉压力梯度，让相邻 cell 的流量尽量守恒。

投影不是“让画面好看”的可选后处理。没有它，源项会不断制造体积误差；有了它，迭代次数不足时仍可能留下残差。可以在 debug view 中显示散度热图，或用一个 reduction 统计最大绝对散度，把视觉结果和数值证据连接起来。

6. 边界、障碍物与 GPU 读写纪律

边界条件决定了流体域的物理含义。墙面常用 no-slip 约束：靠墙的速度分量被置零或按法向反射；密度可以贴着墙流动，也可以按问题定义被吸收。障碍物 mask 需要被扩散、平流、压力和投影共同读取，否则某一阶段把量带过墙，下一阶段再严谨也无法恢复。

回溯位置落到域外时要 clamp 或使用专门的边界 texel；访问一个未定义的纹理坐标可能产生随机密度。压力 stencil 还要区分流体邻居和墙邻居，避免把墙内的压力当成普通未知量。建议给每个 pass 写一条“允许读取哪些字段、允许写哪一个 surface、遇到 mask 怎样处理”的数据契约，再开始优化布局。

历史 GPU 实现用 render-to-texture 和浮点纹理把每轮结果留在显存中。现代实现可以换成计算着色器、storage image 或 GPU buffer，但两个约束仍然一样：读写依赖要用显式 barrier 或分离资源表达，CPU readback 要尽量推迟。把整张速度场每帧读回 CPU 做 debug，会让一个本来适合 GPU 的 pipeline 变成同步瓶颈；更好的做法是只读回小摘要，或使用 GPU 端可视化。

7. 交互实验：把质量旋钮换成预算

实验把一个二维流体 solver 的主要工程旋钮放在一起：平流阶段强调回溯与耗散，投影阶段强调散度与迭代；分辨率平方增长会放大每次纹理读写，Jacobi 迭代会线性增加 pass；障碍物提高了边界正确性但也需要额外 mask 处理；关闭 GPU resident 会显式加入 readback 风险。

可以从 256 × 256、16 次迭代和启用障碍物开始，再逐项改变一个变量。若只提高分辨率，先看带宽和帧时间；若只提高迭代，先看散度热图；若降低黏性，观察细节是否保留但稳定性是否变差。数值是示意，发布前仍需用目标 GPU、真实 shader 和 GPU profiler 校准。

8. 复用清单：从算法到发布实现

把本章方法迁移到一个新的实时模拟时，可以按以下顺序复核：

1. 场的定义：列出每个 cell 的字段、格式、精度和是否需要独立纹理。
2. 时间步顺序：明确源项、扩散、平流、散度、压力、投影和边界的先后关系。
3. 读写契约：每个 pass 只有一个写目标；跨 pass 的依赖通过 ping-pong 和 barrier 表达。
4. 迭代预算：记录 Jacobi 轮数、分辨率、纹理带宽和残差，不用“看起来稳定”替代指标。
5. 边界测试：用静态墙、窄通道和域外回溯测试速度与密度是否穿透。
6. 端到端检查：确认字段是否保持 GPU resident，readback、同步、展示和输入注入是否共同满足帧预算。

小结

• GPU Gems 1 第 38 章把稳定二维流体拆成源项、扩散、平流和投影，每个阶段都能映射为纹理上的并行 pass。
• 半拉格朗日回溯带来很好的稳定性，但插值会造成耗散；Jacobi 迭代能并行近似求解，却需要用 pass 预算换精度。
• 速度、密度、压力和散度应有清晰的数据布局；ping-pong surface 保证输入稳定，不能在同一目标上读写。
• 投影通过散度、压力和压力梯度把速度拉回近似不可压缩状态；边界 mask、clamp 和墙面条件必须贯穿整个 solver。
• 分辨率、迭代次数、障碍物和 GPU resident 策略共同决定交互质量；应结合残差、带宽、同步和视觉巡检做发布判断。