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

Löhner 传感器、缓冲带宽度，以及 flux register 堵住的界面质量泄漏

五千万网格的计算里，真正决定答案的网格有多少个？一道激波、一层剪切层、一张火焰面。数一数，通常只占全体的百分之几。剩下的网格都在光滑区域里光滑地空转。自适应网格加密（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 的唯一理由。

决定切哪里的不是梯度

最先想到的打标判据是 $|\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$ 按同样比例缩小。

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

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 降到 0、再把 regrid every 推到 30，沿着界面就会长出红色单元。这些是传感器在喊“这里该细分”、而网格还没跟上的单元。把 n_buf 提到 2，红色消失，代价是叶子单元数上升。加宽缓冲带的成本，就是这个数字。

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

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

方式	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 内侧的单元不动——那里已经用细网格值的平均覆盖过了。就这一行，把整个计算域的质量找回到机器精度。

下面把开启重通量（refluxing）和关闭重通量的两侧并排跑一遍。

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。把 pulse sigma 调小让脉冲变尖，红色的偏离幅度会变大——而记进账本的 $\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 个百分点。

接着看重通量。在 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

60 步就是 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 赢。

打开 AMR 之前要定的三件事#

第一。 打标判据要做成无量纲的。用归一化的二阶差分，一个阈值就能通吃所有场、所有层级。用原始梯度，就得逐层调参。

第二。 缓冲带宽度按 $\nu N_{\rm regrid}$ 算出来再填。凭感觉填 1 个单元、却每 20 步才 regrid 一次，那么半程计算都是在特征已经漏出 patch 的状态下跑的。

第三。 重通量不是选项。哪怕用的是守恒型格式，coarse-fine 界面也照样破坏守恒。在燃烧、多相流这类质量本身就是答案的问题里，没有 flux register 跑出来的收敛曲线不可信。