ρ 和 u 连末位都对得上，只有应变率虚胖了 45% —

能原封不动跨过层级边界的只有 ρ 和 u。非平衡部分必须改写为原来的 τ_f/(m·τ_c) 倍。

ρ 和 u 连末位都一致，唯独应力鼓了起来#

在格子玻尔兹曼(LBM)代码上接入自适应网格加密(AMR)，就会多出"层级边界"这个新东西。边界处有一层粗网格与细网格互相重叠。这一层里的值必须从一侧誊写到另一侧。

验证誊写是否正确时，人们通常只看密度和速度。可只看这两个，bug 会顺利蒙混过关。把分布函数 $f_i$ 原样复制到细层， $\rho$ 和 $\mathbf{u}$ 依然能对到小数末位。对不上的是应力。本文要讲清这个差距究竟是几倍，以及为什么偏偏是那个倍数。

先把结论摆出来：倍数是 $\tilde\tau_f/(m\tilde\tau_c)$ 的倒数。常见设置下约 1.45 倍，高雷诺数设置下接近 2 倍。

跨越层级必须相同的量，与必须改变的量

层级之间做转换的原则只有一条：守恒量必须保持相同。 $\rho$ 、 $\mathbf{u}$ 、 $p$ 都是物理量，不会因为换了网格就有理由改变。

把加密比记作 $m = \Delta x_c/\Delta x_f$ 。下标 $c$ 指粗网格， $f$ 指细网格。采用声学标度(acoustic scaling)时， $\Delta t$ 也按同一比例缩小。

\Delta x_f = \frac{\Delta x_c}{m}, \qquad \Delta t_f = \frac{\Delta t_c}{m}

这样一来，格子单位速度 $\hat{u} = u\,\Delta t/\Delta x$ 在两个层级上相同。 $\rho$ 也相同，于是平衡分布 $f_i^{\mathrm{eq}}(\rho, \hat{u})$ 在两个层级上是完全同一个值。需要改动的只有偏离平衡的那部分，也就是 $f_i^{\mathrm{neq}} = f_i - f_i^{\mathrm{eq}}$ 。

下面的模拟可以直接动手调。

tau_c0.80m2

Both lattices start from the same sine and advance over the same physical time. With tau rescaled the orange curve sits on the blue one; press tau copied and the fine level suddenly relaxes a different fluid. At step 0 the amplitudes are 1.0000 and 1.0000; the measured flux ratio is 0.0000 against the predicted r = 0.6875. Drag tau_c towards 0.51 and r falls to 1/2 — that is where copying a population costs you a factor of two.

同一块物理区域、同一段物理时间，由两套真实网格并排求解。拖动 tau_c 滑块，观察右侧梯子上的 measured q_f/q_c 是否一直贴着 rescale r。按下 tau copied 按钮，细层就变成了另一种黏度的流体，橙色曲线随即从蓝色曲线上脱开。

一旦把黏度钉住，τ 就得跟着动

物理黏度不能因为换了层级就改变。用格子单位弛豫时间 $\tilde\tau$ 写出的黏度是这样的。

\nu = c_s^2\left(\tilde\tau - \tfrac{1}{2}\right)\frac{\Delta x^2}{\Delta t}

$c_s^2$ 是格子声速的平方， $\Delta x^2/\Delta t$ 是扩散系数的量纲。那个 $-\tfrac{1}{2}$ 是梯形积分留下的痕迹，在LBM 离散化留下的 Δt/2里单独讲过。

在声学标度下 $\Delta x^2/\Delta t$ 变为原来的 $1/m$ 倍。加上 $\nu_f = \nu_c$ 的约束， $\tilde\tau$ 就被牵了出来。

\tilde\tau_f = \frac{1}{2} + m\left(\tilde\tau_c - \frac{1}{2}\right)

