用 18% 的网格算出同样的答案 —— AMR 打标判据与 coarse-fine 重通量

缓冲带窄于 CFL 乘 regrid 次数时，界面会走出 patch。拿掉 reflux 之后，质量漂移的符号会翻转。

五千万网格的计算里，真正决定答案的网格有多少个？一道激波、一层剪切层、一张火焰面。通常只占全体的百分之几。剩下的网格都在光滑区域里光滑地空转。

把 n_buf 降到 0，再把 regrid every 推到 30，界面上会沿途长出红色单元。传感器走在前面，网格跟不上。

threshold0.35n_buf2max level3regrid every14

leaf 0 / 1 · under-resolved 0

set n_buf to 0 and drag “regrid every” up to 30: red cells appear along the layer every cycle and clear the instant the mesh is rebuilt — that is the feature running out of its own patch. Push n_buf back to 2 and the red stops, but watch the leaf-cell count climb. Raising the threshold thins the patch the other way.

把 n_buf 提到 2，红色消失，叶子单元数上升。加宽缓冲带的成本就是这个数。

自适应网格加密（AMR，adaptive mesh refinement · 只在解需要的地方布网格的技术）直接动的就是这个比例。决定在哪里细分的传感器，以及细分之后从界面漏掉的质量。两个要害都落成数字。

六万五千个网格里的一万一千个

先把收益数清楚。取一个 $256^2$ 的计算域，里面装一层 Kelvin–Helmholtz 剪切层（两层存在速度差的流体之间卷起的界面）。均匀网格是 65,536 个单元。用 3 级 AMR 把同一条界面包住，叶子（leaf）单元是 11,776 个。18%。

这个比例不是巧合，是维数定的。 $D$ 维计算域中，界面是 $D-1$ 维的。均匀网格单元数为 $N_{\rm uni} = (L/h)^D$ 时，只覆盖界面所需的单元数是

N_{\rm AMR} \sim \left(\frac{L}{h}\right)^{D-1} = N_{\rm uni}^{(D-1)/D}

其中 $L$ 是计算域尺度， $h$ 是最细网格间距。2D 里指数是 1/2，3D 里是 2/3。也就是说网格加密一倍时，均匀网格的单元数涨 8 倍，AMR 只涨 4 倍。分辨率越高，这笔账越划算。用 AMR 的理由在这里。

把这套网格代进估计式，数就具体了。均匀的 $256^2$ 是 65,536。后面数出来的叶子合计是 11,776， $11776/65536=0.180$ 。2D 里一条直界面的估计 $(L/h)^{1}=256$ ，只占均匀网格的 0.39%。卷起来的剪切层再加上缓冲带和粗叶子，才落到 18%。3D 里把间距减半，均匀网格是 $2^3=8$ 倍，界面一侧的 AMR 是 $2^{2}=4$ 倍。指数是 $(D-1)/D$ ，3D 为 2/3，网格越细，相对均匀网格的比例越往下掉。

决定切哪里的不是梯度

最先想到的打标判据是 $|\nabla \phi| > \epsilon$ 。它失败的原因在于量纲。压力梯度是 Pa/m，密度梯度是 kg/m⁴。换一个场、换一个问题，甚至换一个层级， $\epsilon$ 都得重新调。

Löhner 在 1987 年给出的判据用归一化消掉了这个问题：把二阶差分除以一阶差分绝对值之和。

E_i = \sqrt{ \sum_{d} \frac{\left(\phi_{i+e_d} - 2\phi_i + \phi_{i-e_d}\right)^2}{\Big(|\phi_{i+e_d}-\phi_i| + |\phi_i-\phi_{i-e_d}| + \varepsilon\big(|\phi_{i+e_d}| + 2|\phi_i| + |\phi_{i-e_d}|\big)\Big)^2} }

$e_d$ 是 $d$ 方向的邻居， $\varepsilon$ 是噪声滤波系数（一般取 0.01~~0.05）。分子分母同量纲，所以 $E_i$ 无量纲，取值大致被限制在 $[0,1]$ 内。0.3~~0.4 这一个阈值，对压力管用、对密度也管用，对 level 0 管用、对 level 3 也照样管用。

