烟雾会移动就算流体正确吗？从平流到压力投影拆解实时仿真

时间：2026/04/09

一个二维烟雾程序可能已经“看起来能跑”：密度团会沿速度方向移动，也能渲染出旋涡。但运行几十帧后，烟雾越来越糊；靠近边界时穿墙；投影步骤做了许多次 Jacobi 迭代，速度场的散度仍然不小。

画面会动，只能证明数组在变化，不能证明离散方程、边界条件和缓冲区依赖正确。实时流体仿真同时包含三类问题：连续方程如何拆成离散步骤、每一步怎样验证、这些 stencil 计算如何高效映射到 CPU/GPU。

本文从最容易独立验证的密度平流开始，实现一个标准 C++ 的半拉格朗日（semi-Lagrangian）示例，再把它放回不可压缩流体的一帧管线中，解释散度、压力泊松方程、投影、ping-pong 双缓冲及 GPU 性能边界。

1. 我们实际模拟的是什么？

实时烟雾通常不追踪每个分子，而采用欧拉网格（Eulerian grid）：空间被划成固定网格，每个位置保存场变量。

二维不可压缩流体常维护：

• 速度场 u = (u_x, u_y)；
• 压力场 p；
• 烟雾密度、温度等标量场；
• 障碍物或边界标记。

简化的不可压缩 Navier–Stokes 方程可以写成：

∂u/∂t + (u · ∇)u = -(1/ρ)∇p + ν∇²u + f
∇ · u = 0

其中平流项把速度随流场搬运，黏性项扩散动量，外力 f 加入浮力或用户扰动，压力项负责满足不可压缩约束。实时求解器通常通过操作分裂（operator splitting）把一大步拆成多个较简单步骤，而不是一次求完整方程。

标量密度若只被速度携带、不反过来影响速度，就是被动标量（passive scalar）。本文的最小程序只模拟这一部分，因此能展示“烟雾被搬运”，却还不是一个能自行产生流动的流体系统。

2. 一帧模拟为什么需要多道步骤？

一个常见的概念管线是：

速度 u
  │ 平流 / 黏性扩散 / 外力
  ↓
中间速度 u*（通常不再无散）
  │ 计算 divergence
  │ 解 pressure Poisson equation
  │ 减去 pressure gradient
  ↓
投影后速度 u_new（近似满足 ∇·u = 0）
  │
  └── 平流 density / temperature

不同方法可能调整操作顺序、离散方式或把部分步骤合并。图形学 Stable Fluids 路线强调稳定、规则和交互速度，通常接受较强数值耗散；工程 CFD 则可能更关注守恒、收敛阶和误差控制。不能把某套实时烟雾参数直接当成通用物理求解配置。

3. 平流为什么要从目标位置向后追踪？

若想计算新时刻网格点 x 上的密度，半拉格朗日方法沿速度反方向回到旧时刻位置：

x_previous = x - Δt · u(x)
d_new(x) = sample(d_old, x_previous)

它不是把源单元向前“撒”给邻居，而是让每个目标单元独立询问：形成我的那团流体上一时刻在哪里？这样避免多个源同时写一个目标，也很适合数据并行。

回溯位置通常落在格点之间，需要双线性或更高阶插值。线性插值稳定、实现简单，却会平滑尖锐特征；多次平流后烟雾变糊，往往不是渲染问题，而是数值耗散。

时间步可以无限增大吗？

半拉格朗日方法相对不容易因大时间步直接数值爆炸，但“稳定”不等于“准确”。一次跨越许多网格单元时，速度场沿路径的变化被粗略近似，插值误差和边界误差都会增加。时间步仍应根据速度尺度、网格间距和所需精度测试。

4. 如何写出一个最小可运行的平流版本？

下面使用 6×4 周期网格，网格间距视为 1。初始密度集中在 (2, 2)，恒定速度每秒向右移动一个网格；每步 dt = 0.5，执行两步。双线性插值会让密度形成 0.25, 0.50, 0.25 的分布。

代码使用 C++17，不依赖图形或 CUDA 库：

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>

int wrap(int coordinate, int extent) {
    const int remainder = coordinate % extent;
    return remainder < 0 ? remainder + extent : remainder;
}

