Grid Convergence Index (GCI) Calculator
Grid Convergence Index (GCI) and Richardson Extrapolation Calculator
Three grids, three solutions, and the number a reviewer asks for: the observed order of convergence, the Richardson-extrapolated value and the fine-grid GCI, with the transcendental order equation solved properly rather than approximated, oscillatory and diverging cases detected and refused, and the ratio usually sold as an asymptotic-range check shown for what it is.
Grid convergence index and observed order
A drag coefficient from three 3-D grids of 8.0, 2.7 and 1.0 million cells giving 0.29850, 0.30150 and 0.30673, a formally second-order solver, and the three-grid safety factor of 1.25
Richardson extrapolation, the apparent order, and the Grid Convergence Index
- h
- representative cell size. ASME V&V 20 and Celik et al. define it as the cube root of the mean cell volume in 3-D, or the square root of the mean cell area in 2-D. The domain volume is common to all three grids and cancels from every ratio, so cell counts alone are sufficient
- r
- refinement factor between two grids. Celik et al. recommend it be at least 1.3, which in three dimensions means at least 2.2 times the cell count
- p
- observed, or apparent, order of convergence — the order the solver actually achieved on these grids and this quantity, not the order the scheme was designed for. When the two refinement ratios are equal, q(p) vanishes and p = ln|ε₃₂/ε₂₁| / ln r. When they are not, the equation is transcendental in p
- s
- the sign of ε₃₂/ε₂₁. It is +1 for monotone convergence. A negative value is Celik et al.’s indicator of oscillatory convergence, and in that case the observed order is not defined
- φ_ext
- the Richardson extrapolation of the three solutions to zero cell size. It is an estimate of the exact solution of the discretised equations, not of reality — it says nothing about the turbulence model or the geometry
- e_a
- approximate relative error: the between-grid difference as a fraction of the fine-grid solution. Useful, and routinely mistaken for an error estimate. It is not one — it does not know the order and it has no safety factor
- F_s
- safety factor. Roache recommends 1.25 over three or more grids where the order has been observed, and 3 for two grids with the order assumed
- GCI
- Grid Convergence Index: the error estimate rescaled so that a GCI of 2% on a first-order scheme means the same thing as 2% on a second-order one, which is the whole reason it exists and the reason it can be compared between studies
Worked example
A drag coefficient from three 3-D grids of 8.0, 2.7 and 1.0 million cells giving 0.29850, 0.30150 and 0.30673, a formally second-order solver, and the three-grid safety factor of 1.25
Representative cell sizes first. The domain volume is the same on all three grids, so h is proportional to N to the power minus one third and the volume cancels: r₂₁ = (8.0/2.7) to the power 1/3 = 1.4363 and r₃₂ = (2.7/1.0) to the power 1/3 = 1.3925. Both are above the 1.3 minimum Celik et al. recommend, and — importantly for what follows — they are not equal
The differences: ε₂₁ = φ₂ − φ₁ = +0.00300 and ε₃₂ = φ₃ − φ₂ = +0.00523. Same sign, so convergence is monotone, s = +1, and the ratio ε₃₂/ε₂₁ = 1.7433 is greater than one, so refining is reducing the change rather than increasing it
NOW THE OBSERVED ORDER, AND THIS IS WHERE MOST GCI CALCULATIONS GO WRONG. The closed form p = ln|ε₃₂/ε₂₁| / ln r₂₁ = ln(1.7433)/ln(1.4363) = 1.535 is only valid when the two refinement ratios are equal. They are not. The correct equation carries a term q(p) = ln[(r₂₁^p − s)/(r₃₂^p − s)] which itself contains p, so p appears on both sides and the equation is transcendental. Solving it properly gives p = 1.8708 — the closed form is 18% low here, and it will be low in the same direction whenever r₂₁ > r₃₂
Richardson extrapolation: φ_ext = (r₂₁^p·φ₁ − φ₂)/(r₂₁^p − 1) with r₂₁^p = 1.9686 gives 0.295403. The extrapolation moves the fine-grid answer by −0.003097, which is slightly more than the whole fine-to-medium difference — that is normal at an order below two, and it is the reason a grid study is worth doing even when the two finest grids look close
The error measures. Approximate relative error e_a = |ε₂₁/φ₁| = 1.005%: that is the number people quote as "grid independence to 1%", and it is not an error estimate — it has no order in it and no safety factor. Extrapolated relative error e_ext = |(φ_ext − φ₁)/φ_ext| = 1.048%
THE GRID CONVERGENCE INDEX: GCI₂₁ = Fs·e_a/(r₂₁^p − 1) = 1.25 × 0.010050/0.96862 = 1.297%. One grid coarser, GCI₃₂ = 2.528%. Report the fine-grid figure with the result: Cd = 0.2985 ± 1.3%
What the two-grid convention would have given. With only grids 1 and 2, the formal order of 2 assumed and Roache's Fs = 3, GCI = 3 × 0.01005/(1.4363² − 1) = 2.837%. Two and a bit times larger for the same two solutions — which is the entire value of the third grid, and the reason quoting 1.25 on a two-grid study understates the uncertainty by roughly the same factor
AND THE CHECK EVERYONE REPORTS IS NOT A CHECK. The ratio r₂₁^p·GCI₂₁/GCI₃₂ comes out as 1.010050, reassuringly close to 1. It is close to 1 for any three numbers you like: substitute the definitions and it reduces exactly to |φ₂/φ₁|, which here is 0.30150/0.29850 = 1.010050 to every digit printed. The same is true of φ_ext computed from grids 2 and 3 — it is not a second, independent extrapolation, it is algebraically the same number as the first, and the calculator prints their difference to show it is zero. The real asymptotic-range test is the observed order against the formal order: 1.871 against 2, a ratio of 0.94, which is a genuinely healthy result
Finally, what 1% would cost. GCI falls as h to the power p, and h goes as N to the power minus 1/3, so the cell count needed scales as (GCI ratio) to the power d/p = 3/1.871 = 1.60. Getting from 1.297% to 1.0% needs 1.52 times the cells, about 12.1 million. Getting to 0.5% would need 4.6 times, about 37 million
What the observed order is telling you
| Observed p against formal p | What it usually means | What to do |
|---|---|---|
| p close to the formal order | you are in or near the asymptotic range and the extrapolation is trustworthy | report GCI with Fs = 1.25 and move on |
| p between about 0.6 and 0.9 of formal | limiters, wall functions, skewed cells or a geometric singularity are eroding the order; extremely common on industrial meshes | accept it, use the observed p in the GCI, and expect cells to buy less than the formal order suggests |
| p well below 1 | the quantity may depend on something that is not converging with h at all, or the coarse grid is off the asymptotic curve | check the functional, then add a fourth grid and test the order on the finer triple |
| p above the formal order | not possible in the asymptotic range, so you are not in it. Usually two error terms of opposite sign cancelling on one grid | add a finer grid; do not report the flatteringly small GCI this produces |
| p noisy when a grid is swapped | the between-grid differences are close to the iterative-convergence error | tighten the residual tolerance by two orders and repeat before touching the mesh |
What a factor of refinement costs, and what it buys
| Refinement factor r | Cell-count multiplier in 3-D | GCI falls by, at p = 1 | at p = 2 | at p = 3 |
|---|---|---|---|---|
| 1.15 | 1.52× | 1.15× | 1.32× | 1.52× |
| 1.30 | 2.20× | 1.30× | 1.69× | 2.20× |
| 1.50 | 3.38× | 1.50× | 2.25× | 3.38× |
| 2.00 | 8.00× | 2.00× | 4.00× | 8.00× |
| 3.00 | 27.0× | 3.00× | 9.00× | 27.0× |
Four things a GCI calculation gets wrong, three of which make the answer look better than it is
The Grid Convergence Index turns “I refined the mesh and the answer barely moved” into a number with a defensible meaning, and it is the number a reviewer asks for. It is Roache’s proposal, set out in his 1994 paper and carried into the ASME V&V 20 framework and the procedure Celik and co-authors published as an editorial policy statement for the ASME Journal of Fluids Engineering. All of it is arithmetic on three solutions and three grid sizes. The reason it is worth a page is that four separate things about it are routinely got wrong, and three of them make the reported uncertainty too small.
First: the observed order, not the formal order. The order that goes into the GCI is the order the solver actually achieved on your grids and your quantity — p extracted from the three solutions — not the order the scheme was designed for. On a clean structured mesh with a smooth solution a second-order scheme delivers close to 2. On a mixed polyhedral mesh with limiters, wall functions and a separation point it commonly delivers between 1 and 1.6, because every one of those features is locally lower order. Using 2 when the solver delivered 1.2 understates the GCI by the ratio (r² − 1)/(r^1.2 − 1), which for a refinement factor of 1.3 is a factor of 2.7.
Second: the equation for the observed order is transcendental unless the refinement ratios are equal, and the closed form is not a reasonable approximation to it. With r₂₁ = r₃₂ the whole q(p) term vanishes and p = ln|ε₃₂/ε₂₁| / ln r. With unequal ratios — which is the normal case, because nobody gets exactly 2.2 times the cells twice in a row from a mesher — p sits on both sides of the equation. In the worked example above, ratios of 1.44 and 1.39 differ by only 3%, and the closed form is still 18% low. It is worth arranging equal ratios deliberately if you still can, because then no iteration is needed at all.
How this page solves it, and why not the way the procedure says to. Celik et al. suggest fixed-point iteration, starting from the closed-form value. That iteration is not unconditionally convergent. Its asymptotic contraction factor is 1 − ln r₃₂/ln r₂₁, so it diverges whenever ln r₃₂/ln r₂₁ exceeds 2, that is whenever r₃₂ > r₂₁². That is not an exotic case: r₂₁ = 1.2 with r₃₂ = 2.5 is a perfectly ordinary pair of grid steps and the published iteration runs away from the answer. This page instead rewrites the equation as a single monotone increasing function of p and solves it by thirty bisections on the bracket 0.005 to 32, which cannot diverge because it cannot leave the bracket. The residual of the equation is printed under the answer so you can see it solved, and the calculator says so when it did not — which happens when no positive order satisfies the equation at all, the condition for which is that ε₃₂/ε₂₁ be below ln r₃₂/ln r₂₁.
Third: oscillatory convergence, which is common and is usually reported as if it had not happened. If φ₂ does not lie between φ₁ and φ₃ then the three solutions are not marching towards a limit, the sign s in the equation goes negative and the observed order stops meaning anything. The published equation will still return a number, because it takes an absolute value. So will divergence: if refining the mesh moves the answer more rather than less, ε₃₂/ε₂₁ falls below one, the logarithm goes negative, and the absolute value hands you back a perfectly respectable positive order with the sign quietly discarded. This page detects both and refuses to call either a GCI. For the oscillatory case it reports instead the one honest statement three oscillating values support: half their spread, as a percentage of the fine-grid value. The most frequent innocent explanation for oscillation, incidentally, has nothing to do with the mesh — it is a time-averaged quantity averaged over a different number of shedding cycles on each grid, which you can read about beside the Courant number and time step.
Fourth, and this one is a genuine error in the literature rather than a misreading of it: the ratio everybody reports as the asymptotic-range check is algebraically a tautology. The check is usually written GCI₃₂ ≈ r^p·GCI₂₁, and NASA Glenn’s much-copied tutorial prints the resulting ratio as a headline verification that the solutions are in the asymptotic range. Substitute the definitions. GCI₂₁ = Fs|ε₂₁/φ₁|/(r₂₁^p − 1) and GCI₃₂ = Fs|ε₃₂/φ₂|/(r₃₂^p − 1). The apparent-order equation states that |ε₃₂/ε₂₁| = r₂₁^p(r₃₂^p − 1)/(r₂₁^p − 1). Put them together and every factor cancels except one: r₂₁^p·GCI₂₁/GCI₃₂ = |φ₂/φ₁|, exactly, for any three solutions and any pair of refinement ratios, converging or not. It is close to 1 because φ₂ is close to φ₁, and for no other reason. The same holds for the second extrapolated value: φ_ext from grids 2 and 3 is not an independent estimate to compare against φ_ext from grids 1 and 2, it is identically the same number whenever convergence is monotone. This page prints both pairs side by side, and their differences, so that the point is visible rather than argued. The test that does carry information is the observed order against the formal order, and after that, a fourth grid: compute p from grids 1–2–3 and again from 2–3–4 and see whether they agree.
What a GCI is not. It is an estimate of the discretisation uncertainty in one scalar functional, on the assumption that the error behaves as a single power of h. It says nothing about the turbulence model, the boundary conditions, the fluid properties or the geometry, and a 0.3% GCI on a k-epsilon solution of a separated flow is a precise statement about a wrong answer. It is also specific to the functional: the same three solutions can give 0.5% on integrated lift and 15% on peak wall heat flux, and both are correct. Run the study on the number you are going to report. Finally, it assumes each solution is iteratively converged far tighter than the differences between grids — if your residuals have fallen three orders and the grids differ by 1%, you are measuring the solver, and no amount of GCI arithmetic will tell you so.
Frequently asked questions
Do I really need three grids?
For a GCI with the 1.25 safety factor, yes, because 1.25 is justified by having observed the order rather than assumed it. With two grids you can still report a GCI, using the formal order of your scheme and Roache’s safety factor of 3; this page prints that figure beside the three-grid one so you can see the cost of the missing grid. In the worked example it is 2.84% against 1.30%, a factor of 2.2. What you must not do is use 1.25 with an assumed order, which is the commonest way a published GCI comes out too small.
My observed order came out at 3.4 on a second-order solver. Is that good?
No, it is a warning. A scheme cannot exceed its own formal order in the asymptotic range, so an observed order above it means you are not in the asymptotic range. The usual mechanism is that two error terms of opposite sign happen to cancel on one of the three grids, which flatters the fit and makes the GCI far too small. Add a finer grid and recompute the order on the finest three. The same applies in the other direction to orders above about 1.15 times the formal one, which is where this page starts warning.
What refinement ratio should I use?
At least 1.3 in cell size, which Celik et al. recommend explicitly and describe as based on experience rather than theory. In three dimensions that means at least 2.2 times the cell count per step, and a full factor of two in cell size means eight times the cells. Going finer in steps than 1.3 makes the between-grid differences small enough that iterative-convergence error and round-off start to dominate them, and the observed order becomes wildly sensitive — you will see it swing by a whole unit when you change a residual tolerance. Equal ratios in both steps are worth arranging because they remove the transcendental equation entirely.
How do I refine a mesh by a factor of 1.3 when the mesher only takes a base size?
Scale every length control by the same factor and let the count fall where it falls, then compute h from the counts the mesher actually produced rather than from what you asked for. That is what the representative cell size h = (V/N)^(1/3) is for: it measures the mesh you got. The refinement must be uniform — refining only the boundary layer, or only a refinement region, breaks the single-power-of-h assumption the whole method rests on, and it is the reason a study can produce an observed order of 0.3 on an otherwise well-behaved case. If you keep the near-wall stack fixed while refining the volume, you are not running a grid convergence study; see the first cell height pages for why that stack has to scale too.
Which quantity should the study be run on?
The one you are going to report, and one at a time. Integrated quantities converge faster and more smoothly than local ones: lift and drag on an aerofoil will often show an observed order near 2 while the peak surface heat flux on the same three grids shows 0.8, because the peak sits on a near-singularity where the solution is not smooth and no scheme achieves its formal order. Both GCIs are correct and they are answers to different questions. If you need several, report several.
What GCI is acceptable?
There is no standard answer and this page does not pretend otherwise — ASME V&V 20 is deliberately silent on acceptance levels, because the acceptable discretisation uncertainty depends on what the number is being used for and on how it compares with the other uncertainties in the problem. The bands on this page are editorial guidance, not a standard. A useful test: compare the GCI with the spread you would get from a reasonable change of turbulence model, or with your experimental uncertainty. If the GCI is well below both, more cells are not the thing to spend on next.
Can I use GCI on an unsteady simulation?
On a time-averaged or peak quantity from a statistically stationary run, yes, with two extra cares. The averaging window must be long enough, and the same, on every grid — a different number of shedding cycles between grids produces apparent oscillatory convergence that has nothing to do with the mesh. And the time step should be refined together with the mesh, or held fine enough on all three that temporal error is negligible, otherwise the order you observe is a mixture of the spatial and temporal orders. The same procedure applied to time-step refinement alone gives a temporal GCI, and the arithmetic is identical with Δt in place of h.
Why does the calculator print the residual of the order equation?
Because the equation is transcendental in the general case and this page solves it numerically, so you are entitled to see how well. The residual is the apparent-order equation evaluated at the answer, in log space; it should sit around 1e-9. When it does not, it is because no positive order satisfies the equation for those three values at all — with unequal ratios the smallest ε₃₂/ε₂₁ any positive order can produce is ln r₃₂/ln r₂₁, so a measured ratio below that has no solution. That is information, not a numerical failure, and the page says so rather than printing the bracket edge as an order.
Is Richardson extrapolation the answer, then? Should I just report φ_ext?
No. φ_ext estimates the exact solution of the discretised equations on an infinitely fine grid, which is a different thing from the right answer — it removes discretisation error and nothing else. It is also less robust than the GCI: it depends on p in the numerator and denominator both, so an observed order that is 20% off moves φ_ext by more than it moves the GCI. The defensible way to report is the fine-grid value with the GCI as its uncertainty, mentioning φ_ext as a cross-check. If φ_ext lies outside φ₁ ± GCI, something in the study is wrong.
Related calculators
References
- P. J. Roache, Perspective: A Method for Uniform Reporting of Grid Refinement Studies, ASME Journal of Fluids Engineering 116 (1994), 405–413. The origin of the Grid Convergence Index and of the safety factors: 1.25 where the order has been observed over three or more grids, 3 where two grids force the formal order to be assumed. Cited, not reproduced.
- I. B. Celik, U. Ghia, P. J. Roache, C. J. Freitas, H. Coleman and P. E. Raad, Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications, ASME Journal of Fluids Engineering 130 (2008), 078001. The five-step procedure implemented here: the representative cell size h = (V/N)1/d, the apparent-order equation with its q(p) term, the extrapolation, the two relative errors and the convergence index, together with the recommendation that the refinement factor be at least 1.3 and the statement that a negative ε32/ε21 indicates oscillatory convergence. Clause-level citation only; no table is reproduced.
- H. W. Coleman, F. Stern, A. Di Mascio and E. Campana, The problem with oscillatory behavior in grid convergence studies, ASME Journal of Fluids Engineering 123 (2001), 438–439. A technical brief specifically on the case this page refuses to reduce to a single order; Coleman is also a co-author of the 2008 procedure above, so the oscillatory caveat in it is not an afterthought. Cited by title only.
- ASME V&V 20-2009, Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer. The framework the procedure above sits inside, and the source of the distinction between verification, which this page is about, and validation, which it is not. Copyrighted; cited by title and clause only.
- NASA Glenn Research Center, Examining Spatial (Grid) Convergence, NPARC Alliance CFD Verification and Validation web resource. A US Government work. Its worked example — three grids at a refinement ratio of 2 giving pressure recoveries of 0.97050, 0.96854 and 0.96178 — was used to verify this calculator: the page reproduces p = 1.786170, φext = 0.97130, GCI12 = 0.103083% and GCI23 = 0.356249% to every digit NASA prints.
- Verification performed for this page, not taken from a source (1): the apparent-order equation was solved for exact power-law data φi = φext + A hip0 over 40,000 random combinations of r21, r32 and p0, and the thirty-bisection solver recovered p0 to better than 2×10−5 in every case. The same data also settles the direction of the q(p) quotient, which is printed the wrong way round in at least one circulating transcription of the procedure: with the numerator and denominator swapped the recovered order is wrong by up to a factor of 400.
- Verification performed for this page, not taken from a source (2): the fixed-point iteration the procedure suggests for the apparent order has asymptotic contraction factor 1 − ln r32/ln r21, so it diverges for r32 > r212. It was confirmed to diverge numerically at (r21, r32) = (1.5, 4.0), (1.2, 2.5) and (1.1, 2.0) on data constructed to have an exact second-order solution. This is why the page bisects a monotone function instead.
- Derivation performed for this page, not taken from a source (3): the ratio r21p·GCI21/GCI32, widely reported as an asymptotic-range check, reduces algebraically to |φ2/φ1| for any three solutions and any refinement ratios, and φext computed from grids 2 and 3 is identically φext computed from grids 1 and 2 whenever convergence is monotone. Both were then confirmed numerically over 20,000 random cases to within the solver tolerance. Neither quantity carries information about whether the study is in the asymptotic range.
Setup guidance, not validation. Correlations have ranges of validity and cell-count estimates are order-of-magnitude. A converged simulation is not a correct one. Full disclaimer at calcengines.com/disclaimer/