取 $m=2$ 、 $\tilde\tau_c = 0.8$ ，得到 $\tilde\tau_f = 1.1$ 。关键在于它不是两倍。若把 $\tilde\tau$ 直接乘以 $m$ ， $\nu$ 就对不上了。因为 $\tfrac{1}{2}$ 不参与标度。

非平衡部分上挂着的不只是 τ，还有 Δt#

把 Chapman–Enskog 展开的一阶项写成 Grad 近似，非平衡部分正比于应变率。

f_i^{\mathrm{neq}} = -\frac{w_i \rho \tilde\tau}{c_s^2}\left(\mathbf{c}_i\mathbf{c}_i - c_s^2\mathbf{I}\right) : \hat{\mathbf{S}}

$w_i$ 是权重， $\mathbf{c}_i$ 是格子速度， $\hat{\mathbf{S}}$ 是格子单位应变率。它与物理应变率 $\mathbf{S}$ 的关系是 $\hat{\mathbf{S}} = \mathbf{S}\,\Delta t$ 。同一点处的物理应变率与层级无关，因此 $f^{\mathrm{neq}}$ 正比于 $\tilde\tau\,\Delta t$ 。

r \equiv \frac{f_i^{\mathrm{neq},f}}{f_i^{\mathrm{neq},c}} = \frac{\tilde\tau_f\,\Delta t_f}{\tilde\tau_c\,\Delta t_c} = \frac{\tilde\tau_f}{m\,\tilde\tau_c}

$r$ 就是重标定系数，与 Dupuis 和 Chopard 整理出的形式一致。取 $\tilde\tau_c = 0.8$ 、 $m=2$ ，则 $r = 1.1/1.6 = 0.6875$ 。细层的非平衡量只有粗层的 69%。

把三种传递方式摆进同一张表

在层级边界上把值交出去的做法，工程实践里分成三种。

方式	传递的内容	需要的信息	代价	失效之处
整体复制	$f_i$ 原样	无	最低	$\rho,\mathbf{u}$ 正确，应力虚胖 $1/r$ 倍
宏观量 + Grad 重构	插值 $\rho, \mathbf{u}$ 后重算 $f^{\mathrm{neq}}$	速度梯度	中等	梯度得用有限差分重新求一遍
非平衡插值 + 缩放	$f^{\mathrm{eq}}$ 重算， $f^{\mathrm{neq}}$ 乘 $r$	$\tilde\tau_c, \tilde\tau_f, m$	接近最低	$r$ 里的 $m$ 很容易漏掉

第二种和第三种的结果应当一致。第三种之所以成为工程标准，是因为不必重新求梯度。只要有两个 $\tilde\tau$ 和一个 $m$ ，系数就出来了。

空间插值本身很普通。粗→细在三维里用三线性(trilinear)插值，细→粗取 $2^d$ 个单元的平均。 $d$ 是维数。难的不是插值，而是插完之后到底该给谁乘上 $r$ 。

用 Python 把复原出的应变率量一遍#

在 D2Q9 格子上先定死一个物理应变率，再分别造出两个层级的 $f^{\mathrm{neq}}$ 。然后用细层的 $\tilde\tau_f$ 把应变率复原出来，把重标定过的值和直接复制的值并排放进去比。

CS2 = 1.0 / 3.0
EX = [0, 1, 0, -1, 0, 1, -1, -1, 1]
EY = [0, 0, 1, 0, -1, 1, 1, -1, -1]
W = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
 
 
def equilibrium(rho, ux, uy):
    """D2Q9 平衡分布。层级不同，只要 rho、u 相同就得到相同的值。"""
    out = []
    u2 = ux*ux + uy*uy
    for i in range(9):
        cu = EX[i]*ux + EY[i]*uy
        out.append(rho*W[i]*(1 + cu/CS2 + cu*cu/(2*CS2*CS2) - u2/(2*CS2)))
    return out
 
 
