四叉树 AMR — 让网格只在需要的地方生长

1984年，斯坦福。Marsha Berger 和 Joseph Oliger 在追一道激波。网格大约铺了一百万个单元。真正在做事的，只有激波附近那约1%。其余99%用同样的精度反复重算均匀自由流。当年登在 Journal of Computational Physics 上的论文，就是 Berger–Oliger AMR。同样的习惯现在还留在 OpenFOAM 的 dynamicRefineFvMesh、AMReX、Chombo 和 p4est 里。

拖动白点，网格只在点的周围裂开；threshold 滑块改的是决定是否切开的缩放梯度 $\tau$。

$\tau$ 降到 0.05，叶子会追到高斯尾部；推过 0.4，就只剩近乎均匀的粗网格。按下 add feature 再放一个峰，深处只剩两个核，左上角的读数是相对均匀网格的单元数比例。

均匀网格为什么浪费 90%#

跑过三维可压 RANS 的人熟悉这笔交换。把网格间距放大一倍，单步开销变成八分之一，因为 $2^3=8$。激波、剪切层、分离的梯度也一起被抹平。若为了那条窄带把整个区域都加密，自由流的单元也跟着变小。内存和时间一起爆。

只加密那条带就够了。比例是下面这个估计。

$L$ 是最大加密级别。$V_\text{feature}$ 是需要细网格的体积，$V_\text{domain}$ 是整个区域。$4^{L}$ 的意思是，二维里每深一级，单元数乘四。$L=5$ 时 $4^5=1024$。把开头那个1%当作体积分数，相对最粗网格的单元数是 $0.01\times 1024+0.99=11.23$ 倍。同一级的均匀网格则是粗网格的1024倍。两者相除，AMR 大约只占均匀细网格的1.1%。激波若更薄，$V_\text{feature}/V_\text{domain}\sim h$，这个份额还会再降。所谓个位数百分点，就是这次除法。三维里只把一条边减半，单元数也是8倍。把这8倍铺到整个区域，账单覆盖全域；只铺在窄带里，8倍就停在窄带上。

h-、p-、r- — 自适应的三条路#

自适应网格（adaptive mesh，随解的分布改变单元大小或阶数的网格）按动的是什么来分。

• h-自适应: 缩小单元尺寸 $h$。最常见。说 AMR，通常就是这一路。
• p-自适应: 单元不动，只提高多项式阶 $p$。DGM 和 SEM（谱单元法）里多见。
• r-自适应（relocation）: 单元数和阶数都固定，只把节点挪向梯度大的地方。内存不变。网格畸变还在。

这里只跟 h-AMR，尤其是 二维四叉树。三维里同一想法是八叉树。

四叉树 — 单元拥有自己孩子的数据结构

均匀网格一个数组就够。自适应网格需要 树。一个单元 refine 之后生出四个孩子，孩子指向父亲。孩子再裂，就是孙子。

class QuadCell:
    __slots__ = ("cx", "cy", "size", "level", "children", "parent", "data")
 
    def __init__(self, cx, cy, size, level=0, parent=None):
        self.cx = cx           # 单元中心 x
        self.cy = cy           # 单元中心 y
        self.size = size       # 边长
        self.level = level     # 根 = 0,加密一次 +1
        self.children = None   # 叶子时为 None,加密后为长度 4 的列表
        self.parent = parent
        self.data = 0.0        # 守恒量(如密度)
 
    def is_leaf(self):
        return self.children is None
 
    def split(self):
        h = self.size * 0.5
        q = self.size * 0.25
        offsets = ((-q, -q), (+q, -q), (-q, +q), (+q, +q))
        self.children = [
            QuadCell(self.cx + dx, self.cy + dy, h, self.level + 1, self)
            for dx, dy in offsets
        ]
        for ch in self.children:
            ch.data = self.data    # 用父值初始化(守恒)

split 里用父值初始化孩子，是为了 守恒。父单元的积分 $U \cdot V$ 必须等于四个孩子之和。一个孩子的体积是父亲的四分之一。原样复制时，四个孩子的和是 $4 \times (U \cdot V/4) = U \cdot V$，和父亲相同。这种复制只是一阶插值。生产代码通常再叠上父亲的斜率（slope），做成 MUSCL 型的二阶 prolongation。斜率可以加上。四个孩子的和仍然得等于父亲的积分。和一旦对不上，裂开的单元就在制造或吃掉质量。

切在哪里 — 加密判据

四叉树是空盒子。切哪一个单元，才是加密判据。生产里常见五种。

OpenFOAM 的 dynamicRefineFvMesh 对用户选定的场设阈值。AMReX 把每个单元的布尔函数 tagcells() 留空，由调用者去填。标记打在哪里，树的形状就跟着变。这比数据结构本身更能决定结果。把单元数压下去仍保住同一答案的标记方式，接在标记与 coarse-fine 回流里。

浏览器里的网格只用表的第一行。单元中心的中心差分 $|\nabla u|$ 再乘单元尺寸 $h$。乘积超过 $\tau$ 就 split，否则留作 leaf。

2:1 平衡与 hanging node#

孩子和邻居相遇的面会出问题。一边是大单元，另一边是两个小单元，这个面上就出现 hanging node（悬空节点）。在有限体积里，这个节点会破坏通量平衡。

出路有两条。

1. 强制 2:1 平衡: 把相邻 leaf 的级别差限制在 1 以内。p4est、AMReX、Berger 一系都用这条。
2. non-conforming 插值: 在 hanging face 上另做插值，用两个小单元的值合成大单元的面值。DGM 一系里常见。

