Automatic Derivation of MMS Source Terms (Symbolic Computation)

Theoretical Foundations of MMS and Source Terms

Where It Sits in Code Verification

Writing the governing equation as \( L(u) = 0 \), an arbitrarily chosen manufactured solution \( u_m \) generally won't satisfy it. So define the residual as a source term and add it to the equation:

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

By construction, \( u_m \) is the exact solution of the modified equation. Feed the solver the source term \( s \) and the boundary/initial conditions derived from \( u_m \), then confirm that the error between the numerical solution and \( u_m \) decreases at the theoretical order under mesh refinement.

Pass/Fail by Observed Order of Accuracy

Measuring an error norm \( E_h \) on a mesh sequence with refinement ratio 2, the observed order of convergence is

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

For a second-order discretization, \( p_{obs} \to 2 \) — and if it doesn't get there, something in the discretization, boundary implementation, or source handling is defective. That is MMS's power: it replaces the subjective "the answer looks plausible" with the objective "the order comes out or it doesn't."

Criteria for Choosing Manufactured Solutions

The manufactured solution needs no physical meaning, but some choices maximize its verification power:

• Smooth enough — \( C^\infty \) functions (combinations of trig and exponentials are standard) so the formal order is observable
• Exercises every term — pick solutions for which no derivative term vanishes; constants or linear fields cannot detect diffusion-term bugs
• Nontrivial variation in every direction — different wavenumbers per direction in multi-D catch direction-swapping bugs
• Respect physical constraints — quantities that must be positive (turbulence variables, density) should get positive manufactured solutions so positivity limiters stay inactive

Deriving Source Terms by Symbolic Computation

Why Automated Derivation Is Mandatory

For 1-D heat conduction you can derive the source term by hand. Substitute a 3-D manufactured solution into the compressible Navier-Stokes equations, however, and the source term runs to hundreds of terms — manual differentiation cannot be done without mistakes. And if the source term itself is wrong, "bug in the code" and "bug in the source term" become indistinguishable, so verification collapses. Symbolic derivation with a CAS is the de facto prerequisite for MMS. Whether SymPy, Mathematica, or Maple, the procedure is the same:

1. Define the manufactured solution \( u_m(\mathbf{x},t) \) symbolically
2. Apply the governing operator verbatim (\( \partial_t \), \( \nabla\cdot \), \( \nabla^2 \), …)
3. Simplify, then emit solver-language code via code generation

A Worked Example — Transient Heat Conduction

For the heat equation \( \partial T/\partial t - \alpha \nabla^2 T = s \) with the manufactured solution \( T_m = \sin(\pi x)\cos(\pi y)\,e^{-t} \), substitution gives the closed form directly:

$$ 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} $$

Boundary conditions come from the same solution: Dirichlet boundaries get the boundary values of \( T_m \), Neumann boundaries get \( -k\,\partial T_m/\partial n \). Forgetting to make the boundary conditions consistent with the manufactured solution is the single most common MMS implementation error.

Code Generation in Practice — Cancellation and CSE

Emitting raw symbolic expressions can hurt both runtime and accuracy. Three standing rules: ① apply common subexpression elimination (CSE) so trig calls aren't evaluated hundreds of times; ② simplify forms that subtract large nearly-equal terms before emitting (avoids catastrophic cancellation); ③ verify the source term itself by comparing the generated code against an independent evaluation — another CAS, or high-order numerical differentiation of \( u_m \) at random points. Step ③ is "verifying the verifier" and is often skipped; a random-point finite-difference cross-check catches the vast majority of derivation mistakes.

Applying It in Practice

SymPy Derivation and Code Output

A free, self-contained toolchain: from derivation to C/Fortran output.

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)   # manufactured solution

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

print(sp.ccode(s))          # C expression (for UDFs / user subroutines)
print(sp.fcode(s))          # Fortran expression
s_num = sp.lambdify((x, y, t, alpha), s)                    # Python checker

For Navier-Stokes-scale expressions, decompose with sp.cse(s) before emitting — the generated code becomes dramatically shorter and faster.

Standard Structure of a Verification Campaign

1. Mesh sequence — at least 3 levels, preferably 4; uniform refinement (ratio 2) keeps order evaluation simple
2. Error norms — \( L^2 \) as primary, \( L^\infty \) alongside (local order loss shows up in \( L^\infty \) first)
3. Separate space and time — freeze the time step (or use steady MMS) when measuring spatial order, and vice versa
4. Tighten iterative convergence — drop residuals at least two orders below the discretization error, or you're measuring iteration error
5. Order plot — \( \log h \) vs. \( \log E \) with the theoretical-slope reference line

Watch the Numerical Integration of the Source Term

In FEM the source term enters the load vector through element-level quadrature. A high-wavenumber manufactured solution makes the source steep, and default quadrature orders can introduce integration errors that dominate discretization error and pollute the observed order. Remedies: raise the quadrature order one or two levels and confirm results don't change, or lower the solution's wavenumber. Finite-volume codes hit the same issue via the cell-centroid-times-volume approximation for strongly varying sources.

Tooling and Source-Term Injection

Derivation-Side Tools

Injection Paths into Solvers

Research Frontiers

Extending MMS to Complex Physics

The frontier is making order verification work for ever more complex physics. Manufactured solutions for turbulence-model transport equations (Spalart-Allmaras, k-ω families) require careful positivity and production/dissipation balance — dedicated solution sets have accumulated in the literature and in MASA. Interface tracking in multiphase flow, stiff reacting source terms, and moving-boundary/ALE formulations remain active MMS research targets.

Non-Smooth Problems and the Theory of Order Loss

With shocks, even high-order schemes drop to first-order global convergence — that is theory, not a bug. Modern practice verifies "does the order emerge in smooth regions" and "is the captured discontinuity width as designed" separately, using weighted or region-split norms. Knowing the cases where the order is supposed not to appear prevents failing a healthy code.

Wiring MMS into CI

Contemporary numerical-code development embeds MMS order tests as automated regression tests: every commit runs two coarse mesh levels, computes the observed order, and fails the build outside a tolerance band (e.g., theoretical order ±0.2). With derivation-to-codegen scripted, the tests track equation and discretization changes cheaply — mechanically preventing the "order silently degraded at some point" class of regressions.

Troubleshooting

Diagnosis Table When the Observed Order Fails

The Culprit Is Usually the Boundary

MMS never proves a code "correct," but it makes the falsifiable property "converges at this order" systematically testable — the strongest tool in code verification. With derivation automated, its maintenance cost is small and it compounds into a regression-test asset. Related: MMS overview, troubleshooting observed convergence rates, mesh convergence verification.