Method of Manufactured Solutions

The method of manufactured solutions chooses a smooth field first, substitutes it into the continuum equations, and treats the resulting nonzero forcing as part of a new exact problem. Because the manufactured field can activate every derivative, coupling, coordinate term, boundary closure, and mesh interface, it can expose implementation defects that a symmetric physical solution never excites.

Required background. Probability and Inference for Gravity supplies residual and convergence summaries.

Helpful background. Formulations of General Relativity supplies the coupled continuum systems to be exercised; Differential Equations for Relativity supplies consistency and boundary-value concepts; Computational Relativity supplies the volume-wide workflow; and Mathematical Relativity supplies regularity and domain hypotheses.

Manufacture the continuum problem

Write the target equations schematically as

$$\mathcal L[u]=0.$$

Choose a sufficiently smooth, nontrivial field $u_{\mathrm M}$ on the full spacetime domain and define

$$f_{\mathrm M}=\mathcal L[u_{\mathrm M}].$$

The verified code then solves $\mathcal L[u]=f_{\mathrm M}$ with initial and boundary data obtained from $u_{\mathrm M}$. The numerical solution should approach $u_{\mathrm M}$ at the designed order. This verifies implementation of the forced equations; it does not validate the unforced physical model.

For a simple $1+1$ scalar wave,

$$\partial_t^2u-\partial_x^2u=f_{\mathrm M},$$

on the periodic interval $0\le x\lt1$ for $0\le t\le1$, choose

$$u_{\mathrm M}(t,x) =\sin(2\pi x)\cos(3t) +0.1\cos(4\pi x)\sin(2t).$$

Direct substitution gives

$$\begin{aligned} f_{\mathrm M}={}& \left(4\pi^2-9\right)\sin(2\pi x)\cos(3t)\\ &+0.1\left(16\pi^2-4\right)\cos(4\pi x)\sin(2t). \end{aligned}$$

The initial data are $u(0,x)=\sin(2\pi x)$ and $\partial_tu(0,x)=0.2\cos(4\pi x)$, with periodic matching for the field and derivatives. This example has nonzero time and space derivatives, but it does not test curved geometry, nonlinear coupling, or physical boundary closures. A production verification family must be designed around the actual operators.

A variable coefficient that exposes a missing derivative

Consider instead the dimensionless problem

$$u_{tt}-\partial_x\!\left(a(x)u_x\right)=f, \qquad a(x)=1+x, \qquad 0\le x\le1.$$

Choose $u_{\mathrm M}=e^{-t}(1+x+x^2)$. The product rule gives

$$\begin{aligned} \partial_x(a u_x) &=a u_{xx}+a' u_x=e^{-t}(3+4x),\\ f_{\mathrm M} &=e^{-t}(x^2-3x-2). \end{aligned}$$

Use $u(0,x)=1+x+x^2$, $u_t(0,x)=-(1+x+x^2)$, the left Dirichlet condition $u(t,0)=e^{-t}$, and the right Robin condition $u_x(t,1)+u(t,1)=6e^{-t}$. These data match the manufactured field at the initial corners. The coefficient $a$ is positive, so the continuum wave speed is $\sqrt{a(x)}$ and the time-step restriction must account for its maximum $\sqrt2$.

A code that replaces $\partial_x(a u_x)$ by $a u_{xx}$ misses $a'u_x=e^{-t}(1+2x)$. This term does not vanish under refinement. Evaluating the incorrect continuum operator on the exact field leaves precisely that residual, making the fixture sensitive to a defect that a constant-coefficient test cannot detect. This is still a scalar test; it makes no claim to exercise Einstein constraints or curved-spacetime tensor couplings.

Excite every implemented path

Use fields whose components and derivatives do not vanish accidentally. In a tensor system, assign distinct analytic functions to independent components so an index swap cannot survive by symmetry. In curvilinear coordinates, avoid choices that cancel connection terms. For coupled matter–geometry equations, make every source and flux component nonzero while preserving required algebraic constraints.