有一点要留神：这个传感器不是特征探测器，而是分辨率探测器。如果界面在该层级上已经被 4~5 个单元解开了，二阶差分就变小，传感器随之熄火。这是好性质——只细分到够用为止，然后自动停手。但如果你想要“激波一律加密到最高层级”，光靠它做不到。这时不要去调大 $\varepsilon$ 的噪声滤波项，而应该另外挂一条物理判据（例如 $\Delta p / p > 0.1$ ）做 OR。

缓冲带是买给下一次 regrid 的保险#

只细分传感器点亮的单元，下一步就会崩。因为网格不是每一步都重建的。regrid 通常每 4~20 步跑一次，这期间界面一直在移动。把打标单元向外膨胀 $n_{\rm buf}$ 个单元，理由就在这里。

需要的缓冲带宽度可以直接算出来。把层级 $\ell$ 上 $N_{\rm regrid}$ 步内特征移动的距离换算成单元数：

n_{\rm buf} \;\ge\; \frac{|u|_{\max}\,\Delta t_\ell\,N_{\rm regrid}}{h_\ell} \;=\; \nu\,N_{\rm regrid}

$\nu$ 是该层级的 CFL 数。CFL 取 0.4、每 10 步 regrid 一次，最少需要 4 个单元。这个关系与层级无关——因为用了子循环（subcycling）之后， $\Delta t_\ell$ 和 $h_\ell$ 按同样比例缩小。

写成下限就是 $n_{\rm buf}\ge 0.4\times 10=4$ 。代码里的缓冲带 2 比这条下限窄。再窄，下一次 regrid 之前界面就会走出打标区。那些红单元，把文首拨盘的 n_buf 调到 0 就能看见。

块、单元、补丁——切分单位决定什么

同样的打标结果，按什么单位去细分，会让单元数和代码复杂度分道扬镳。

方式	refine 单位	数据结构	过度 refine	代表实现
Block-based	固定大小的块（ $8^3$ 等）	块 octree	大	PARAMESH, FLASH
Cell-based	单个单元	单元级 tree	无	OpenFOAM `hexRef8`, RAGE
Patch-based	任意大小的矩形 patch	box 列表	小	Chombo, BoxLib/AMReX

块方式的数据结构最简单，缓存局部性好。代价是一个块里哪怕只有一个被打标的单元，整块都要细分。单元方式完全没有过度细分，但每次找邻居都得遍历树。OpenFOAM 属于这一类——dynamicRefineFvMesh 是引擎，hexRef8 是切刀（把一个六面体切成 8 个），refinementHistory 是用来回退的履历。只把需要的格子长出来的四叉树用的也是单元这个单位。不会因为一个标记就把整块切开。patch 方式介于两者之间，矩形内部可以照常跑结构网格循环，对向量化有利。

同一个界面上有两个答案

真正的麻烦从这里开始。看层级 $\ell$ 与 $\ell+1$ 相接的那个面。粗网格一侧它是一个面，细网格一侧是 $r^{D-1}$ 个面（ $r$ 是加密比 refinement ratio，通常为 2）。再加上用了子循环，时间步长也不一样。

\Delta t_\ell = \frac{\Delta t_0}{r^{\ell}}

粗网格走一步的工夫，细网格走 $r$ 次。于是穿过那一个面，粗网格算出的通量 $F^c$ ，和细网格在 $r^{D-1}$ 个面上分 $r$ 次实际推出去的通量，是两个不同的数。两者之差

\delta F_f = \frac{1}{r^{D}}\sum_{s=1}^{r}\sum_{f'\subset f} F^{s}_{f'} \;-\; F^{c}_{f}

这就是界面凭空造出的质量。不该有的质量冒出来，本该有的质量消失掉。哪怕你用的是守恒型格式也一样——守恒性只在单个网格内部成立，两套网格相接的地方并不保证。

记在账本上，一次性结清

解法很简单。粗步开始时把 $F^c_f$ 记进账本（flux register）。细网格跑 $r$ 次的过程中，把实际通量累加进同一本账。粗步结束时，把差额 $\delta F_f$ 退还给 patch 外侧的粗单元。

\phi_c \;\leftarrow\; \phi_c \;\mp\; \frac{\Delta t_c}{h_c}\,\delta F_f

