Grid Convergence Index: an error bar on a FEM result

Every finite element result contains discretisation error, and everyone knows it. What is usually missing is a number: how much error, in this result, on this mesh.

“I refined it and it did not change much” is the usual answer, and it is not one. The Grid Convergence Index turns the same three runs that produce that impression into a quantitative uncertainty band, using arithmetic that fits on a page.

Every number in this article is recomputed from its inputs by a script that runs on each build of this site.

Three meshesPermalink to “Three meshes”#

Solve the same problem on three systematically refined meshes — each element size a constant factor smaller than the last — and record the quantity of interest. Take element sizes of 1.00, 0.50 and 0.25 m, so the refinement ratio is $r = 2$, with results

$$f_3 = 15.620, \qquad f_2 = 15.324, \qquad f_1 = 15.240$$

numbered from coarse to fine. The differences are

$$\varepsilon_{32} = f_3 - f_2 = 0.296, \qquad \varepsilon_{21} = f_2 - f_1 = 0.084$$

The observed order of convergencePermalink to “The observed order of convergence”#

If the error behaves as $C h^p$, then successive differences fall by $r^p$, and the exponent can be recovered from the results themselves:

$$p = \frac{\ln\left|\varepsilon_{32}/\varepsilon_{21}\right|}{\ln r} = \frac{\ln 3.5238}{\ln 2} = 1.8171$$

This is the single most informative number in the exercise, and it is worth reading before anything else. It is the observed order, not the theoretical one, and comparing the two is a check on the whole model:

• Close to the theoretical order for the element — 2 for a quadratic displacement field on a smooth solution — means the mesh is fine enough to be in the asymptotic range and the refinement is doing what it should.
• Much lower means something is limiting the convergence: a singularity at a re-entrant corner, a discontinuity in material properties, a contact condition, or a solution that is simply not smooth. The rate is then a property of the problem, not a defect.
• Much higher, or negative, means the results are not converging monotonically, and the extrapolation below does not apply. Three points that go up, down and up again carry no information about a limit.

Here 1.8171 against a theoretical 2 is a mild reduction, consistent with a solution that is nearly but not perfectly smooth.

The extrapolated value and the bandPermalink to “The extrapolated value and the band”#

Richardson extrapolation estimates the value at zero element size:

$$f_{\text{ext}} = f_1 + \frac{f_1 - f_2}{r^p - 1} = 15.240 - \frac{0.084}{2.5238} = 15.2067$$

and the Grid Convergence Index converts the remaining difference into a reported uncertainty, with a safety factor of 1.25 for a three-grid study:

$$\mathrm{GCI}_{21} = \frac{1.25\,\left|\varepsilon_{21}/f_1\right|}{r^p - 1} = 0.27\%, \qquad \mathrm{GCI}_{32} = 0.96\%$$

So the answer is reported as 15.240 with a numerical uncertainty of 0.27%, which is 0.042 in the units of the quantity. The fine mesh sits 0.22% from the extrapolated value, and the 1.25 factor is what turns that estimate into something closer to a band than to a best guess.

The check that makes it trustworthyPermalink to “The check that makes it trustworthy”#

The whole procedure assumes the meshes are in the asymptotic range, where the error really is dominated by the $h^p$ term. That assumption is testable from the same three numbers:

$$\frac{\mathrm{GCI}_{32}}{r^p\,\mathrm{GCI}_{21}} = 0.9945$$

A value near 1 means the two grid pairs tell a consistent story about the convergence rate, and the extrapolation is defensible. A value of 0.6 or 1.5 means they do not, and the honest report is that the meshes are too coarse for an extrapolation — which is useful information, obtained cheaply, and much better than a confident number with no basis.

Why it belongs in every benchmarkPermalink to “Why it belongs in every benchmark”#

A verification benchmark compares a computed value against a reference. Without an error estimate on the computed value, agreement to 1% and disagreement by 1% are indistinguishable claims: the reader cannot tell whether the code is right and the mesh coarse, or the mesh fine and the code wrong.

With one, the comparison becomes a statement with two numbers in it. That is what makes a benchmark reproduced to 0.03% inside an uncertainty band of 0.21% a stronger claim than “close agreement” — the deviation is inside the numerical noise, and that is a conclusion rather than an impression.

Three things are worth remembering when doing it:

• Refine systematically. The refinement ratio must be uniform across the whole mesh, or $h$ has no single meaning. Refining locally, where the answer is, invalidates the exercise even though it improves the answer.
• Refine enough. A ratio of 1.1 leaves the differences buried in solver tolerance and floating-point noise. Aim for 1.3 at minimum, and 2 is comfortable.
• Choose the quantity first. The order of convergence differs between displacements, stresses and reactions on the same run, because they are different functionals of the same field. The GCI applies to the quantity you report, not to the analysis.

And the last point is the one worth keeping. This estimates discretisation error only. It says nothing about whether the constitutive model, the initial stress state or the geometry were right — it answers “are the equations solved accurately”, not “were they the right equations”.