先把仿真跑起来:激波与烟羽
1. 这一章讲什么
三件事: 一次仿真的「一步」不是一个公式,是一条五件事的流水线; 这五件事分别在算什么、各自的输入输出长什么样; 以及哪一步最贵——这一点单独拎出来,因为后面所有「让网络去救求解器」的做法, 下手的地方就是那一步。
它在全书链条里的位置:这是第 02 章的兑现。 第 02 章说「切成格子、一步步往前推」, 但没说一步里发生了什么。这一章把两个真实的仿真从头跑到尾,每一步都带真数字。
为什么值得花一整章在「怎么正着跑一遍仿真」上: 后面说「把求解器搬进训练过程」,搬的就是这条流水线; 说「让求解器交出梯度」,是给这五步各配一个求导的办法。 没有这一章,后面十几章的「求解器」只是一个词。
2. 顶层全景
这一章跑两个仿真,由简到繁。
① 一维 Burgers(热身) ② 二维浮力烟羽(全书的主力场景)
128 个格子上的一串数 32×40 个格子上的两个场
│ │
┌────┴────┐ ┌────────┴────────┐
│ 扩散一次│ │ 1. 平流烟雾 │ 烟被流带走
│ (抹平) │ │ 2. 加浮力 │ 热的地方往上顶
├─────────┤ │ 3. 平流速度 │ 速度被自己带走
│ 平流一次│ │ 4. 扩散 │ 黏性把速度抹平
│ (搬运) │ │ 5. 压力投影 ★ │ 强制每格净流出为零
└────┬────┘ └────────┬────────┘
│ 重复 32 次 │ 重复 N 次
▼ ▼
中心鼓出一个陡坎 烟柱上升、撞顶壁、往右弯
(激波) 最大速度 0.463 → 5.484
这张图在讲什么: 一步仿真就是几个小算子按固定顺序各作用一次。 把复杂的方程拆成几个各管一件事的算子轮流上,这种做法叫算子分裂(operator splitting)。 打星号的那一步是整条流水线里最贵的,它是本章的落点。
下面出现的所有数字都是书里那两个例子真跑出来的输出,来源逐条标在脚注里。
3. 热身:一维的一根线,跑 32 步
先看设置,它简单到可以整个记住。
一根线段,从 −1 到 +1,切成 128 个格子;每个格子里存一个数,表示那一处的速度。 时间上跑 1 个时间单位,分成 32 步,所以每步 Δt = 1/321。
边界用的是周期边界(periodic boundary): 线段的右端和左端接在一起,从右边出去的东西从左边回来。 这是最省事的一种边界条件——不用操心边缘上该发生什么。
初始状态是一条正弦波:u(x) = −sin(πx)。 摊开看是这样:
线段左半边的速度是正的(往右冲),右半边是负的(往左冲)。
两股东西相向而行,而中间没有任何机制拦着——第 02 章说过,这种情况会撞出激波。
书打印了初始状态第 10 到 14 个格子的值,可以直接抄下来当起点:
0.4929 / 0.5350 / 0.5758 / 0.6152 / 0.6533——一串正数,在往右冲2。
每一步只做两件事,顺序固定:
- 扩散一次: 把速度抹平一点点。用的算法叫中心差分——
一格的变化由它左右两个邻居共同决定(左邻加右邻减去两倍自己)。
「显式」的意思是只用当前这一步已知的值算出下一步,不解方程组,便宜;
黏性系数 ν 取
0.01/(128π),是个很小的数,所以抹得很轻1; - 平流一次: 把速度顺着流速搬走。这一步下一节单独讲。
跑完 32 步之后,书打印了前五个格子的值:
0.0057 / 0.0172 / 0.0286 / 0.0401 / 0.05153。
把这两串数对比着读,一整个演化就出来了: 起点那一带原本是 0.49 到 0.65 的大正值,一秒之后前五格只剩 0.006 到 0.05—— 边缘的速度几乎被抽干了,因为那些东西全都被搬到中间去了。
而中间发生了什么:左边的正速度往右冲、中心右边的负速度往左冲,撞在一起, 在中心堆出一个几乎垂直的陡坎4。这就是第 02 章那个激波,现在它是真跑出来的。
另起一处的走查(一维 Burgers): 128 格 × 32 步 → 每步「中心差分扩散 + 半拉格朗日平流」→ 第 10–14 格从
0.4929…0.6533掉到前五格的0.0057…0.0515,中心鼓出一个陡坎。 这一处的每个数都是书里的真实打印输出。
4. 半拉格朗日平流:反过来问「这一格的东西是从哪儿飘来的」
先看一个会出事的做法。 要把某一格的东西往下游搬,最自然的想法是: 算出这一格的东西这一步会飘多远,然后把它加到目标格子里去。
这个做法的麻烦在于:飘的距离通常不是格子宽度的整数倍。 东西会落在两格中间,你得把它劈开分给两个邻居; 劈得不好,总量就会一点点多出来或者少掉,几步之后整个仿真就飞了。
半拉格朗日平流把这个问题反过来问: 不问「我这一格的东西要去哪」,而是问 「我这一格现在应该有的值,上一步是从哪个位置飘过来的」。
做法只有两步: ① 从这一格的中心出发,顺着速度往回退 Δt 的距离, 落到某个位置上(通常不在格子中心);② 在旧的那一层数据上采样那个位置的值 (落在两格之间就按距离加权取两边的平均),把采到的值直接写进这一格。
为什么这样就稳:因为每一格都只是「去旧数据里取一个值」,不涉及往外发东西。 不管速度多大,取到的值总是旧数据里已有的值的加权平均—— 它永远不会造出比旧数据更大的值来,所以炸不了。
书对这个算子的定语正是「一阶的、稳定的近似」5。 「一阶」是它的代价:精度只有一阶,东西被搬得有点糊—— 这个糊后面会变成一个真问题,第 10 章那个「让小网络去救粗糙求解器」救的就是它。
5. 二维烟羽:一步里的五件事
换到全书的主力场景:一团热烟从底部偏左的位置冒出来,往上升。
设置: 横向 32 格、纵向 40 格,覆盖的物理区域横 0 到 80、纵 0 到 100; Δt = 1.5,黏性 ν = 0.01。 入流是一个半径 10、圆心在 (30, 15) 的圆——注意 30 偏在左边(横向中点是 40), 这个偏心一会儿会有后果6。
要跟踪两个场: 速度场,和一个表示「这儿有多热」的标记场 (书用它代替真的算密度,就是第 02 章那个布森涅斯克近似)。
一步 step() 里,五个算子按这个顺序各作用一次7:
| # | 做什么 | 用大白话说 |
|---|---|---|
| 1 | 半拉格朗日平流,作用在标记场上 | 烟被当前的流带走,然后在入流处补一圈新烟 |
| 2 | 按标记场的大小加一份竖直向上的力 | 热的地方往上顶 |
| 3 | 半拉格朗日平流,作用在速度场上 | 速度被它自己带走(这就是那个非线性项) |
| 4 | 显式扩散 | 黏性把速度场抹平一点 |
| 5 | 压力投影 | 强制「每格净流出为零」 |
第 3 步值得停一下:速度被自己搬运。 这一句就是第 02 章说的「非线性」的具体样子—— 搬运的工具和被搬运的东西是同一个场,所以你没法把它拆成「输入乘一个系数」。
6. 交错网格:速度不存在格子中心
先看一个会出事的做法。 最自然的存法是:每个格子的中心存一个烟雾浓度、 再存一个速度(横向分量和纵向分量各一个数)。
麻烦出在 算散度的时候。 散度是「这一格的净流出」, 它问的是穿过格子边界的流量——而你手上只有格子中心的速度。 只能拿两个中心的速度去估边界上的值,估出来的散度不准, 于是「强制散度为零」这一步就变成了强制一个不准的量为零。
交错网格(staggered grid)把速度挪了个位置: 烟雾浓度存在格心,速度的横向分量存在格子的左右两条竖边上、 纵向分量存在上下两条横边上——正好就是流量真正穿过的地方8。
代价是各个量的个数对不上,而书把这个数原样打印了出来: 标记场是 32 × 40;而速度的横向分量是 31 × 40、纵向分量是 32 × 399。 少的那一格是因为 32 个格子之间只有 31 条内部竖边。
这个「对不上」不是瑕疵,它是交错网格的标志。 你以后在任何流体代码里看到「x 分量比 y 分量少一行」,就是这件事。
7. 压力投影:整条流水线上最贵的一步
结论先行:这一步就是把第 02 章那条硬约束真的强制执行一次, 而它是全书所有「混合方法」下手的地方。
先看现象。 走完前四步之后,速度场不满足「每格净流出为零」—— 平流会把东西挤到一起,浮力会凭空往上顶,这两件事都不管守恒。 所以每一步的最后必须补一刀:把当前这个速度场,改成一个最接近它、 但每一格净流出确实为零的速度场。
「投影」这个词就是从这儿来的: 你在一大堆可能的速度场里, 挑出那个满足约束的、离你手上这个最近的。 就像把空间里一个点垂直落到一个平面上,落点是那个平面上离它最近的点。
怎么落下去?靠一个中间量:压力。 具体做法是 先解出一个压力场,再用这个压力场的梯度去改速度10。 压力在这里不是「物理上真实的压强」那么直观的东西,它更像一个中间变量: 它的作用就是提供一份处处合适的推力,把多余的净流出推平。
贵就贵在「解出那个压力场」这一步,它要解一个泊松方程。 泊松方程是这样一种方程:某一格的未知量,由它周围邻居的未知量共同决定—— 每一格都这么写,于是所有格子上的未知量扭成了一整个方程组。
关键在这个「一整个」。 前面四步都是局部的: 一格的新值只用到它自己和几个邻居的旧值,算完就完。 而泊松方程是全局的:随便改动任何一格,原则上会牵动全场每一格。 这类「一处变动瞬间牵动全场」的方程叫椭圆型方程—— 这个分类到第 04 章选网络架构时会变成一条硬判据。
书的原话很干脆:压力投影这一步「通常是上面 那个序列里计算上最贵的一步」10。
边界上怎么办?这个烟羽用的是封闭容器: 速度在四壁上被摁成零 (这类「直接规定边界上的值」的条件叫狄利克雷边界), 压力则规定它在边界方向上的变化率为零11。
主走查第 1 步: 一个 32×40 的状态(格心存烟雾浓度、竖边存横向速度、 横边存纵向速度)进入
step()→ 五个算子依次作用 → 第 5 步解一个牵动全部 1280 格的方程 → 吐出 Δt = 1.5 之后的新状态。
8. 十帧真数:烟羽是怎么撞到顶壁往右弯的
书连着调了十次 step(),把每一帧的最大速度打印了出来12:
| 帧 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|---|
| 最大速度 | 0.463 | 0.897 | 1.410 | 2.041 | 2.928 | 3.839 | 4.527 | 4.868 | 5.131 | 5.484 |
这串数怎么读:头几帧几乎在翻倍(0.463 → 0.897 → 1.410 → 2.041), 后几帧涨幅明显收住了(4.527 → 4.868 → 5.131 → 5.484,每帧只涨三四个百分点)。
前半段是浮力在持续加速: 热的地方一直被往上顶,速度一路累积。 后半段涨不动了,因为烟柱已经顶到了容器的上壁—— 上壁的速度被摁成零(第 7 节那个狄利克雷边界),往上冲的势头无处可去。
然后是这一节最有意思的一句:烟羽撞到顶壁之后会往右弯。 原因是入流的圆心在 x = 30,而横向区域是 0 到 80——它偏左13。 往上冲的柱子一撞顶就得往两边摊开,而左边空间小、右边空间大,于是整体朝右倾。
这句话的分量在于:那个弯是从五个算子里长出来的,没有任何一行代码写着「往右弯」。 你能在方程里指出「往右弯」这个行为在哪一项吗?指不出来。 它是这五个算子反复作用 N 次之后涌现的结果。
主走查第 2 步: 上面那个状态被反复送进
step()十次 → 最大速度 0.463 → 5.484(涨了约 12 倍),头几帧翻倍、后几帧收住 → 烟柱顶到上壁 → 因为入流偏左(x = 30 落在 0–80 里),整体朝右弯。
9. 这一章在全书里到底是干什么用的
别把它当成一段工具教程,它是后面十几章的名词表。
| 后面会说 | 说的其实是这一章的什么 |
|---|---|
| 「把求解器搬进训练环路」 | 把第 5 节那条五步流水线搬进训练环路 |
| 「让求解器交出梯度」 | 给那五个算子各配一个「输入变一点、输出变多少」的算法 |
| 「用一个粗糙的低保真求解器 + 网络修正」 | 把格子数调小、把第 4 节那个一阶平流的糊留着,让网络去补 |
| 「硬约束」 | 第 7 节那个投影动作 |
| 「反问题」 | 不给 t = 0 的状态,只给第 20 帧的形状,反推初始速度 |
10. 作者的判断与证据
书里给了证据的(全是真实运行输出):
| 说法 | 证据 |
|---|---|
| 中心会形成激波 | 32 步之后画出来的图,以及边缘速度被抽干的那串数34 |
| 烟羽最大速度的十帧演化 | 逐帧打印:0.463 → 5.48412 |
| 交错 网格上各分量个数不同 | 打印出来的形状:32×40 / 31×40 / 32×399 |
| 压力投影是最贵的一步 | 书直接陈述,没有给计时数据10 |
作者的判断、书没给数据的:
- 「压力投影通常是最贵的一步」——这一句是本章的落点,而书没有给任何计时。 它是数值流体领域的常识性说法,但在这本书里它是一句断言。
判断(我们的,不是书里的): 这句断言我们认为成立,理由不在书里而在结构上—— 前四步是局部的、每格只看几个邻居,一次扫过全部格子就完事; 而第五步要解一个把所有格子拴在一起的方程组,通常得靠迭代法反复扫很多遍。 一个是扫一遍,一个是扫很多遍,量级差在这里。 如果错,会错在: 如果格子数很少(比如这里只有 1280 格), 迭代法几步就收敛,压力投影可能并不显著地贵; 这句话真正成立的是大规模三维仿真,那才是它被反复强调的场景。 判据是:格子数越多,这一步占的比例越高。
11. 边界与局限
- 这两个例子都极小。 128 个格子和 1280 个格子,在真实的工程仿真里是玩具规模 (三维湍流动辄上千万个单元)。它们的作用是把机制亮清楚,不是给性能参考;
- 书没有讲这些算子内部怎么实现。 半拉格朗日平流的采样细节、 泊松方程用哪种迭代法解,书都交给了它用的那个仿真框架;
- 第 4 节那个「一阶平流会糊」这一句是我们补的——书只说它是「一阶的、稳定的」, 没有展开说一阶意味着什么代价;
- 激波在这里只是「看见了」,没有被量化。 书没有给激波位置或陡度的数;
- 这一章完全没有神经网络。 它是纯数值仿真,是后面一切的地基。
12. 可带走的
- 一步仿真 = 几个算子按固定顺序各作用一次。 把复杂方程拆成几个各管一件事的算子, 这个做法叫算子分裂;
- 烟羽那条流水线的五步要背下来: 平流标记场 → 加浮力 → 平流速度 → 扩散 → 压力投影。 后面十几章说的「求解器」就是这五步;
- 半拉格朗日平流反过来问「这一格的东西上一步从哪儿来」,所以它稳。 代价是一阶精度,东西会被搬糊;
- 交错网格:浓度在格心,速度在格子的边上。 各分量个数对不上是它的标志,不是 bug;
- 压力投影 = 把速度场投影到「每格净流出为零」的那一套里去。 做法是解出一个压力场,再用它的梯度改速度;
- 它是整条流水线里最贵的一步,因为它要解一个牵动全场的方程组—— 前四步都只看邻居,只有它是全局的;
- 「一处变动瞬间牵动全场」的方程叫椭圆型,这个分类在第 04 章会变成选网络的判据;
- 复杂行为不写在任何一行代码里。 烟羽撞顶往右弯,是五个算子反复作用的涌现结果—— 而它的起因只是入流的圆心偏左了 10 个单位。
13. 原文地图
| 主题 | 原书章 | 原文位置 |
|---|---|---|
| Burgers 的格子数、步数、黏性、初始状态 | 3.6.2 Importing and loading phiflow | text/02-p21-40.txt:658(搜「N=128 cells」) · text/02-p21-40.txt:671(搜「NU = 0.01」) |
| 初始状态第 10–14 格的值 | 3.6.2 同上 | text/02-p21-40.txt:705(搜「Velocity tensor entries 10 to 14」) |
| 两个算子:显式扩散 + 半拉格朗日平流 | 3.6.3 Running the simulation | text/02-p21-40.txt:718(搜「stable first-order approximation」) |
| t = 1.0 时前五格的值 | 3.6.3 同上 | text/02-p21-40.txt:732(搜「New velocity content at t=1.0」) |
| 中心形成激波、两个速度包相撞 | 3.6.4 Visualization | text/03-p41-60.txt:28(搜「shock developing in the center」) |
| 烟羽的方程、浮力项与边界条件 | 3.7 Navier-Stokes Forward Simulation | text/03-p41-60.txt:88(搜「simple buoyancy model」) · text/03-p41-60.txt:92(搜「Dirichlet boundary conditions」) |
| 格子数、Δt、入流位置 | 3.7.2 Setting up the simulation | text/03-p41-60.txt:110(搜「40 × 32 cells」) · text/03-p41-60.txt:119(搜「radius=10」) |
| 一步里的五个算子 | 3.7.2 同上 | text/03-p41-60.txt:147(搜「advect.semi_lagrangian(smoke」) · text/03-p41-60.txt:166(搜「advected the smoke field」) |
| 压力投影最贵、解泊松方程 | 3.7.2 同上 | text/03-p41-60.txt:171(搜「computationally most expensive step」) |
| 交错网格上各分量个数不同 | 3.7.3 Datatypes and dimensions | text/03-p41-60.txt:226(搜「non-uniform shape」) |
| 十帧的最大速度 | 3.7.4 Time evolution | text/03-p41-60.txt:238(搜「Computed frame 0, max velocity」) |
| 入流偏左、撞顶往右弯 | 3.7.4 同上 | text/03-p41-60.txt:258(搜「curve towards the right」) |