流体仿真实战

发表于 分类于 ， 阅读次数：

流体仿真实战

时间：2026/04/09

关键词：Eulerian grid、advection、diffusion、projection、pressure solve、ping-pong、texture object、CUDA
核心目标：从工程角度理解一个实时流体模拟在 GPU 上通常由哪些步骤组成，以及每一步为什么这样设计。

1. 先建立问题模型

实时流体仿真通常不是在“跟踪每一个分子”，而是在求解离散化后的场：

  • 速度场 u(x, y[, z])
  • 压力场 p
  • 密度/烟雾浓度 d
  • 温度、外力等附加场

在图形学实时仿真里，最常见的是 Eulerian 网格法：

  • 空间被划分成规则网格
  • 每个网格单元保存局部状态
  • 关注的是“这个位置上的流体状态”，而不是“某个粒子去哪了”

2. 一帧流体模拟通常做什么

一个经典的 2D/3D 实时烟雾模拟步骤，通常长这样：

  1. 对速度场做平流（advection）
  2. 对密度/温度做平流
  3. 加入外力（force / impulse）
  4. 计算散度（divergence）
  5. 解压强泊松方程（pressure solve）
  6. 用压力梯度做投影（projection）
  7. 得到无散速度场
  8. 渲染密度或其他可视化结果

这条管线里最关键的思想是：

  • 平流负责“把东西搬着走”
  • 压力投影负责“让速度场满足不可压缩约束”

3. 为什么很多实时流体用 Stable Fluids

Jos Stam 的 Stable Fluids 思想之所以经典，是因为它强调：

  • 数值上更稳定
  • 能容忍较大的时间步长
  • 非常适合图形实时场景

它的代价通常是：

  • 会更粘、更耗散
  • 细节可能被抹平

所以工程上通常要在：

  • 稳定性
  • 速度
  • 视觉效果

之间做平衡。

4. 网格、边界与数据布局

4.1 最常见的数据

二维时经常有：

  • vel_x
  • vel_y
  • pressure
  • divergence
  • density

4.2 常见存储方式

对 GPU 更友好的方式通常是：

  • 扁平连续数组
  • 纹理/表面对象
  • ping-pong 双缓冲

例如二维数组常展平成：

1
idx = y * width + x;

4.3 为什么常用 SoA 思路

不同字段通常分开存：

  • 一个数组/纹理存 vel_x
  • 一个数组/纹理存 vel_y
  • 一个数组/纹理存 pressure

这比把所有字段塞进一个大结构更容易：

  • 控制带宽
  • 做批量 kernel
  • 针对单一字段优化访问

5. 平流（Advection）是什么

直觉上，平流就是：

按速度场把某个量“搬运”到下一个时刻。

对于密度场来说，可以理解成：

  • 烟雾随着流速移动

对于速度场本身来说：

  • 速度也会被自身流动搬运

5.1 半拉格朗日（Semi-Lagrangian）平流

实时图形里很常见的做法是“反向追踪”：

  • 对当前网格点 (x, y)
  • 沿速度反方向回溯到上一个时刻的位置
  • 从旧场中插值取样

公式直觉：

1
phi_new(x) = phi_old(x - dt * u(x))

5.2 为什么它适合 GPU

因为每个网格点都能独立算：

  • 回溯位置
  • 插值取样
  • 写回新值

这非常适合 data-parallel kernel。

6. 为什么纹理对象在流体里常见

流体平流特别需要“按浮点坐标插值采样旧场”。
这正是纹理对象很擅长的事情。

6.1 纹理对象的实际优势

  • 支持硬件插值
  • 支持边界寻址模式
  • 访问接口清晰
  • 对二维/三维场很自然

因此在 CUDA 流体项目里，经常会看到：

  • 旧场绑定成 texture
  • 新场写到 output buffer / surface

6.2 一个常见工程模式

  • src 作为只读纹理
  • dst 作为可写 buffer
  • kernel 结束后交换 src/dst

这就是典型的 ping-pong。

7. Ping-Pong 双缓冲为什么几乎是标配

平流、扩散、压力迭代这类步骤，通常不能一边读旧值一边原地覆盖写新值，否则会污染后续计算。

最常见做法：

  • A 存旧场
  • B 存新场
  • 一步计算后交换 A/B

例如：

1
2
3
4
advect:    vel0 -> vel1
swap
project:   vel1 -> vel0
swap

双缓冲的代价是多占一份显存，但好处是：

  • 读写依赖清晰
  • kernel 更容易写对
  • 更适合并行

8. 散度（Divergence）为什么要算

不可压缩流体通常要求：

1
div(u) = 0

但经过加力、平流等步骤后，速度场往往不再满足这个条件。
所以需要先计算散度，看看当前速度场“有多不守恒”。

二维中心差分的直觉形式：

1
div = (u_right - u_left + v_up - v_down) * 0.5 / h

算出来的散度场会成为后面压力方程的右端项。

9. 压力求解（Pressure Solve）是核心瓶颈之一

为了让速度场重新无散，通常要解一个泊松方程：

1
∇²p = div

这里：

  • p 是压力
  • div 是上一阶段算出的散度

9.1 为什么这是难点

因为它不是一个局部一步完成的更新，而是一个全局耦合问题。
实时实现中常用近似迭代法，例如：

  • Jacobi
  • Gauss-Seidel
  • Red-Black Gauss-Seidel
  • Multigrid（更复杂也更强）

9.2 为什么很多教程先用 Jacobi