def grad_neq(rho, tau, sxy):
    """用 Grad 近似造出的非平衡部分。sxy 是该层级的格子单位剪切应变率。"""
    out = []
    for i in range(9):
        q_xy = EX[i]*EY[i]            # Q_i 的 xy 分量（对角项与 sxy 配不上对）
        out.append(-(W[i]*rho*tau/CS2) * 2.0 * q_xy * sxy)
    return out
 
 
def recover_strain(fneq, rho, tau, dt):
    """从非平衡矩里复原出物理单位的剪切应变率。"""
    pi_xy = sum(fneq[i]*EX[i]*EY[i] for i in range(9))
    return -pi_xy / (2.0*rho*CS2*tau*dt)
 
 
def tau_on_level(tau_c, m):
    """黏度钉住时，细层应当具有的弛豫时间。"""
    return 0.5 + m*(tau_c - 0.5)
 
 
def nu_physical(tau, dx, dt):
    return CS2*(tau - 0.5)*dx*dx/dt
 
 
rho, ux, uy = 1.0, 0.05, 0.0
s_phys = 0.004          # 物理剪切应变率 —— 与层级无关的量
m = 2                   # 加密比
dx_c, dt_c = 1.0, 1.0
dx_f, dt_f = dx_c/m, dt_c/m
 
feq_c = equilibrium(rho, ux, uy)
feq_f = equilibrium(rho, ux, uy)
print("max |feq_c - feq_f| = %.3e" % max(abs(a-b) for a, b in zip(feq_c, feq_f)))
print()
print("%5s  %6s  %8s  %8s  %7s  %11s  %14s  %6s" % (
    "tau_c", "tau_f", "nu_c", "nu_f", "r_meas", "tf/(m*tc)", "copied/true", "err%"))
for tau_c in [0.51, 0.55, 0.60, 0.80, 1.20, 2.00]:
    tau_f = tau_on_level(tau_c, m)
    fneq_c = grad_neq(rho, tau_c, s_phys*dt_c)
    fneq_f = grad_neq(rho, tau_f, s_phys*dt_f)
    r = fneq_f[5]/fneq_c[5]
    s_ok = recover_strain(fneq_f, rho, tau_f, dt_f)
    s_bad = recover_strain(fneq_c, rho, tau_f, dt_f)
    print("%5.2f  %6.3f  %8.5f  %8.5f  %7.4f  %11.4f  %14.4f  %6.1f" % (
        tau_c, tau_f, nu_physical(tau_c, dx_c, dt_c), nu_physical(tau_f, dx_f, dt_f),
        r, tau_f/(m*tau_c), s_bad/s_ok, (s_bad/s_ok - 1)*100))
 
print()
worst = grad_neq(rho, 0.51, s_phys)
print("mass moment of f^neq     = %.3e" % sum(worst))
print("momentum moments of f^neq = %.3e, %.3e" % (sum(worst[i]*EX[i] for i in range(9)),
                                            sum(worst[i]*EY[i] for i in range(9))))

max |feq_c - feq_f| = 0.000e+00
 
tau_c   tau_f      nu_c      nu_f   r_meas    tf/(m*tc)     copied/true    err%
 0.51   0.520   0.00333   0.00333   0.5098       0.5098          1.9615    96.2
 0.55   0.600   0.01667   0.01667   0.5455       0.5455          1.8333    83.3
 0.60   0.700   0.03333   0.03333   0.5833       0.5833          1.7143    71.4
 0.80   1.100   0.10000   0.10000   0.6875       0.6875          1.4545    45.5
 1.20   1.900   0.23333   0.23333   0.7917       0.7917          1.2632    26.3
 2.00   3.500   0.50000   0.50000   0.8750       0.8750          1.1429    14.3
 
mass moment of f^neq     = 0.000e+00
momentum moments of f^neq = 0.000e+00, 0.000e+00

三件事一次性摆在了眼前。平衡分布在两个层级上完全相同。 $\nu_c$ 与 $\nu_f$ 两列连位数都吻合。测得的 $r$ 与闭式 $\tilde\tau_f/(m\tilde\tau_c)$ 一致。