Boundary conditions are part of the implementation. Derive Dirichlet, Neumann, Robin, constraint-preserving, or characteristic boundary data from the manufactured field and forcing. Test corners and patch boundaries. For adaptive mesh refinement, place smooth structure across refinement interfaces and move it through regridding so prolongation, restriction, synchronization, and reflux operations are exercised.

The manufactured solution must possess enough derivatives for the formal order being tested. A discontinuous or barely differentiable fixture cannot demonstrate high-order convergence. Conversely, a single low-frequency trigonometric mode may be too easy and can hide dispersion or aliasing. Sweep wavelength relative to grid spacing without leaving the resolved regime.

Salari and Knupp give practical construction and boundary guidance in Salari and Knupp 2000, Chapter 3, §§3.3–3.4, Chapter 4, and Chapter 5, §5.2. Roache emphasizes the logical role of manufactured solutions in Roache 2002, pp. 4–10.

Keep the oracle independent

Generating $f_{\mathrm M}$ with the same symbolic expressions or helper functions used in the solver can reproduce the same sign or index error. Derive the forcing in an independent algebra system or hand calculation, export a simplified expression, and test sample values against a second evaluator. Keep conventions—metric signature, curvature sign, index positions, and units—explicit at the interface.

Compare numerical and exact fields in $L_1$, $L_2$, and $L_\infty$ norms on at least three resolutions. Monitor constraints and boundary residuals separately. The global $L_2$ rate can look correct while a lower-order boundary defect remains localized; $L_\infty$ and region-restricted norms reveal it.

Mutation is a decisive check. Reverse one forcing sign, interchange two indices, omit one source term, change a stencil coefficient, or corrupt an AMR prolongation weight in an isolated implementation. The manufactured test should fail with the expected norm pattern. Restore the term and recover convergence. A test that survives a relevant mutation has not verified that path.

Benchmark and reproduction steps

The scalar-wave fixture above supplies a transparent benchmark, augmented by the Apples-with-Apples linear and gauge waves for a relativistic system and an independent physical benchmark for the final observable. Run at no fewer than three grids, make the estimated iterative solution error smaller than discretization differences, and add a finer grid to confirm the observed order. Establish that error through conditioning estimates or tolerance refinement, as explained in solution verification; a raw residual is not a solution-error bound. Refine boundary and AMR components with the interior rather than holding one error source fixed.

The environment record should pin the continuum equations, manufactured fields and forcing generator, independent evaluator, code and dependency revisions, compiler, precision, meshes, time step, boundary data, AMR schedule, tolerances, norms, and analysis scripts. Outputs include exact and numerical fields, forcing samples, residual histories, regional norms, orders, and mutation outcomes.

A clean rerun would regenerate the forcing from versioned symbolic input, compare it at fixed sample points with the independent evaluator, solve the benchmark over the refinement sequence, and repeat after each deliberate mutation. This page specifies the sequence without claiming an absent executable ran it.

Reject the code-verification claim if important terms vanish for every fixture, forcing and solver share an unchecked implementation, boundary or AMR regions are excluded, convergence misses the designed order without explanation, or a consequential mutation survives. The central limitation is scope: manufactured smooth solutions verify exercised code paths, not the physical validity or shock behavior of a different regime.

Check your understanding

A fourth-order evolution code achieves fourth-order $L_2$ convergence for a periodic manufactured wave, but its production problem uses mesh refinement and an outer Robin boundary. Has the production implementation been verified?

Solution

Only the periodic interior path has been verified. Construct a smooth manufactured field whose derived Robin data are nonzero and whose structure crosses refinement interfaces and regridding events. Measure boundary, interface, and global norms separately, then mutate the boundary sign and one prolongation coefficient. The production claim remains incomplete until those paths fail under mutation and recover the expected convergence after restoration.

References

• Roache, Patrick J. 2002. “Code Verification by the Method of Manufactured Solutions.” Journal of Fluids Engineering 124, 4–10. DOI.
• Salari, Kambiz, and Patrick Knupp. 2000. Code Verification by the Method of Manufactured Solutions. Sandia Report SAND2000-1444. DOI.