4. 代码验证的方法与实践

代码验证回答的是一个很具体的问题：我写的求解器，是否按我声称的精度在解我写的方程？它不关心物理模型选得对不对，只关心离散方程、时间推进、线性求解和程序实现是否自洽。因此代码验证的核心工具不是实验数据，而是可精确控制的数学解与可量化的误差收敛率。本节围绕网格收敛指数（GCI）、制造解方法（MMS）、精确解对比、收敛阶验证、时间步收敛性和迭代收敛容差设定展开，并给出一个可复现的实操案例。

4.1 收敛阶验证：代码验证的主线

设离散误差 e 随网格尺度 h 按 e ≈ C h^p 衰减，其中 p 为观测收敛阶。若格式理论精度为 p_t，则代码验证的基本判据是 p ≈ p_t。常用做法是取三套系统性加密的网格 h₁ < h₂ < h₃，加密比 r = h₂/h₁ = h₃/h₂，用同一离散格式计算同一问题，得到解 f₁、f₂、f₃，则

p = ln[(f3 - f2)/(f2 - f1)] / ln(r)

这里 f 可以是某个关注量，如中心点温度、出口平均速度、升力系数或积分热流。注意：p 只有在进入渐近收敛区后才接近理论值；若三套网格都太粗，或误差被其他误差源污染，p 会明显偏离。因此收敛阶验证不是一次计算，而是“加密—观察—再加密”的过程。

4.2 网格收敛指数 GCI

GCI 由 Roache 提出，用于估计细网格解相对于“网格无关解”的离散误差带。在得到 p 后，用 Richardson 外推估计精确解：

f_exact ≈ f1 + (f1 - f2)/(r^p - 1)

细网格上的相对误差估计为

GCI_fine = Fs * |(f2 - f1)/f1| / (r^p - 1)

其中 F_s 为安全因子：当三套网格均处于渐近区且 p 与理论值接近时取 1.25；若 p 不可靠或网格不足，取 3.0。工程上常报告“f₁ ± GCI”，例如“出口平均温度 352.4 K，GCI = 0.8%”。若 GCI 大于可接受阈值（如 5%），说明网格仍需加密，或该关注量对网格过于敏感。

4.3 制造解方法 MMS

精确解对比的困难在于：有解析解的方程往往太简单，无法覆盖真实代码中的变系数、源项、边界条件和多物理耦合。MMS 的思路是反过来：先人为指定一个光滑的制造解，例如

T(x,y) = sin(πx) sin(πy) + 0.3 x^2 y

把它代入控制方程，得到非零的源项 S；再把制造解代入边界条件，得到 Dirichlet/Neumann 边界值。然后修改代码，加入该源项并使用该边界条件，运行后与制造解逐点比较。由于制造解是已知的，误差可以精确计算，收敛阶可以严格验证。MMS 的关键操作是：选择足够光滑的制造解（保证各阶导数存在）、推导源项、实现源项与边界、在多个网格上计算误差范数、拟合 p。它特别适合验证变系数扩散、对流项、源项和耦合项的实现是否正确。

4.4 时间步收敛性与迭代容差

瞬态问题中，总误差由空间离散误差、时间离散误差和迭代误差叠加。时间步收敛性验证的做法与网格收敛类似：固定空间网格，取三套时间步 Δt₁ < Δt₂ < Δt₃，计算同一时刻的关注量，用相同公式估计时间收敛阶。若时间格式为一阶隐式，理论 p_t = 1；二阶 Crank-Nicolson 或 BDF2 则 p_t = 2。若观测阶远低于理论值，常见原因是时间步太大、非线性迭代未收敛或空间误差主导。

迭代收敛容差设定需要区分三层：线性求解器容差、非线性迭代容差和时间步内迭代容差。原则是：迭代误差应比离散误差小一个数量级以上，否则会污染收敛阶。例如，若关注量约为 350 K，目标离散误差约 0.1 K，则非线性残差容差可设为 1e-6～1e-8（相对残差），线性求解器容差再紧 1～2 个数量级。容差过松会导致 p 虚高或虚低，容差过紧则浪费计算时间。建议做一次“容差敏感性检查”：把容差收紧 10 倍，若关注量变化远小于 GCI，则当前容差足够。

4.5 代码验证案例实操

以一个二维稳态导热求解器为例，验证其扩散项离散是否达到二阶精度。取制造解 T = sin(πx)sin(πy)，计算域 [0,1]×[0,1]，热导率 k = 1，源项由 -∇·(k∇T) = S 推出 S = 2π² sin(πx)sin(πy)。边界取制造解值。采用五点差分格式，取 N = 20、40、80 三套网格，计算中心点温度 T(0.5,0.5) 与 L2 误差范数。

表中数据为教学示例，用于说明判读方式。观测阶 p ≈ 2.00，与五点差分理论二阶精度一致，说明扩散项实现正确。进一步用 GCI 估计细网格误差：取 r = 2、p = 2、F_s = 1.25，得 GCI ≈ 1.25 × |(7.8e-4 - 1.95e-4)/1.95e-4| / (2² - 1) ≈ 0.42%，说明 80×80 网格的中心点温度已具有约 0.4% 的离散误差带。

把上述流程固化为回归测试后，代码验证就从“一次性检查”变成“持续保障”。每次修改求解器，先跑 MMS 看 p 是否仍接近理论值，再看 GCI 是否稳定；若 p 下降或 GCI 突增，优先检查最近改动的离散项、边界实现和迭代容差，而不是急于调物理参数。这样，代码验证才能真正为后续的模型确认提供可信的数值基础。