最后两行正是本文标题的出处。 $f^{\mathrm{neq}}$ 的零阶矩和一阶矩精确为零。无论是复制还是重标定， $\rho$ 和动量都纹丝不动。出错的只有二阶矩这一处。

τ 越贴近 0.5，复制的代价越是翻倍#

从上往下读误差那一列，方向就显出来了。 $\tilde\tau_c = 2.0$ 时是 14%。 $\tilde\tau_c = 0.51$ 时是 96%。写成闭式， $r$ 的取值范围就清楚了。

r = \frac{1/2 + m(\tilde\tau_c - 1/2)}{m\,\tilde\tau_c}

$\tilde\tau_c \to 1/2$ 时 $r \to 1/2$ ， $\tilde\tau_c \to \infty$ 时 $r \to 1$ 。黏度大的时候，复制了也看不出破绽。黏度小的时候，应力直接翻倍。

麻烦在于，人们上 AMR 的理由大多正是高雷诺数。 $\tilde\tau$ 会被压到 0.5 附近来用。复制这个 bug 炸得最狠的区域，恰好就是最想用 AMR 的区域。症状也容易误导。质量和动量守恒，所以并不发散，只在层级分界线上留下一层薄薄的涡量。

把 $m$ 调大会更糟。 $m=4$ 、 $\tilde\tau_c = 0.8$ 时 $\tilde\tau_f = 1.7$ 。于是 $r = 1.7/3.2 = 0.531$ ，复制误差跳到 88%。

多加一个层级，账单按 m^(d+1) 开出来#

重标定只发生两次：从粗层下沉到细层的 explosion，以及从细层上浮回粗层的 coalescence。这两者之间就是普通的 collide-and-stream。

下面这只时钟可以一步一步走完一个周期。

m2dwell42

Use next phase to walk the cycle one gate at a time: explosion, 2 fine sub-steps, coalescence. The two dashed lines are the only moments populations cross levels — everything between them is ordinary collide-and-stream. Raise m or switch to d = 3 and the work factor climbs as m^(d+1) = 8; currently at cycle 0, phase explosion.

按 next phase，确认 explosion → $m$ 次子步 → coalescence 的先后顺序。绿色和紫色虚线，是分布函数跨越层级的仅有两个时刻。把 m 和 d 调大，右侧账本里的 work factor 就按 $m^{d+1}$ 往上走。

这个指数才是 AMR 设计中真正的约束。三维下取 $m=2$ ，一块补丁就贵了 16 倍。所以细层放在哪里、放多少，对性能的影响远大于格式的选择。

那么重叠层为什么只能是一层

重叠层是两个层级的节点同时存在的那一层单元。为什么只要一层？因为流动步一次只挪一格。细层每走一个子步，从外部涌进来的分布函数正好来自一格之外。把那一格填上就够了。

这个视角和边界节点丢失的分布函数是同一套。无论是壁面还是层级边界，第一件事都是数清"流动之后哪些 $f_i$ 是空的"。在壁面上由几何形状给出答案，在层级边界上则由加密比 $m$ 决定。

改用单元中心(cell-centered)网格，重叠层比节点中心(node-centered)更好处理。没有重合的节点，归属权清晰，并行剖分时的通信对象就落成一份单元清单。代价是粗→细插值中位置会错开半格，插值模板必须照此重新摆放。

无论是新写一套 AMR 界面还是读别人的代码，三行就能诊断完。 $\tilde\tau_f$ 是不是 $\tfrac{1}{2} + m(\tilde\tau_c - \tfrac{1}{2})$ ？ $f^{\mathrm{neq}}$ 有没有乘上 $r$ ？ $r$ 的分母里有没有 $m$ ？漏得最多的是第三条。

ρ 和 u 连末位都对得上，只有应变率虚胖了 45% —— LBM 网格加密的非平衡重标定