float sample_periodic(
    const std::vector<float>& field,
    int width,
    int height,
    float x,
    float y) {
    const int x0 = static_cast<int>(std::floor(x));
    const int y0 = static_cast<int>(std::floor(y));
    const int x1 = x0 + 1;
    const int y1 = y0 + 1;
    const float tx = x - static_cast<float>(x0);
    const float ty = y - static_cast<float>(y0);

    const auto at = [&](int px, int py) {
        const int wrapped_x = wrap(px, width);
        const int wrapped_y = wrap(py, height);
        return field[static_cast<std::size_t>(wrapped_y * width + wrapped_x)];
    };

    const float bottom =
        (1.0F - tx) * at(x0, y0) + tx * at(x1, y0);
    const float top =
        (1.0F - tx) * at(x0, y1) + tx * at(x1, y1);
    return (1.0F - ty) * bottom + ty * top;
}

void advect_density(
    const std::vector<float>& source,
    std::vector<float>& destination,
    int width,
    int height,
    float velocity_x,
    float velocity_y,
    float dt) {
    for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            const float previous_x =
                static_cast<float>(x) - dt * velocity_x;
            const float previous_y =
                static_cast<float>(y) - dt * velocity_y;

            destination[static_cast<std::size_t>(y * width + x)] =
                sample_periodic(
                    source, width, height, previous_x, previous_y);
        }
    }
}

int main() {
    constexpr int width = 6;
    constexpr int height = 4;
    std::vector<float> source(
        static_cast<std::size_t>(width * height), 0.0F);
    std::vector<float> destination(source.size(), 0.0F);
    source[static_cast<std::size_t>(2 * width + 2)] = 1.0F;

    for (int step = 0; step < 2; ++step) {
        advect_density(
            source, destination,
            width, height,
            1.0F, 0.0F, 0.5F);
        source.swap(destination);
    }

    std::cout << std::fixed << std::setprecision(2);
    for (int x = 0; x < width; ++x) {
        std::cout << source[static_cast<std::size_t>(2 * width + x)]
                  << (x + 1 == width ? '\n' : ' ');
    }

    const double mass = std::accumulate(
        source.begin(), source.end(), 0.0);
    std::cout << "mass = " << mass << '\n';
}

编译运行：

clang++ -std=c++17 -O2 -Wall -Wextra -pedantic advection.cpp -o advection
./advection

预期输出：

0.00 0.00 0.25 0.50 0.25 0.00
mass = 1.00

这个特例中，恒定平移、周期边界和对称线性插值恰好保持总和。不能据此断言一般半拉格朗日平流质量守恒；空间变化速度、复杂边界和重复插值都可能改变总质量与峰值。

5. 这段代码为什么必须使用两份场？

每个目标格点都可能读取旧场中的四个邻居。若直接写回 source，前面更新的值会污染后面格点的采样：

// 错误思路：读写同一 field
field[index] = sample_periodic(field, ...);

结果会依赖遍历顺序，GPU 上还会形成并发读写竞态。ping-pong 双缓冲把角色分开：

step 0: source A --read--> compute --write--> destination B
        swap(A, B)
step 1: source B --read--> compute --write--> destination A

swap 只交换 vector 管理的缓冲区，不复制全部网格。GPU 实现同样通常交换设备指针、纹理对象角色或资源句柄。

主例的速度是常量。真实求解器会从速度场插值得到回溯速度，并可能使用 midpoint/RK2 等积分改善轨迹。速度分量放在格心还是面心，也会影响采样位置与离散一致性。

6. 边界条件为什么会改变整个问题？

示例使用周期边界：离开右侧会从左侧返回。这适合环形域测试，不适合有墙的烟雾盒。

常见边界包括：

边界	典型含义	需要处理的内容
周期边界	两侧首尾相连	索引 wrap
固壁、无穿透	流体不能穿过墙	法向速度与压力边界
无滑移壁面	壁面切向速度也受约束	黏性与速度 ghost cell
开放边界	流体可流出或流入	防反射与流入条件

简单地把越界坐标 clamp 到边缘，会隐含一种数值边界行为，不等于自动满足物理壁面。平流、散度、压力和投影必须使用相互一致的边界离散，否则压力解再多迭代也可能无法得到正确无散场。

障碍物还会把规则矩形域变成带掩码的区域。固体内部、流体—固体界面和狭窄通道需要单独定义 stencil，不能只在渲染时盖住密度。

7. 为什么平流之后还要压力投影？

外力、平流和数值误差会得到中间速度 u*，它通常不满足不可压缩条件 ∇·u = 0。投影法把它分解并减去压力梯度部分：

∇²p = (ρ / Δt) ∇·u*
u_new = u* - (Δt / ρ) ∇p

符号和尺度会随方程移项、网格间距与代码约定变化，但压力方程、梯度更新和散度离散必须成套一致。