符号由该单元位于面的哪一侧决定。patch 内侧的单元不动——那里已经用细网格值的平均覆盖过了。就这一行，把整个计算域的质量找回到机器精度。

$r=2$ 、 $D=2$ 时 $r^{D}=4$ 。细面 $r^{D-1}=2$ 个，再加上 $r=2$ 次子循环，都要除以这个 4，才能和一条粗通量对齐。少了这个除法， $\delta F$ 的量级就错了。账本记了、也还了，质量仍然对不上。

用代码数一数单元数和质量漂移

先看打标。给剪切层的一张快照套上 Löhner 传感器，加上缓冲带，再按块堆叠层级，最后数叶子单元。

import numpy as np
 
N_EFF, BLOCK, MAX_LEVEL = 256, 4, 2   # 256^2 等效分辨率, 每块 4x4 单元, 层级 0~2
EPS_L, THRESH, N_BUF = 0.02, 0.35, 2  # 噪声滤波 / 打标阈值 / 缓冲带 (单位: 单元)
 
 
def shear_layer(n, t=1.35):
    """Kelvin-Helmholtz 卷起的快照 —— 不用求解器, 只造一个场。"""
    x = (np.arange(n) + 0.5) / n
    xx, yy = np.meshgrid(x, x, indexing='ij')
    warp = 0.06 * np.sin(2 * np.pi * xx + t) + 0.025 * np.sin(4 * np.pi * xx - 2 * t)
    return np.tanh((yy - 0.5 - warp) / 0.012)
 
 
def shift(a, d, ax):
    """邻居引用: x 为周期, y 为零梯度 —— 不在壁面处造出假跳变。"""
    return np.roll(a, d, ax) if ax == 0 else np.pad(a, 1, mode='edge')[1:-1, 1 + d:a.shape[1] + 1 + d]
 
 
def lohner_sensor(f):
    """归一化的二阶差分。无量纲, 所以一个阈值通吃所有层级。"""
    e2 = np.zeros_like(f)
    for ax in (0, 1):
        p, m = shift(f, -1, ax), shift(f, 1, ax)
        num = np.abs(p - 2.0 * f + m)
        den = np.abs(p - f) + np.abs(f - m) + EPS_L * (np.abs(p) + 2 * np.abs(f) + np.abs(m))
        e2 += (num / np.maximum(den, 1e-30)) ** 2
    return np.sqrt(e2)
 
 
def grow(mask, width):
    """缓冲带: 保证下一次 regrid 之前特征跑不出 patch。"""
    for _ in range(width):
        out = mask.copy()
        for ax in (0, 1):
            out |= shift(mask, 1, ax) | shift(mask, -1, ax)
        mask = out
    return mask
 
 
def tag_blocks(f, thresh, n_buf):
    """单元传感器 -> 单元级缓冲带 -> 只要含一个打标单元的块就 refine。"""
    tagged = grow(lohner_sensor(f) > thresh, n_buf)
    nb = f.shape[0] // BLOCK
    return tagged.reshape(nb, BLOCK, nb, BLOCK).any(axis=(1, 3))
 
 