强制的循环很短。拿到加密候选，若邻居里有低两级以上的 leaf，就把该邻居也压进队列。直到队列空掉。

def enforce_2to1(roots, refine_queue):
    """确保 refine_queue 中所有单元与邻居满足 2:1"""
    pending = list(refine_queue)
    queued = set(id(c) for c in pending)
    while pending:
        cell = pending.pop()
        for nb in face_neighbors(cell, roots):
            if nb.is_leaf() and cell.level - nb.level >= 1:
                if id(nb) not in queued:
                    pending.append(nb)
                    queued.add(id(nb))
    return [c for c in (refine_queue + pending) if c.is_leaf()]

face_neighbors 返回共享同一面的 leaf。走的是树，复杂度 $O(\log N)$。

用 Python 看 1D 自适应网格#

把完整的二维四叉树写出来会很长。同一个想法放在 一维二叉树 上走一步就够。放上高斯脉冲，只对超过阈值的区间二等分。

import numpy as np
 
def gaussian_pulse(x, x0=0.6, sigma=0.04):
    return np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
 
def refine_indicator(cell_lo, cell_hi, fn, h_scale=1.0):
    """缩放一阶导 = 单元两端值的差的绝对值"""
    u_lo = fn(cell_lo)
    u_hi = fn(cell_hi)
    return abs(u_hi - u_lo) * h_scale
 
def adaptive_bisect_1d(fn, tau=0.05, max_level=8):
    """按阈值递归二等分根 [0, 1]"""
    leaves = []
    def recurse(lo, hi, level):
        ind = refine_indicator(lo, hi, fn)
        if level >= max_level or ind < tau:
            leaves.append((lo, hi, level))
            return
        mid = 0.5 * (lo + hi)
        recurse(lo, mid, level + 1)
        recurse(mid, hi, level + 1)
    recurse(0.0, 1.0, 0)
    return leaves
 
cells = adaptive_bisect_1d(gaussian_pulse, tau=0.05, max_level=7)
print(f"leaf cells: {len(cells)}")
print(f"对均匀第 7 级: {len(cells) / 2**7 * 100:.1f}%")
# 输出示例: leaf cells: 41, 对均匀第 7 级: 32.0%

一个高斯只占区域里很窄的一条。leaf 是 41 个，相对均匀第 7 级是 32.0%。两个脉冲时，深处只有两个核，中间保持粗糙。

41 个叶子和 128 个单元之间#

把输出的 41 除以 $2^7=128$，得到 0.3203。注释里的 32.0% 就是这个数收成一位小数的百分比。128 减去 41 是 87。若均匀铺到第 7 级，大约 68% 的单元不会被做成那么小。代码里 $\sigma=0.04$，只占区间 $[0,1]$ 长度的 4%。梯度超过 $\tau=0.05$ 的地方局限在这附近，平坦自由流两端的值几乎没有差，二等分会早早停下。一维切一次，孩子是两个。二维四叉树是四个，所以同样窄的带子，单元数之比会更陡地跟着前面的 $4^{L}$ 走。把 $\tau$ 抬到 0.4，原先落在 0.05 和 0.4 之间的尾部就不再继续裂。

并行陷阱 — load balancing 与 SFC#

分布式内存还会撞上另一堵墙。若某个 rank 抱住所有深树，这个 rank 就决定整个 wall clock。单元数的平衡（load balancing）已经破了。

把树排成一维再切开即可。用 空间填充曲线（space-filling curve, SFC） — Morton（Z-order）或 Hilbert — 做索引，空间上相邻的单元在一维里也保持接近。p4est 的 Morton 索引是把坐标的比特交错插进去。

def morton_2d(ix, iy, level):
    """把整数坐标 (ix, iy) 转为 Morton(Z-order)索引"""
    code = 0
    for b in range(level):
        code |= ((ix >> b) & 1) << (2 * b)
        code |= ((iy >> b) & 1) << (2 * b + 1)
    return code
 
# 4x4 网格的 Z-order
for iy in range(4):
    for ix in range(4):
        print(morton_2d(ix, iy, 2), end=" ")
    print()
# 0  1  4  5
# 2  3  6  7
# 8  9 12 13
# 10 11 14 15

按 Morton 索引排序，再按 rank 数切开，相邻单元落在同一个 rank，通信量就下来了。正方形区域的角上仍会出现跳跃。Hilbert 曲线避开这种跳跃，但索引更贵。实务里 Morton 仍是事实上的标准。

标错场的话，32% 也仍然贵#

树的形状差不多已经定了。孩子四个，再用 Morton 切开，邻居就会进同一个 rank。4×4 上顺序从 0、1、4、5 开始，是因为坐标比特被逐位插进索引。钱漏在下一步。标记打在密度上，误差却在剪切层里，那么即使下到第 7 级，41 个 leaf 也会堆在错误的带子上。32% 属于基准对准 $x_0=0.6$、$\sigma=0.04$ 的时候。基准场一换，式子里的 $V_\text{feature}$ 就换，相对粗网格的 11.23 倍也跟着动。

也可以选择丢掉 2:1。把邻居级别差限在 1，hanging node 在一个面上最多对着两个小单元。放开这个限制，每个面都得单独写 non-conforming 插值。面上的通量容易和父亲的积分 $U\cdot V$ 对不齐。空间切完以后，refine 放在时间推进的哪一行，仍然没定。那一行能把误差放大到什么程度，写在LBM AMR 的递归时间推进里。