Solution Verification and Convergence

Solution verification estimates numerical error in one computed result after the implementation has been tested. A familiar convergence slope is not enough: iterative error must already be smaller, at least three meaningfully separated resolutions must sample the same continuum problem, and the data must show an asymptotic regime rather than a fitted wish. Richardson analysis and a grid-convergence index are conditional tools for that regime.

Required background. Probability and Inference for Gravity supplies uncertainty summaries and interpretation of finite samples.

Helpful background. Formulations of General Relativity supplies the continuum variables and gauges being compared; Differential Equations for Relativity supplies consistency and convergence concepts; Computational Relativity supplies the shared numerical workflow; and Mathematical Relativity distinguishes numerical evidence from continuum convergence theorems.

Error components before extrapolation

Separate algebraic or iterative error, discretization error, roundoff, input interpolation, and model error. Make the estimated algebraic error in the solution or reported observable much smaller than the change between the two finest discretizations. A residual measures how nearly the discrete equations are satisfied; converting it into solution error requires their conditioning. Residual and solution norms can also have different units, so comparing their raw magnitudes is not a stopping rule.

For an invertible linear system $AU_*=b$ and approximation $\widehat U$, define $r=b-A\widehat U$. Then

$$U_*-\widehat U=A^{-1}r, \qquad \|U_*-\widehat U\|\le\|A^{-1}\|\,\|r\|,$$

using compatible vector and induced matrix norms. In a dimensionless example, $A=\operatorname{diag}(\varepsilon,1)$ and $r=(\varepsilon,0)^{\mathsf T}$ give unit solution error however small positive $\varepsilon$ becomes. For a nonlinear solve, the inverse Jacobian gives a local sensitivity estimate when the iterate is sufficiently close to a regular solution. Use such an estimate, or demonstrate that tighter solver tolerances change the observable negligibly compared with grid refinement; residual control alone does not prove this. The distinction between residual, backward error and forward error is developed in Higham 2002, § 7.1, Theorems 7.1–7.2.

Define one scalar quantity of interest $\phi$ or a norm on fields at a common physical time and frame. If waveforms are aligned, fix the alignment rule before measuring convergence and also show an unaligned diagnostic. If adaptive meshes are compared, define an effective resolution only when the refinement hierarchy changes systematically; counting total points is not a mesh spacing.

The standardized numerical-relativity tests in Babiuc et al. 2008, §§3–8 illustrate why constraints, solution error, and runtime behavior can converge differently. Report the quantity that supports the scientific claim, not only the norm that looks cleanest.

Observed order from three grids

Let grid 1 be finest and grid 3 coarsest, with spacings $h_1<h_2<h_3$. Prefer approximately uniform refinement ratios

$$r_{21}=\frac{h_2}{h_1}, \qquad r_{32}=\frac{h_3}{h_2},$$

larger than about $1.3$ so differences are not dominated by noise. For equal ratio $r$ and monotone convergence, an observed order is

$$p=\frac{\ln\left| \dfrac{\phi_3-\phi_2}{\phi_2-\phi_1} \right|}{\ln r}.$$

The extrapolated limit based on the fine pair is

$$\phi_{\mathrm{ext}} =\frac{r^p\phi_1-\phi_2}{r^p-1}.$$

These formulas follow by assuming $\phi(h)=\phi_*+Ch^p+o(h^p)$ with one nonzero leading coefficient $C$ over the sampled range. Subtracting adjacent values removes the unknown limit $\phi_*$; taking their ratio removes $C$. For unequal ratios, solve

$$\frac{\phi_3-\phi_2}{\phi_2-\phi_1} =\frac{h_3^p-h_2^p}{h_2^p-h_1^p}$$

for a positive $p$ under the same single-leading-term assumption. Taking absolute values does not repair an alternating-sign sequence: such a sequence does not fit this model with one fixed $C$. Oscillatory differences, near-zero denominators, shocks, mesh-interface motion, or mixed error orders can make a single $p$ meaningless. Celik and collaborators state a disciplined reporting procedure in Celik et al. 2008, “Preliminaries,” recommended Steps 1–5, Eqs. (1)–(10).

Do not force $p$ to the formal scheme order by rescaling curves. A lower stable order may reveal boundaries, interpolation, limiters, or nonsmooth data. An implausibly high order often signals cancellation or an unresolved leading term. Add a finer run and predict its difference before seeing it.

Grid-convergence index and its limits

For a nonzero scalar, a fine-grid relative grid-convergence index can be written