def leaf_cells(field):
    """从粗层级逐级向下扫, 只数实际参与求解的单元。"""
    counts, live = [], None
    for lev in range(MAX_LEVEL + 1):
        n = N_EFF >> (MAX_LEVEL - lev)
        f = field.reshape(n, N_EFF // n, n, N_EFF // n).mean(axis=(1, 3))
        nb = n // BLOCK
        child = np.zeros((nb, nb), bool) if lev == MAX_LEVEL else tag_blocks(f, THRESH, N_BUF)
        live = np.ones((nb, nb), bool) if live is None else live
        counts.append(int((live & ~child).sum()) * BLOCK * BLOCK)
        live = np.kron(live & child, np.ones((2, 2), bool))
    return counts
 
 
counts = leaf_cells(shear_layer(N_EFF))
total, uniform = sum(counts), N_EFF * N_EFF
for lev, c in enumerate(counts):
    print(f'  level {lev}  h = 1/{N_EFF >> (MAX_LEVEL - lev):<3d}  leaf cells = {c:6d}')
print(f'  AMR total     = {total}')
print(f'  uniform 256^2 = {uniform}  ->  {100 * total / uniform:.1f} % of the cells')

输出是这样的。

  level 0  h = 1/64   leaf cells =   3280
  level 1  h = 1/128  leaf cells =   1520
  level 2  h = 1/256  leaf cells =   6976
  AMR total     = 11776
  uniform 256^2 = 65536  ->  18.0 % of the cells

把 N_BUF 依次改成 0、1、2、4，总数会走成 9,616 → 10,864 → 11,776 → 14,224。两格缓冲带的标价是 2,160 个单元，相对均匀网格是 3.3 个百分点。

从 0 加到 2 的 2,160 个，大约是 9,616 的 22%。加到 4 格是 14,224，比没有缓冲带多 4,608 个单元。按层级看，最细的 6,976 占全部 11,776 的 59%。第 0 层的 3,280 是用粗间距盖住光滑区域的那一份。18% 不是界面单独的数，而是这三层加在一起。

接着看重通量。在 2D 标量对流上叠两个层级，加进子循环，然后称一称总质量。

import numpy as np
 
NC, R, CFL, NSTEP = 48, 2, 0.4, 60   # 粗网格 / 加密比 / CFL / 粗步数
BOX = (12, 28, 16, 32)               # patch 边角 (按粗网格索引)
U, V = 1.0, 0.6
 
 
def upwind_faces(f, dx, dy):
    """周期块所有 x·y 面上的 donor-cell 通量。返回 (Fx, Fy)。"""
    fx = U * (f if U > 0 else np.roll(f, -1, 0))          # 面 i 是单元 i 的左侧
    fy = V * (f if V > 0 else np.roll(f, -1, 1))
    return np.roll(fx, 1, 0), np.roll(fy, 1, 1)
 
 
def march_block(f, fx, fy, dt, dx, dy):
    return f - dt / dx * (np.roll(fx, -1, 0) - fx) - dt / dy * (np.roll(fy, -1, 1) - fy)
 
 
def gaussian_patch(n, x0, y0, s):
    c = (np.arange(n) + 0.5) / n
    xx, yy = np.meshgrid(c, c, indexing='ij')
    return np.exp(-((xx - x0) ** 2 + (yy - y0) ** 2) / s ** 2)
 
 
def two_level_run(reflux):
    i0, i1, j0, j1 = BOX
    dx, dxf = 1.0 / NC, 1.0 / NC / R
    dt = CFL * dx / (abs(U) + abs(V))
    dtf = dt / R
 
    coarse = gaussian_patch(NC, 0.32, 0.42, 0.09)
    fine = np.kron(coarse[i0:i1, j0:j1], np.ones((R, R)))   # patch 起始时与粗网格值保持一致
 
    for _ in range(NSTEP):
        cfx, cfy = upwind_faces(coarse, dx, dx)
        # 粗网格“以为”自己穿过 patch 边界的通量
        edge_c = {'lo_x': cfx[i0, j0:j1].copy(), 'hi_x': cfx[i1, j0:j1].copy(),
                  'lo_y': cfy[i0:i1, j0].copy(), 'hi_y': cfy[i0:i1, j1].copy()}
        coarse = march_block(coarse, cfx, cfy, dt, dx, dx)
 
        edge_f = {k: np.zeros_like(v) for k, v in edge_c.items()}
        for _ in range(R):                                   # 子循环: 粗网格 1 步对应细网格 R 步
            g = np.zeros((fine.shape[0] + 2, fine.shape[1] + 2))
            g[1:-1, 1:-1] = fine
            g[0, 1:-1] = np.repeat(coarse[i0 - 1, j0:j1], R)  # ghost: 把粗值以分段常数注入
            g[-1, 1:-1] = np.repeat(coarse[i1, j0:j1], R)
            g[1:-1, 0] = np.repeat(coarse[i0:i1, j0 - 1], R)
            g[1:-1, -1] = np.repeat(coarse[i0:i1, j1], R)
            gfx, gfy = upwind_faces(g, dxf, dxf)
            fine = march_block(g, gfx, gfy, dtf, dxf, dxf)[1:-1, 1:-1]
            for k, s in (('lo_x', gfx[1, 1:-1]), ('hi_x', gfx[-1, 1:-1]),
                         ('lo_y', gfy[1:-1, 1]), ('hi_y', gfy[1:-1, -1])):
                edge_f[k] += s.reshape(-1, R).mean(axis=1) / R   # 沿面方向与时间方向取平均
 
        coarse[i0:i1, j0:j1] = fine.reshape(i1 - i0, R, j1 - j0, R).mean(axis=(1, 3))
 
        if reflux:                                            # 账本结算
            coarse[i0 - 1, j0:j1] -= dt / dx * (edge_f['lo_x'] - edge_c['lo_x'])
            coarse[i1, j0:j1] += dt / dx * (edge_f['hi_x'] - edge_c['hi_x'])
            coarse[i0:i1, j0 - 1] -= dt / dx * (edge_f['lo_y'] - edge_c['lo_y'])
            coarse[i0:i1, j1] += dt / dx * (edge_f['hi_y'] - edge_c['hi_y'])
 
    mask = np.ones((NC, NC), bool)
    mask[i0:i1, j0:j1] = False
    return (coarse * mask).sum() * dx * dx + fine.sum() * dxf * dxf
 
 
m0 = gaussian_patch(NC, 0.32, 0.42, 0.09).sum() / NC ** 2
for tag, on in (('reflux off', False), ('reflux on ', True)):
    m = two_level_run(on)
    print(f'  {tag}:  mass = {m:.12f}   drift = {(m - m0) / m0:+.3e}')

  reflux off:  mass = 0.025702960114   drift = +1.006e-02
  reflux on :  mass = 0.025446894886   drift = +0.000e+00

把 pulse sigma 调小，红色偏离变大，账本里的 $\delta F$ 也按同样的幅度变大。

ratio r4CFL0.60pulse sigma0.050

step 0 · drift 0.00e+0 vs 0.00e+0

the pulse is harmless until it reaches the pink face. From that step on the red trace walks away from zero and never comes back, while the green one stays pinned at machine zero. Sharpen the pulse (sigma down) or raise r and the red excursion grows — the register grows with it, because it is exactly the same number with the opposite sign.

脉冲碰到粉色面之后，只有红色曲线离开 0，绿色贴着 0。关掉重通量时的 drift +1.006e-02 就是那条红线。

60 步就是 1%。把输出里的两个质量相减，差值大约是 $2.56\times 10^{-4}$ 。再用 on 一侧的质量去除，就和 drift +1.006e-02 那个 1% 对上。而且这个值并不单调增长。第 15 步是 −1.3%，第 30 步是 −0.77%，第 60 步是 +1.0%。团块每次进出 patch，符号就翻一次。跑收敛性测试时看到这种曲线，人往往会怀疑格式本身，真凶其实是没做重通量的那个界面。

负载均衡——用一条曲线把块排成一队

AMR 上并行以后会冒出新问题。每次 regrid，各进程手里的块数都不一样。抱着界面的那个 rank 单元数暴涨 5 倍，只有光滑区域的 rank 原地不动。

标准处方是空间填充曲线（space-filling curve, SFC · 用一条线把多维网格串起来遍历的曲线）。用 Morton（Z-order）或 Hilbert 曲线给所有叶子块编上一维索引，再把这条队伍按 rank 数均分。靠曲线的局部性，相邻索引大体上在空间上也相邻，于是均分同时就是通信量很小的划分。p4est 和 AMReX 走的就是这条路。用图划分器（ParMETIS、Zoltan、Scotch）划分质量更好，但每次 regrid 都得重跑一遍，开销比 SFC 大得多。AMR 属于网格频繁变化的一类，所以大多数情况下 SFC 赢。

拿掉重通量，收敛曲线会把格式当成凶手

打标保持无量纲。归一化的二阶差分让一个阈值跨过场和层级。原始梯度得逐层重调。缓冲带按 $\nu N_{\rm regrid}$ 来算。CFL 0.4、10 步，就是 4 个单元。凭感觉填 1 个单元、却每 20 步才 regrid，很大一部分步数是在特征已经漏出 patch 的情况下跑完的。

把重通量当成可选项，coarse-fine 界面就会造出质量。第 60 步的 +1% 还不单调。符号来回翻的 drift 很容易被算到对流格式头上。把 refine 的次序挪了一行、误差就变大的 LBM 时间推进是同一类。次序错了，单元切得再细，界面上的答案仍然是错的。