MMS源项的自动推导（符号计算）

MMS与源项的理论基础

在代码验证中的定位

把控制方程写作 \( L(u) = 0 \)，任意选取的制造解 \( u_m \) 一般不满足它。于是把残差定义为源项、加进方程：

$$ s(\mathbf{x}, t) := L(u_m), \qquad L(u) = s $$

按构造，\( u_m \) 是改造后方程的精确解。把源项 \( s \) 和由制造解导出的边界·初始条件交给求解器，确认数值解与 \( u_m \) 的误差随网格细化按理论阶数下降即可。

用观测阶数判定合格与否

在细化比为2的网格序列上测误差范数 \( E_h \)，观测收敛阶数为

$$ p_{obs} = \frac{\ln(E_{2h}/E_h)}{\ln 2} $$

二阶精度的离散化应有 \( p_{obs} \to 2 \)；达不到就说明离散化、边界实现或源项处理某处有缺陷——这是一条明确的合格判据。MMS的威力就在这里：把"答案看起来还行"的主观确认，换成"阶数出得来/出不来"的客观判定。

制造解的选取准则

制造解完全不需要物理意义，但有一些能最大化验证能力的选取准则：

• 足够光滑——为观测到名义阶数，用 \( C^\infty \) 级函数（三角函数与指数的组合是惯例）
• 让方程的每一项都起作用——选任何偏导项都不为零的解；常数或线性场检不出扩散项的bug
• 每个方向都有非平凡变化——多维时各方向取不同波数，可捕捉方向弄反的bug
• 顾及物理约束——湍流量、密度等必须为正的变量要配正的制造解（避免正值性限制器介入）

用符号计算推导源项

为什么自动推导是必需的

一维热传导的源项还能手算。但把三维制造解代入可压缩Navier-Stokes方程，源项会膨胀到几百项的规模——手工求导不可能不出错。而源项本身要是错了，"代码的bug"与"源项的bug"就无法区分，验证不再成立。因此用计算机代数系统（CAS）自动推导，是MMS事实上的前提条件。SymPy、Mathematica、Maple的流程都一样：

1. 把制造解 \( u_m(\mathbf{x},t) \) 定义为符号表达式
2. 原样施加控制方程的微分算子（\( \partial_t \)、\( \nabla\cdot \)、\( \nabla^2 \)……）
3. 化简（simplify）后，用代码生成功能输出为求解器语言

推导实例——非定常热传导

对热传导方程 \( \partial T/\partial t - \alpha \nabla^2 T = s \)，取制造解 \( T_m = \sin(\pi x)\cos(\pi y)\,e^{-t} \)，代入即得闭式源项：

$$ s = \frac{\partial T_m}{\partial t} - \alpha \nabla^2 T_m = \left(2\pi^2 \alpha - 1\right)\sin(\pi x)\cos(\pi y)\,e^{-t} $$

边界条件同样由制造解导出：Dirichlet边界给 \( T_m \) 的边界值，Neumann边界给 \( -k\,\partial T_m/\partial n \)。忘记让边界条件与制造解保持一致，是MMS实现中最高频的错误。

代码生成的实务——消去误差与公共子表达式

把符号表达式原样输出，运行时间和数值精度都可能出问题。三条常备对策：①公共子表达式消除（CSE），把重复的三角函数求值合并；②形如"两个大项相减"的表达式先simplify整理再输出（避免灾难性消去）；③用两种以上独立方法（另一套CAS，或对 \( u_m \) 做高阶数值微分）交叉核对生成代码，对源项本身做验证。第③条是"对验证工具的验证"，最容易被省略；仅在若干随机点与数值微分对比，就能抓住绝大多数推导错误。

实务应用流程

SymPy推导与代码输出

免费即可闭环的标准配置是SymPy。从推导到C/Fortran输出的骨架：

import sympy as sp

x, y, t, alpha = sp.symbols("x y t alpha")
T_m = sp.sin(sp.pi * x) * sp.cos(sp.pi * y) * sp.exp(-t)   # 制造解

L = sp.diff(T_m, t) - alpha * (sp.diff(T_m, x, 2) + sp.diff(T_m, y, 2))
s = sp.simplify(L)                                          # 源项

print(sp.ccode(s))          # 输出C表达式（UDF/用户子程序用）
print(sp.fcode(s))          # 输出Fortran表达式
s_num = sp.lambdify((x, y, t, alpha), s)                    # Python核对用函数

NS级的长表达式先用 sp.cse(s) 分解公共子表达式再输出，生成代码会短几个量级、快几个量级。

验证算例的标准构成

1. 网格序列——至少3档，最好4档；均匀细化（比值2）让阶数评估最简单
2. 误差范数——以 \( L^2 \) 为主、\( L^\infty \) 并记（局部降阶先出现在 \( L^\infty \)）
3. 时间与空间分离——测空间阶数时把时间步固定得足够小（或用定常MMS），反之亦然
4. 收紧迭代收敛——迭代残差要比离散化误差低两个量级以上，否则量到的是迭代误差
5. 阶数图——横轴 \( \log h \)、纵轴 \( \log E \)，叠加理论斜率参考线报告

注意源项的数值积分

有限元中源项经单元数值积分进入荷载向量。制造解取了高波数，源项会变得陡峭，默认积分阶数的积分误差可能盖过离散化误差、污染观测阶数。对策：把积分阶数提高一两档并确认结果不变，或降低制造解的波数。有限体积法里"格心值×体积"的近似对剧烈变化的源项也有同样问题。

工具支持与源项注入方法

推导侧的工具

求解器侧的注入手段

前沿研究动态

向复杂物理的扩展

MMS研究的前线在于"多复杂的物理还能把阶数验证做成"。湍流模型（Spalart-Allmaras、k-ω等）输运方程的制造解，难点是正值性与生成/耗散项的平衡设计——文献和MASA中已积累了专用解集。多相流界面追踪、刚性化学反应源项、动边界/ALE格式等含非光滑或强非线性的体系，仍是活跃的MMS研究对象。

非光滑问题与降阶的理论

含激波的可压缩流动中，即使高阶格式，全局阶数也会理论性地降到一阶。此时MMS进化为"光滑区域阶数是否出得来""间断捕捉宽度是否符合设计"的分区验证，并研究跨间断的误差范数处理（加权范数、区域分割评估）。不知道"阶数出不来才是正常"的情形，会把健康的代码误判为不合格。

接入CI（持续集成）

当代数值代码开发把MMS阶数测试作为自动回归测试嵌入CI已成标准做法：每次提交用两档粗网格计算观测阶数，超出容许带（如理论阶±0.2）就让构建失败。推导到代码生成都脚本化后，跟随方程或离散化变更的成本很小，可机械化地防住"不知何时阶数掉了"这类退化。

故障排查

观测阶数不达标时的诊断表

降阶的罪魁多半在边界

MMS不能证明"代码正确"，但它把"以此阶数收敛"这一可证伪性质变成了可系统检验的对象——这是代码验证的最有力手段。推导自动化之后维护成本很低，还能沉淀为回归测试资产。相关文章：MMS概述、MMS收敛阶数的故障排查、网格收敛性验证。