MMS源项的自动推导——故障排查指南

MMS源项自动导出（符号计算） — 故障排除指南

源项的bug与"求解器的bug"无法区分

前提整理：用数学式处理系统（SymPy等）从制造解 \( u_m \) 与算子 \( L \) 推导源项 \( s = L(u_m) \)，以生成代码（C/Fortran/UDF）形式注入求解器——查出这条链上的缺陷是本页的范围。推导的基础流程见整合版。

第一件事——随机点上与数值微分交叉核对

推导出的源项，可以用与符号微分独立的数值方法确认。在域内取约10个随机点 \( \mathbf{x}_i \)，比较①生成代码给出的 \( s(\mathbf{x}_i) \) 与②对制造解施加算子的高精度有限差分数值评估。中心差分步长取 \( h \approx 10^{-5} \) 时，相对误差应在 \( 10^{-8} \) 量级一致——哪怕一个点超差，就说明推导或代码生成有误。这个核对几分钟就能写好，几乎必然抓住符号错误、系数错误和漏项。把它设为MMS运行前的强制关卡。

症状1：核对不一致——按项分解锁定元凶

整体核对不一致时，把算子拆项逐一核对。例如对流扩散反应方程：

$$ s = \underbrace{\frac{\partial u_m}{\partial t}}_{s_1} + \underbrace{\mathbf{a}\cdot\nabla u_m}_{s_2} - \underbrace{\nabla\cdot(k\nabla u_m)}_{s_3} + \underbrace{r\,u_m}_{s_4} $$

把 \( s_1 \)～\( s_4 \) 分别与数值微分对比，坏项立现。高频元凶如下：

症状2：推导没错，可生成代码慢·坏

NS级方程的生成式可达数千次运算，代码生成环节有自己的一类问题：

• 编译不过·超行长——撞上Fortran续行上限或编译器复杂度限制。用 sympy.cse() 分解为中间变量链输出
• 运行极慢——同一个三角函数被求值几百次；照样CSE，另把 pow(x, 2) 改写为乘法
• 微小不一致（相对1e-10量级）——simplify改变了运算顺序带来的舍入差；无实害，把合格线定在1e-8左右即可
• 因消去而大幅不一致——生成了"大项相减"的形式；用 simplify／factor 消解相消后重新输出

症状3：代码没错，可在求解器里不起作用

UDF·用户子程序的注入环节，有与推导无关的坑：

1. 评估坐标搞错——该给单元中心的位置传了节点坐标（或反之）。有限体积用格心、FEM用积分点才对
2. 体积规约——求解器要的是"单位体积源"还是"单元总量"？网格加密一倍源项效果变四倍，就是规约错了
3. 时间没传对——非定常源项拿到的是上一步的时间，表现为一阶时间误差
4. 作用域遗漏——多区模型有的区没挂UDF；误差的空间分布在区界处不连续

此处最强的诊断是可视化误差的空间分布 \( |u_h - u_m| \)：推导错误铺满全域；注入错误带着结构——条带、分区、边界。

有章法的推进——逐项启用

源项流水线的质量清单

• 通过了随机点数值微分核对（相对误差≤1e-8为参考线）
• 做了"求解器方程形式（符号·系数·坐标系·守恒/非守恒）vs 推导式"的对应表
• CSE之后重新跑过核对（验证生成环节）
• 确认了源项评估位置（格心/积分点）与体积规约
• 可视化了误差空间分布（有结构⇒注入错；铺满全域⇒推导错）
• 记录了逐项启用的履历（到哪一级还出阶数）

相关文章：MMS源项的自动推导（整合版）、NS方程MMS的故障排查、收敛阶数的故障排查。