因为它简单、规则、容易并行：

  • 每次迭代只读旧压力
  • 写新压力
  • 适合 GPU 大规模并行

代价是收敛速度不算快。

10. 投影（Projection）在做什么

解出压力场后，用压力梯度修正速度：

1
u' = u - ∇p

直觉上就是：

  • 从原速度里减掉“造成体积膨胀/压缩”的那部分

投影之后，速度场更接近无散状态，也就更符合不可压缩流体假设。

这一步通常是：

  1. 读取速度
  2. 读取相邻压力
  3. 做梯度差分
  4. 写回新速度

11. 边界条件不能忽略

很多新手代码“能跑但看起来不对”，问题往往在边界。

常见边界处理包括：

  • 固壁边界
  • 开放边界
  • 周期边界
  • 障碍物边界

例如固壁边界时，速度在法向方向通常要满足不穿透约束。

如果边界没处理好，可能会出现：

  • 烟雾穿墙
  • 边缘数值爆炸
  • 压力解不稳定

12. CUDA 上通常如何拆 kernel

一个典型实时流体项目不会把全部逻辑塞进一个 kernel，而是拆成多步：

  • advect_velocity
  • advect_density
  • add_force
  • compute_divergence
  • jacobi_pressure
  • project_velocity
  • render_density

这样做的好处：

  • 每步逻辑清晰
  • 便于调试
  • 便于替换某一步算法

代价是：

  • kernel launch 次数变多
  • 中间结果要反复读写全局内存

所以工程优化通常是在“模块清晰”和“减少中间带宽”之间折中。

13. 一个典型的 CUDA 数据流

可以把一帧想成这样：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
velocity_src --advect--> velocity_dst
swap

density_src --advect--> density_dst
swap

velocity_src --add_force--> velocity_dst
swap

velocity_src --divergence--> divergence

pressure_src --jacobi--> pressure_dst
swap
... 迭代多次 ...

velocity_src + pressure_src --project--> velocity_dst
swap

density_src --render--> framebuffer

其中压力求解阶段往往会迭代很多轮，是核心带宽热点之一。

14. 纹理对象在这里最常见的落点

14.1 平流时

因为需要按浮点坐标回溯并插值取样旧场。

14.2 可视化时

例如把密度场映射到图像。

14.3 有时不用纹理也能做

完全可以手写双线性插值，但纹理对象通常更方便，也更贴近 GPU 图像管线思维。

15. 性能瓶颈通常在哪

15.1 不一定是算术，而是带宽

流体很多步骤都是 stencil / 邻域访问：

  • 读周围多个点
  • 写一个点

因此常见瓶颈是：

  • 全局内存带宽
  • 多次中间场读写

15.2 压力迭代轮数

Jacobi 每多迭代一轮，都意味着：

  • 再扫一遍整个压力场

所以视觉质量和帧率经常直接在这里博弈。

15.3 边界和条件判断

边界格子处理如果分支太多，也会影响 warp 执行效率。

16. 常见优化方向

16.1 减少全局内存往返

例如：

  • 更少的中间缓冲
  • 更合理的 kernel 融合

16.2 用 shared memory 做局部邻域缓存

对于 stencil 类计算，局部 tile 常能减少重复全局加载。

16.3 让访问尽量规则

例如：

  • 扁平数组
  • 合理线程块布局
  • 连续内存访问

16.4 区分“视觉上够用”和“数值上更准”

实时图形里很多时候不是求学术最优，而是：

  • 足够稳定
  • 足够快
  • 视觉上自然

17. 为什么这类项目特别适合做综合练习

流体仿真把很多高性能主题都串起来了：

  • 数据布局
  • 邻域访存
  • ping-pong 双缓冲
  • GPU kernel 设计
  • 纹理对象
  • 迭代求解器
  • 可视化与物理折中

所以它既是一个图形学题目，也是一个很典型的并行计算工程题。

18. 初学实现建议

如果你第一次做，比较稳的路线是：

  1. 先做 2D，不要一开始上 3D
  2. 先做密度平流，再加速度场
  3. 再补散度、压力、投影
  4. 压力求解先用 Jacobi 跑通
  5. 最后再考虑纹理对象、shared memory、性能优化

因为最难的不是“写出一个 kernel”，而是：

  • 整条数据流逻辑正确
  • 每步边界一致
  • 数值上稳定

19. 一页总结

实时流体仿真的工程骨架，可以压缩成下面几句：

  1. 用规则网格存速度、压力和密度场
  2. 用平流把量沿速度场搬运
  3. 用散度和压力解修正速度场
  4. 用投影得到近似不可压缩流体
  5. 用 ping-pong 双缓冲保证每一步读写分离
  6. 用纹理对象和连续布局提高 GPU 访问效率

如果只记一个最核心的工程关系，那就是：

平流负责“搬”，压力投影负责“稳”。

20. 建议继续补充的相关主题

和本篇衔接最紧密的内容：

  1. CUDA texture / surface object 细节
  2. Staggered grid（MAC grid）
  3. Jacobi / Gauss-Seidel / Multigrid 对比
  4. Vorticity confinement 与视觉细节增强
  5. 3D 烟雾与体渲染

21. 参考资料

  1. Jos Stam, Stable Fluids
    https://www.dgp.toronto.edu/public_user/stam/reality/Research/pdf/GDC03.pdf

  2. GPU Gems: Fast Fluid Dynamics Simulation on the GPU
    https://developer.nvidia.com/gpugems/gpugems/part-vi-beyond-triangles/chapter-38-fast-fluid-dynamics-simulation-gpu