典型离散流程为：

1. 用相邻速度计算每格散度；
2. 根据散度求压力泊松方程；
3. 用相邻压力计算梯度；
4. 从中间速度减去梯度；
5. 再次计算散度，验证是否降低。

最后一步不能省略。画面自然并不能证明投影正确，投影前后的散度范数才是直接诊断信号。

8. Jacobi 迭代次数为什么不能写成固定魔法数字？

二维五点 stencil 的 Jacobi 更新大致是：

p_new(i,j) =
    (p_left + p_right + p_bottom + p_top - h² · rhs(i,j)) / 4

每轮必须读旧压力、写新压力，所以同样需要 ping-pong。Jacobi 规则、易并行，却对低频误差收敛较慢；网格增大后，固定 40 次迭代未必达到相同精度。

更可靠的停止依据是压力方程残差：

r = rhs - ∇²p

可以监测 L2 或最大范数，并结合绝对/相对容差与最大迭代数。实时图形应用有时为了固定帧预算采用固定轮数，也应离线测量该轮数在不同分辨率和场景下的残差，而不是把历史教程参数当成普遍结论。

Gauss–Seidel、red-black、共轭梯度和 multigrid 各有并行性、收敛和实现复杂度取舍。规则 GPU 网格上的 multigrid 常能更有效处理多尺度误差，但不应在基础版本尚未验证时过早引入。

压力为什么可能差一个常数？

对某些封闭域 Neumann 边界，压力只确定到加法常数；给所有压力加同一个值不会改变梯度。求解器可能需要固定一个参考压力或移除均值，右端项也要满足兼容条件。只比较压力绝对值可能误判，应重点检查残差和投影后速度。

9. 网格位置为什么会影响压力与速度耦合？

把压力和两个速度分量都放在格心，称为 collocated grid。若离散不慎，可能产生棋盘格压力等奇偶解耦问题。

MAC staggered grid 常把：

pressure：单元中心
u_x：左右面中心
u_y：上下边中心

这样散度自然使用穿过单元面的通量，压力梯度也位于速度位置。代价是索引、边界和插值更复杂。教学原型可从格心布局开始理解流程，但工程实现应明确变量位置，不能只写“二维数组”而忽略采样坐标。

10. Stable Fluids 的“稳定”有什么边界？

半拉格朗日平流与隐式扩散使方法能容忍相对大的时间步，不容易像简单显式格式那样迅速爆炸。这种稳定性换来明显数值耗散：涡旋和密度细节会被插值抹平。

稳定不意味着：

• 满足质量或动量严格守恒；
• 任意大 dt 都准确；
• 任意边界处理都正确；
• 固定 Jacobi 轮数已经收敛；
• 低分辨率结果代表真实流体。

MacCormack/BFECC 等高阶修正可减少平流耗散，却可能产生过冲，需要 limiter；vorticity confinement 能在视觉上补回旋涡，但它添加的是模型力，不是恢复被丢失的真实解。每项增强都要用数值指标和视觉目标共同评估。

11. 如何系统验证每个步骤？

不要等整套烟雾渲染出来才调试。为每个算子构造可预测输入：

平流测试

• 零速度：新场应等于旧场；
• 常量速度 + 周期边界：密度应按已知方向平移；
• 常量场：平流后仍为常量；
• 监测总质量、最小值与最大值。

散度与投影测试

• 常量速度场的内部散度应接近零；
• 人工构造有散速度，投影后散度范数应明显下降；
• 压力残差应随迭代减少；
• 边界法向速度应满足约束。

收敛与稳健性测试

• 缩小 dt 或细化网格，结果应呈合理收敛趋势；
• 每帧检查 NaN/Inf；
• 记录质量、动能、最大速度和残差时间序列；
• CPU 参考实现与 GPU 实现比较容差。

“没有崩溃”不是数值测试。另一方面，实时视觉算法未必追求严格守恒，验收阈值应由实际目标定义并记录。

12. 为什么 GPU 很适合这些步骤，又常常受带宽限制？

平流、散度、Jacobi 和投影大多对每个格点执行相似 stencil，天然具有数据并行性。扁平连续数组或二维 CUDA 资源能让相邻线程访问相邻位置。

但每个格点通常只做少量算术，却读取多个邻居并写一个结果，算术强度不高。Jacobi 每多一轮都要再读写整个压力场，因此常受全局内存带宽限制。

优化方向包括：