$$\operatorname{GCI}_{21} =F_s\, \frac{|\phi_1-\phi_2|} {|\phi_1|(r_{21}^{p}-1)},$$

where the safety factor $F_s$ and procedure are declared. This is an engineering estimate conditioned on the assumed error expansion and sampled regime. It is not a confidence interval derived from random sampling, nor does it include model, gauge, boundary, or calibration uncertainty.

A complete three-grid calculation

Suppose a dimensionless observable takes these values:

Grid	Spacing $h$	Observable $\phi$
finest, 1	0.05	2.0075
middle, 2	0.10	2.0300
coarsest, 3	0.20	2.1200

The differences are $0.0225$ and $0.0900$, so $r=2$, $p=\log_2(4)=2$, and $\phi_{\mathrm{ext}}=(4\times2.0075-2.03)/3=2$. The fine-grid absolute error estimate is $0.0225/(4-1)=0.0075$. With a declared safety factor $F_s=1.25$, the relative GCI is $1.25(0.0075)/2.0075\simeq0.00467$, or $0.467\%$.

These synthetic data follow $\phi(h)=2+3h^2$ exactly. A fourth grid at $h=0.025$ is therefore predicted to give $2.001875$. In a real calculation that prediction is a test, not a value to impose on the code. A small extra term or roundoff can spoil the fit, and the safety factor does not turn an untested expansion into a rigorous bound.

For fields, compute $L_1$, $L_2$, and $L_\infty$ norms after conservative restriction or projection to a common representation. State the measure and normalization: for example, a uniform-grid root-mean-square error is $[N^{-1}\sum_j e_j^2]^{1/2}$, whereas the discrete integral norm is $[h\sum_j e_j^2]^{1/2}$. They differ by the square root of the domain length. Local defects can disappear in an average norm, while $L_\infty$ can be dominated by a moving feature. For a waveform, report amplitude and phase errors mode by mode and verify that time and frame choices do not set the slope.

The periodic scalar-wave experiment supplies an actual four-grid calculation with separate spatial and temporal errors. Its analytic solution lets you compare the extrapolation logic with the true error before applying it to a black-hole simulation whose continuum result is unknown.

The benchmark suite should include a manufactured scalar wave with known continuum error, plus robust-stability, linear-wave, and gauge-wave tests. A physical black-hole or matter observable supplies an independent check. Manufactured smooth data verify the formal regime; a shock or mesh-refinement boundary tests the lower-regularity production regime.

Rerun and rejection rules

The environment record should pin equations, input checksums, code and dependency revisions, compiler, precision, solver tolerances, every mesh, refinement criterion, boundary treatment, output time, norms, and analysis scripts. Save iterative residual histories and raw $\phi_i$, not only the fitted $p$ and GCI.

A clean rerun would solve the scalar-wave benchmark on at least three grids with ratios above $1.3$, repeat the fit after adding a fourth finer grid, and emit residuals, norms, observed order, extrapolated value, and GCI. It would then repeat the logic on the claimed production observable. No absent executable is claimed here.

Reject the error estimate if estimated algebraic solution or observable errors are comparable to grid differences, grid ratios are too small, the three solutions do not represent identical continuum data, observed order changes materially after adding a finer run, convergence oscillates without a suitable model, or the result depends on forced rescaling. The principal failure mode is fitting an asymptotic formula before the asymptotic regime has begun.

Check your understanding

At spacings $h$, $2h$, and $4h$, a quantity is $1.0004$, $1.0016$, and $1.0064$. The differences suggest second-order convergence. However, the nonlinear-solver uncertainty in each value is estimated as $0.002$. May Richardson extrapolation support a four-digit result?

Solution

No. The fine–medium difference is only $0.0012$, smaller than the iterative uncertainty, so the apparent ratio of differences need not represent discretization error. Tighten the solver until its residual contribution is well below the grid differences, recompute all grids, and add a finer run to test the predicted second-order behavior. Only then evaluate the extrapolation and GCI.

References

• Babiuc, Maria C., et al. 2008. “Implementation of Standard Testbeds for Numerical Relativity.” Classical and Quantum Gravity 25, 125012. arXiv.
• Celik, Ismail B., et al. 2008. “Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications.” Journal of Fluids Engineering 130, 078001. DOI.
• Higham, Nicholas J. Accuracy and Stability of Numerical Algorithms. 2nd ed. Philadelphia: Society for Industrial and Applied Mathematics, 2002. DOI.