MMS for 2-D Steady Heat Conduction — Troubleshooting Guide

2-D Steady Heat Conduction — The Right "First Problem" for MMS

A steady heat-conduction MMS is the ideal entry point to code verification, because the source term can be worked out by hand. Take the unit square with a constant conductivity \( k \) and pick the manufactured solution \( T_m = \sin(\pi x)\sin(\pi y) \), and you get

$$ s = -k\,\nabla^2 T_m = 2\pi^2 k \,\sin(\pi x)\sin(\pi y) $$

in a single line. The boundary can be \( T_m = 0 \) (Dirichlet) all the way around, so no symbolic computation is needed at all. The shortest learning path is to absorb the shape of a verification study — mesh sequence, error norm, observed order — on this one fully checkable problem, and only then move on to elasticity and fluids. Even so, this simple problem has several places to get stuck; this page clears them symptom by symptom.

The Five-Step Standard Procedure

1. Define the source term \( s(x,y) \) in the solver as a volumetric heat source (units W/m³, consistent with the units of \( k \))
2. Apply Dirichlet boundaries \( T = 0 \) all around (in general, the boundary values of \( T_m \))
3. Solve on three mesh levels (for example 10×10, 20×20, 40×40)
4. Compute the error norm \( E_h = \|T_h - T_m\|_{L^2} \) at each level
5. Compare the observed order \( p = \ln(E_{2h}/E_h)/\ln 2 \) against the element's theoretical order (2 for linear elements, 3 for quadratic — in the \( L^2 \) norm)

Symptom 1: The Error Never Gets Small — Obviously a Different Solution

Symptom 2: It Converges, but the Order Falls Short of Theory

• How the spatial distribution of the source is supplied — if the tool only accepts a single constant volumetric heat source over a region, the source becomes an element-wise constant approximation and the accuracy is capped at first order. Supply it as a coordinate-dependent expression, a user function, or a sufficiently fine interpolated table
• How the error is measured — is the exact solution evaluated at the same location as the numerical one, and does the \( L^2 \) norm carry volume weights? (Details in the convergence-order troubleshooting guide)
• When Neumann edges are mixed in — the sign of the heat flux \( q = -k\,\partial T_m/\partial n \) (outward positive or inflow positive: follow the tool's convention), and the quadrature accuracy for values that vary along the edge
• The error goes to "zero" with quadratic elements — a polynomial manufactured solution of degree two or lower is reproduced exactly by quadratic elements, so no order can be measured. Use a trigonometric form (the \( T_m \) on this page is safe)

Symptom 3: It Stopped Matching Once k(T) Became Temperature-Dependent

The natural next step is to introduce \( k(T) \), and the easy thing to overlook is that the source-term derivation changes. Correctly,

$$ s = -\nabla\cdot\left(k(T_m)\,\nabla T_m\right) = -k(T_m)\nabla^2 T_m - k'(T_m)\,|\nabla T_m|^2 $$

and forgetting the second term (the \( k' \) term) is the standard mistake. From here on it is safer to switch to symbolic derivation (SymPy or similar). The problem is also nonlinear now, so drive the iteration (Picard/Newton) two decades below the discretization error — a shallowly converged nonlinear run will happily destroy the order you were getting in the linear case.

Symptom 4: The Error Stops Falling on Fine Meshes

This simple problem is usually where you first meet the "error floor". The causes are ① the convergence criterion of the linear solver (for iterative methods), ② round-off (once the error reaches the \( 10^{-10} \) range), and ③ approximations mixed into the evaluation of the source term or boundary values (a coarse interpolation table, for instance). Diagnosis takes two moves: look at the absolute value of the error, and tighten the convergence criterion by one decade to see what changes. Going all the way to "refine until you hit the floor, identify what the floor is made of, and report it" is what raises a verification report from adequate to complete.

Symptom 5: The Commercial Tool Won't Accept a Coordinate-Dependent Source

Whether MMS is feasible depends less on the solver itself than on "can you define a coordinate-dependent volumetric heat source?" The routes in the major tools: Ansys Mechanical = function definition (Function, or %table% and coordinate functions in APDL); Abaqus = an Analytical Field or a DFLUX subroutine; COMSOL = type the expression directly (the easiest); OpenFOAM = codedSource/fvOptions; Fluent = a UDF (DEFINE_SOURCE). Even in tools where no expression can be entered from the GUI, a fine spatial table plus interpolation serves as a substitute (make it fine enough that the interpolation error does not pollute the order — see Symptom 2). Establish this "MMS route in my own tool" once, and every verification study afterwards reuses it.

Checklist and Next Steps

• Checked the sign convention, coefficient, and units of the source term against the solver's own way of writing the equation
• Verified every boundary assignment visually (remember: unassigned means insulated)
• Confirmed how the error norm is measured (matching locations, volume weights)
• Computed the observed order over three or more levels and compared it against the theoretical order
• Identified what the error floor is made of (convergence criterion / round-off)
• Next steps: extended in stages to a version with Neumann/Robin boundaries mixed in, a nonlinear k(T) version, and an unsteady version (with the temporal order measured separately)

Related: MMS overview (consolidated), convergence-order troubleshooting, 2-D elasticity MMS (the next step).