• 使用 SoA 或合适向量类型匹配实际字段访问；
• 让线程布局与行连续方向一致，形成合并访存；
• 通过 texture object 获得二维寻址、缓存和硬件插值能力；
• 对产生复用的 stencil 使用 shared-memory tile；
• 减少不必要的中间场和 host/device 往返；
• 用更快收敛的压力求解器减少整场扫描。

纹理插值的坐标约定、归一化模式、寻址模式与精度依 CUDA 配置，必须用小网格测试验证。手写双线性插值更透明，但可能增加指令；两者应在相同边界语义下比较。

kernel 融合为什么不是越多越好？

融合能减少 launch 和中间读写，却可能增加寄存器、共享内存和控制流，降低 occupancy；有些阶段之间还存在全网格依赖，不能安全放进普通单 kernel 而省掉全局同步。先用 Nsight 找出带宽与启动瓶颈，再融合语义允许的相邻步骤。

13. CPU 版本怎样迁移到 CUDA？

正确路线通常是先让 CPU 小网格版本通过算子测试，再逐步迁移：

1. 用连续数组固定字段布局和边界函数；
2. 每个算子改为独立 kernel，保留相同输入输出；
3. 每迁移一步就与 CPU 结果比较；
4. 保持 ping-pong 指针角色清楚；
5. 最后才引入纹理、shared memory、异步 streams 或 kernel 融合。

CUDA kernel launch 是异步的。调试阶段要在启动后检查即时错误，并在合适同步点捕获执行错误。生产管线可使用同一 stream 的有序依赖、事件和少量必要同步，避免每个 kernel 后都全设备同步。

3D 网格把内存和计算量从 N² 提高到 N³，压力求解与体渲染成本迅速增大。不要在二维离散、边界和验证尚未可靠时直接扩展三维。

14. 常见误区如何纠正？

误区一：烟雾形状自然，方程就正确

不对。渲染、插值和外力可以掩盖散度或边界错误。必须检查残差、散度、质量与已知解。

误区二：半拉格朗日无条件稳定，所以 dt 随便选

不对。它可能不爆炸，却会在大步长下严重失真。稳定性与准确性是不同要求。

误区三：Jacobi 固定迭代 40 次就够

不对。所需轮数随网格、边界、初值和容差变化。固定预算也要测残差和投影后散度。

误区四：平流可以原地更新

通常不可以。旧场要被所有目标点一致读取，应使用双缓冲或证明安全的专用算法。

误区五：投影后散度一定精确为零

不对。离散误差、边界、压力残差和浮点误差都会留下非零散度；目标是满足明确容差。

误区六：GPU 版本只需把双重循环改成 kernel

不对。还要管理全局依赖、资源生命周期、带宽、同步、错误传播和 CPU/GPU 结果差异。

15. 什么时候适合使用这套方法？

规则网格 Stable Fluids 路线适合：

• 游戏、交互艺术和可视化中的烟雾、染料、火焰基础效果；
• 规则二维/三维域，追求稳定与实时吞吐；
• 能接受一定数值耗散和近似压力解；
• 希望利用 GPU stencil 并行。

它不应未经评估直接用于：

• 需要严格质量/动量守恒的工程分析；
• 高雷诺数湍流的可信预测；
• 复杂自由表面、强多相耦合或精细边界；
• 需要可验证误差界和网格收敛报告的安全关键模拟。

这些问题可能需要有限体积、有限元、粒子方法、自由表面跟踪和经过验证的专业求解器。

16. 总结

回到开头，烟雾会移动只证明平流路径产生了某种结果。完整不可压缩流体还要让速度经过压力投影满足离散无散约束，并让所有算子使用一致网格、边界和时间尺度。

最重要的结论是：

1. 半拉格朗日通过反向追踪和插值完成平流，稳定但会产生数值耗散。
2. stencil 更新通常不能原地进行，ping-pong 双缓冲保证所有格点读取同一旧状态。
3. 压力泊松解的质量应由残差和投影后散度验证，迭代次数不是魔法常量。
4. 边界条件与变量在网格上的位置是离散方程的一部分，不是渲染细节。
5. GPU 能并行处理规则网格，却常受带宽与全场迭代限制，优化必须建立在正确参考实现上。

实践中先做本文这样的 4×6 小测试：零速度、常量平移、周期 wrap 都得到可手算结果后，再添加散度与投影。每增加一个算子，就增加一个数值不变量；这样比对着烟雾画面猜 bug 快得多。

参考资料

• Jos Stam, Stable Fluids
• NVIDIA GPU Gems：Fast Fluid Dynamics Simulation on the GPU