引言
激波间断、裂纹尖端应力奇异性、焊接熔池的陡温度梯度、天体物理中的引力坍缩——多尺度问题的共同特征是：需要高分辨率的区域往往只占计算域的百分之一甚至千分之一。均匀加密网格在这种问题上极其浪费：把全局网格加密到局部所需分辨率，单元数可能增加数百倍，而真正”用到”这些单元的区域不到零头。自适应网格细化（Adaptive Mesh Refinement, AMR）让网格疏密跟随解的特征动态调整，是求解多尺度问题的使能技术。本文梳理 AMR 的误差指示器、块结构与单元树两类实现路线、粗-细界面处理与并行负载均衡，并给出一个稳态边界层问题完整 solve–mark–refine 自适应循环的 Python 实现。
理论基础
后验误差估计与细化判据
AMR 的第一步是回答”哪里需要加密”。常用判据按成本从低到高有三类：
- 特征量指示器：对物理量的梯度或二阶差分设阈值，如 $\eta_i = h_i,|\nabla u|_i$ 或 $\eta_i = |u_{i+1} - 2u_i + u_{i-1}|$。实现最简单，对流体力学激波、温度前沿等”可见”特征非常有效，但它度量的是解的变化剧烈程度而非真实误差；
- Richardson 外推估计：在同一区域用 $h$ 与 $h/2$ 两套网格各算一次，对 $p$ 阶格式有
$$e_h \approx \frac{u_h - u_{h/2}}{2^p - 1}$$
这给出定量的局部误差，代价是约 1/8 的额外计算；
- 残差型/对偶加权估计：有限元框架下由单元残差与通量跳跃加权得到严格误差界，理论最完善但实现最重。
得到误差分布后，常用 Dörfler 标记策略：按 $\eta_i$ 降序标记最少的单元，使其误差贡献之和占总量的固定比例（如 50%），兼顾精度与单元数控制。
块结构 AMR 与单元树 AMR
工程级 AMR 有两条实现路线。Berger-Colella 块结构 AMR 把标记单元聚类为若干矩形（块），块按固定细化率 $r = 2$ 或 $4$ 整体加密并层层嵌套，要求细层被粗层”恰当嵌套”（proper nesting），层间留出过渡带；单元树（cell-based）AMR 则对每个单元独立二分（一维）、四分（二维 quadtree）或八分（三维 octree），并施加 2:1 平衡约束——任意相邻单元的层级差不超过 1，保证网格质量。两条路线的对比如下：
| 维度 | 块结构 AMR | 单元树 AMR |
|---|---|---|
| 数据结构 | 矩形块集合，块内规则 | 树（quadtree/octree） |
| 求解器兼容 | 可直接复用结构网格求解器 | 需非结构/树网格求解器 |
| 细化贴合度 | 块内含冗余单元（聚类损失） | 最贴合误差分布 |
| 缓存效率 | 高（块内连续存储） | 低（指针跳转） |
| 典型软件 | AMReX、SAMRAI | p4est、deal.II |
粗-细界面与时间推进
AMR 的数值难点集中在界面。粗细单元交界处（hanging node 或 coarse-fine boundary）需要两组操作：粗到细的延拓（prolongation，插值生成细网格边界值）与细到粗的限制（restriction，把细网格解平均回粗网格）。对守恒律还需回流修正（refluxing）：用细网格界面的高精度通量替换粗网格通量，
$$F^{\mathrm{coarse}}_{i+1/2} \leftarrow \frac{1}{r}\sum_{j=1}^{r} F^{\mathrm{fi\ne}}_{j}$$
否则全局守恒性会被破坏。时间推进有两种策略：全层同步（时间步由最细层决定，粗层被迫小步推进）与分层子循环（subcycling，细层以 $\Delta t/r$ 走 $r$ 步追赶粗层），后者节省约一半计算量但要求层间时间插值。大规模并行时，各处理器持有的细单元数差异悬殊，常用空间填充曲线（Morton/Hilbert 序）剖分计算域，把几何上邻近的单元分给同一进程，兼顾负载均衡与通信局部性。
数值方法与实现
以下代码实现一个完整的 solve–mark–refine 自适应循环：求解稳态对流-扩散方程 $a u' - \varepsilon u'' = 0$（Péclet 数 $Pe = aL/\varepsilon = 100$，解在 $x=1$ 处形成宽度约 $\varepsilon/a = 0.01$ 的指数边界层），在单元中心有限体积网格上以迎风通量处理对流项，以单元相邻界面落差的最大值作为误差指示器，逐轮二分加密（最多 6 级，局部分辨率提高 $2^6 = 64$ 倍），并与均匀网格对比：
import numpy as np
def exact(x, Pe=100.0):
"""精确解: a u' - eps u'' = 0, u(0)=1, u(1)=0（Pe = aL/eps）"""
return (1.0 - np.exp(Pe * (x - 1.0))) / (1.0 - np.exp(-Pe))
def solve_on_mesh(leaves, Pe=100.0, a=1.0, uL=1.0, uR=0.0):
"""单元中心有限体积：迎风对流 + 中心扩散，Dirichlet 边界"""
eps = a / Pe
xc = leaves.mean(axis=1)
h = leaves[:, 1] - leaves[:, 0]
n = len(leaves)
A = np.zeros((n, n)); b = np.zeros(n)
for i in range(n):
if i < n - 1: # 右界面 i+1/2
d = 0.5 * (h[i] + h[i+1])
A[i, i] += a + eps / d; A[i, i+1] -= eps / d
else: # 右边界 Dirichlet
A[i, i] += a + 2*eps/h[i]; b[i] += 2*eps/h[i]*uR
if i > 0: # 左界面 i-1/2
d = 0.5 * (h[i] + h[i-1])
A[i, i] += eps / d; A[i, i-1] -= (a + eps / d)
else: # 左边界 Dirichlet
A[i, i] += 2*eps/h[i]; b[i] += (a*uL + 2*eps/h[i]*uL)
return xc, h, np.linalg.solve(A, b)
def amr_loop(base=32, max_level=6, thresh=0.02, niter=12):
"""solve-mark-refine 自适应循环：求解 -> 标记 -> 二分细化"""
leaves = np.array([[i/base, (i+1)/base] for i in range(base)])
for _ in range(niter):
xc, h, u = solve_on_mesh(leaves) # solve
du = np.abs(np.diff(u)) # mark: 界面落差
eta = np.zeros(len(u))
eta[:-1] = np.maximum(eta[:-1], du)
eta[1:] = np.maximum(eta[1:], du) # 单元取相邻界面最大值
lev = np.round(np.log2((1.0/base) / h)).astype(int)
mark = [int(i) for i in np.nonzero(eta > thresh)[0] if lev[i] < max_level]
if len(mark) == 0:
break
mset = set(mark) # refine: 二分
new = [list(l) for j, l in enumerate(leaves) if j not in mset]
for i in mark:
xl, xr = leaves[i]; xm = 0.5*(xl + xr)
new.append([xl, xm]); new.append([xm, xr])
new.sort(key=lambda e: e[0])
leaves = np.array(new)
return leaves
leaves = amr_loop()
xc, h, u = solve_on_mesh(leaves)
h_min = (leaves[:, 1] - leaves[:, 0]).min()
l1_amr = np.sum(np.abs(u - exact(xc)) * h)
print(f"AMR: {len(leaves)} 单元, 最细 h = {h_min:.5f}, L1 误差 {l1_amr:.5f}")
for n in (32, 512, 2048):
un = np.array([[i/n, (i+1)/n] for i in range(n)])
xu, hu, uu = solve_on_mesh(un)
print(f"均匀 {n:4d} 单元: L1 误差 {np.sum(np.abs(uu - exact(xu)) * hu):.5f}")运行结果非常说明问题：自适应循环从 32 个粗单元出发，指示器自动把加密集中到边界层所在的 $x \in [0.9,, 1]$ 窄带内，最终 85 个单元达到 $L_1$ 误差 $3.8\times10^{-4}$——优于 512 均匀网格的 $9.1\times10^{-4}$，逼近 2048 均匀网格的 $2.4\times10^{-4}$，而单元数仅为后者的 1/24。误差残余主要来自对流项的一阶迎风精度（数值黏性把边界层额外抹宽约 $h \cdot Pe/2$），这也是工业代码在层内改用高阶通量的动机。教学实现还省略了两个工程要点：只细化不粗化（动态问题需按滞后阈值对离开特征区的单元做 coarsening），以及未显式施加 2:1 平衡约束与守恒回流修正。
工程案例分析
问题：激光选区熔化熔池温度场模拟
某 316L 不锈钢选区熔化工艺，激光光斑直径 80 μm、扫描速度 800 mm/s，熔池温度梯度达 $10^7$ K/m。对比三种网格策略（以 1 μm 均匀网格全场结果为基准）：
| 网格策略 | 单元数（峰值） | 峰值温度偏差 | 相对计算成本 |
|---|---|---|---|
| 均匀粗网格（16 μm） | 4.8×10⁵ | 12% | 1 |
| 均匀细网格（1 μm） | 1.2×10⁸ | 基准 | ≈250 |
| 动态 AMR（4 级，跟踪热源） | 3.6×10⁶ | 1.5% | ≈9 |
AMR 方案把细化区域限制在光斑周围约 0.5 mm 的移动窗口内，误差指示器取温度梯度与热流密度的组合量，每 20 个时间步重构一次网格；配合 2:1 平衡约束与界面回流修正，热输入全局守恒误差低于 0.3%。以约 1/28 的成本取得 1.5% 的峰值温度精度，使参数扫描与工艺窗口优化在工程时间内可行。需要注意的是动态 AMR 的网格重构本身有约 5% 的开销，热源移动越快、重构越频繁，越要采用增量式重构（仅局部重划分）而非全量重建。
结论与展望
AMR 的核心逻辑是”让网格分辨率跟随误差分布”，其完整技术栈包括：后验误差指示器决定哪里加密、块结构或单元树数据结构组织层级、延拓/限制/回流修正保证界面精度与守恒、子循环时间推进与空间填充曲线负载均衡支撑大规模并行。实践建议是：特征明显的问题先用梯度指示器加单元树快速见效；守恒严格的问题（爆炸、天体物理）必须上回流修正；并行规模超过数百进程后，负载均衡与重网格开销会成为新的瓶颈。当前前沿是把 AMR 与机器学习结合——用训练好的代理模型预测下一时刻的细化区域，或用强化学习在线调优标记阈值，进一步压缩重网格成本。
发表回复