Grid Convergence Index (GCI) calculator
Estimate the discretisation uncertainty of a CFD result from a grid convergence study on three meshes.
Settings
1Grids and solutions
The example is a 2D mesh. Is yours 3D?
Fine, medium and coarse follow the grid sizes, not the row order.
Paste a block from Excel or a Simcenter STAR-CCM+ table into any cell: one row per mesh, or one row per quantity.
Convergence plot
Show all values
Show the calculation with your numbers
2Result
3Report
4Next mesh for a target uncertainty
Simcenter STAR-CCM+ Base Size for the target
Scale the Base Size of the fine mesh in the automated mesh operation and keep the other sizes relative to the Base Size, so that only one parameter changes.
- Base Size
- Estimated cell count
- Set by
When the run is done, add it as a fourth grid. The calculator then uses the three finest grids and checks the order against the coarsest. To run the meshes in one go, make the Base Size a parameter and let a Sweep study in Design Manager change it. Our post on how to run Design Manager projects covers the set-up.
The estimate covers the discretisation error of the quantities you enter, not modelling or iteration errors. Volupe provides the calculator free of charge and without warranty.
The grid convergence index (GCI) puts a band on the discretisation error of a CFD result, an estimate of how far the value on the finest of three systematically refined meshes may lie from the value at zero cell size. Roache proposed it in 1993 and published it in 1994. Celik et al. (2008) wrote down the procedure the ASME Journal of Fluids Engineering recommends, with GCI = 1.25 ea/(rp − 1), where ea is the relative change between the two finest meshes, r the refinement ratio and p the observed order.
How to use the GCI calculator
1. Grids and solutions. The calculator opens with the first column of Celik’s Table 1. Clear it and enter your meshes as cell counts or Base Sizes, one row per mesh and one column per quantity, or paste a block from a spreadsheet. Choose 2D or 3D by the directions you refine in. A one-cell-thick mesh refined in-plane is 2D.
2. Result. Each quantity gets φ1 ± U, the value to report, with its class of convergence and Celik’s GCI next to it. The checks under it tell you what to look at, such as an observed order far from the nominal one or too few digits.
3. Report. Copy the result as a paragraph, a table, LaTeX or CSV, with every input, the rule and the factor of safety. The link restores the whole study for a colleague.
4. Next mesh. If U is too large for what the result is used for, enter a target. The calculator gives the Base Size and cell count of the next mesh, with a range.
What the GCI tells you
The GCI puts a band on the discretisation error, which is the difference between your fine-mesh result and the exact solution of the same equations. It says nothing about whether those equations describe the real flow. A wrong turbulence model or a domain that is too small can converge just as nicely, towards the wrong answer, with a small GCI.
A band of ± 0.5 % means that further refinement is unlikely to move the result by more than about half a per cent.
A grid study belongs to verification, which Stern et al. (2001) define as assessing the numerical uncertainty of a simulation. Validation, the comparison with benchmark experimental data, assesses the modelling uncertainty and is a separate step (see verification and validation).
Mesh independence study: how many meshes do you need?
A mesh independence study runs the same case on successively finer meshes and checks how much the result still changes. It is also called a grid independence study, a mesh convergence study, a grid convergence study, a mesh sensitivity study or a grid refinement study. The names are used for the same work, and the GCI is one recommended way of reporting it.
You need at least three meshes. With two you only see that the result is sensitive to the mesh, not that it converges, and you have to assume the order of convergence. That is why the factor of safety Roache recommends is 3 for two meshes and 1.25 for three or more.
The study in five steps
1. Make a fine, a medium and a coarse mesh by changing one global size by the same ratio, above 1.3.
2. Keep the prism layers, physics, solver settings and geometry the same.
3. Run each mesh until the residuals have dropped three to four orders of magnitude and the quantities you report have stopped changing.
4. Enter the meshes and the results in the calculator above.
5. If U is too large for what the result is used for, run the finer mesh that step 4 suggests.
Grid 1 fine
cell size h1
N1 ≈ 530
Grid 2 medium
h2 = 1.41 h1
N2 ≈ 260
Grid 3 coarse
h3 = 2 h1
N3 ≈ 130
N ∝ h−Dcells per grid, D = 2 or 3
| per level | grid 3 to grid 1 | |
|---|---|---|
| 2D | × 2 | × 4 |
| 3D | × 2.83 | × 8 |
Systematic refinement and the ratio
Systematic means that you change one global parameter, such as the Base Size in Simcenter STAR-CCM+, by a constant factor between the levels, instead of remeshing with different settings each time. On an unstructured mesh the refinement itself should still be structured, and geometrically similar cells are preferable (Celik et al. 2008).
Use a refinement ratio above 1.3 between neighbouring meshes. That is Celik’s advice, which the paper says comes from experience, not from a formal derivation. The NASA NPARC tutorial accepts ratios down to 1.1, but the calculator flags anything below 1.3. A ratio of √2 ≈ 1.41 multiplies the cell count by about 2.8 per level in 3D (2 in 2D), and most of the cost of an 8 M / 2.8 M / 1 M cell study sits in the fine mesh.
Keep everything else constant
Prism layers. Keep the number of layers, the stretching and the total thickness the same on all three meshes, so that only the core mesh changes. In Simcenter STAR-CCM+ the Prism Layer Total Thickness can be set relative to the Base Size, and then it changes with every Base Size. Set it as an absolute value for the study. You then size the first cell once, and the y+ calculator gives its height for a target y+. Note that the GCI then covers the core mesh only, not the resolution normal to the wall.
Physics and solver settings. Use the same turbulence model, boundary conditions, discretisation schemes, time step and convergence criteria on every mesh.
Geometry and domain. The CAD, the extent of the domain and every geometric feature stay the same. Only the characteristic cell size changes.
Prefer integral quantities
Use quantities that depend on a large part of the mesh, such as a drag or lift coefficient, a pressure drop, a mass flow or an area-averaged temperature. Local quantities converge much more slowly on coarser meshes (Phillips and Roy 2014). Quantities converge at different rates, and you should check each one you will report. The calculator takes several at once and names the one with the largest uncertainty.
Converge each run first
Before you trust the three results, check that each run converged. Celik et al. ask for the normalised residuals to fall by at least three, and preferably four, orders of magnitude, with the iteration error at least an order of magnitude below the discretisation error. Residuals alone are not enough, as Christoffer Johansson writes in our post on residuals in Simcenter STAR-CCM+. We watch the quantities we report until they stop changing, and average them over the last part of the run if they fluctuate.
The GCI method and formula, step by step
The calculator follows the five steps of Celik et al. (2008) with five changes. It solves for the observed order with a bracketed method, classifies by an existence test instead of the bands of R, lets you choose the factor-of-safety rule and, for oscillating results, reports half the range of the solutions when that is larger (Stern et al. 2001). When the digits you entered leave the rule’s choice open, it also takes the larger U. Mesh 1 is always the finest, h1 < h2 < h3. With cell counts as the input, D sets the ratios in (1) and with them the observed order.
h = (V/N)1/D, r21 = h2/h1, r32 = h3/h2 (1)
ε21 = φ2 − φ1, ε32 = φ3 − φ2, s = sign(ε32/ε21) (2)
The changes between the meshes. s = −1 means that the solution alternates, which the paper takes as a sign of oscillatory convergence.
p ln r21 = ln|ε32/ε21| + ln[(r21p − s)/(r32p − s)] (3)
The observed order p, Celik’s equation 3. The last term vanishes when r21 = r32. The paper writes the right-hand side inside an absolute value, which keeps p positive. We drop it and test (4) first, because the absolute value can give a positive p where none fits.
ε32/ε21 > ln r32/ln r21 if s = +1, |ε32/ε21| > 1 if s = −1 (4)
The existence test. Equation (3) has a positive root only when (4) holds. With equal ratios it is 0 < R < 1 or −1 < R < 0, where R = ε21/ε32.
Fixed-point iteration, which the paper suggests for (3), can fail to converge when r32 is larger than r21, and then returns a p that does not satisfy (3). The calculator checks (4) first, brackets the root, halves the bracket until p is fixed and accepts it only if it satisfies (3). In a million exact test cases, monotonic and oscillatory with equal and unequal ratios, it gave no wrong result.
Richardson extrapolation
Richardson extrapolation combines the solutions on two meshes into an estimate of the value at zero cell size. Richardson used it in a 1911 paper on finite differences, applied to the stresses in a masonry dam. With the observed order from (3), the procedure extrapolates from the two finest meshes.
φext21 = (r21pφ1 − φ2)/(r21p − 1) (5)
Richardson extrapolation. φext32 follows in the same way from the medium and coarse meshes.
Note that the extrapolation is only accurate in the asymptotic range, when the meshes are fine enough for the error to fall as hp (Phillips and Roy 2014). Across shocks and other discontinuities it does not hold locally, as the NPARC tutorial points out. For oscillatory results the calculator shows φext only for traceability, and when (4) fails it shows none.
ea21 = |(φ1 − φ2)/φ1|, eext21 = |(φext21 − φ1)/φext21| (6)
U = Fs|ε21|/(r21p − 1) (7)
The band on the fine-mesh result in the units of φ, and U/|φ1| gives it in per cent. With Fs = 1.25 and the observed p this is Celik’s GCIfine21 = 1.25ea21/(r21p − 1).
Which Fs and which p go into (7) depends on the rule, and the calculator has three. Whichever you choose, it also shows Celik’s value, since that is the Journal of Fluids Engineering procedure and the value to compare with published studies.
The default, Oberkampf and Roy (2010). If the observed order is within 10 % of the nominal order pf, (7) uses Fs = 1.25 with pf in place of p. Otherwise it uses Fs = 3 with min(max(0.5, p), pf). The order is then never taken above the nominal one or below 0.5, and U more than doubles at the edge of the band. We use the rule as Phillips (2012) writes it out.
The paper’s rule, Celik et al. (2008). Fs = 1.25 with the observed p, whatever its value. For an oscillating quantity the calculator reports the larger of this and half the range, and the headline can then exceed the paper’s value. The card’s Celik line always shows Celik’s GCI, for column 3 of the example 0.472 %, which the paper rounds to 0.5 %.
The Xing and Stern (2010) rule
The factor of safety depends on P = p/pf, with Fs = 1.6P + 2.45(1 − P) for P ≤ 1 and Fs = 1.6P + 14.8(P − 1) for P > 1. The observed p, but not below 0.5, goes into (7). We take the coefficients from Phillips (2012), who also adds the 0.5 floor for consistency.
hnew/h1 = (Utarget/U)1/p, Nnew = N1(h1/hnew)D (8)
The next mesh, from U ∝ hp with the p that set U. The calculator also gives the range you get with the smallest and the largest of that p, the observed p and the nominal order.
Assumptions: systematic refinement, iteratively converged solutions and meshes in the asymptotic range. Percentages are relative to |φ1| unless you give a reference value, the better choice when φ1 is close to zero.
Read the results
Convergence type
The calculator sorts each quantity into one of four classes (monotonic convergence, oscillatory convergence, no positive order and divergent) from the signs of ε21 and ε32 and the existence test (4). The ITTC procedure, based on Stern et al. (2001), uses bands of the convergence ratio R = ε21/ε32 instead, and they assume equal refinement ratios. Say the drag coefficient is 0.32000, 0.34048 and 0.36000 on 4 096 000, 1 000 000 and 455 166 cells (ratios 1.60 and 1.30). R is 1.05 and the bands would say divergent, yet these values converge monotonically with p = 1.5.
Monotonic convergence
Same sign, shrinking: 0 < R < 1
Report: U from the rule, with Celik’s GCI next to it
Oscillatory convergence
Sign changes, shrinking: −1 < R < 0
Report: U = max(Urule, ½(φmax − φmin))
No positive order
Same sign, not shrinking: R ≥ 1
Report: No GCI; a rough band for reference
Divergent
Sign changes, not shrinking: R ≤ −1
Report: No uncertainty estimate
Unequal ratios: R is not the test; same-sign changes converge only if ε32/ε21 > ln r32 / ln r21.
| Class | Test | What it means | What the calculator reports |
|---|---|---|---|
| Monotonic convergence | s = +1 and (4) holds. With equal ratios 0 < R < 1 | The changes have the same sign and shrink fast enough for a positive order | U from the rule, with Celik’s GCI next to it |
| Oscillatory convergence | s = −1 and (4) holds, so |ε32| > |ε21|. With equal ratios −1 < R < 0 | The solution alternates from mesh to mesh, and the changes shrink | The larger of U and ½(φmax − φmin). A fourth mesh helps |
| No positive order | s = +1 and (4) fails. With equal ratios R ≥ 1 | The changes have the same sign but do not shrink fast enough for any p > 0. With a constant ratio the solution moves away from mesh to mesh | No GCI. A rough band with Fs = 3 and p = 0.5 (Oberkampf and Roy) for reference. Check iterative convergence and the mesh, or add a finer mesh |
| Divergent | s = −1 and (4) fails, so |ε32| ≤ |ε21|. With equal ratios R ≤ −1 | The solution alternates and the changes do not shrink | No uncertainty. Check iterative convergence, mesh topology and boundary conditions |
Observed order and the factor of safety
Compare the observed order p from (3) with the nominal order of the scheme, pf, which is 2 for a second-order scheme. Agreement between the two is a good sign that the meshes are in the asymptotic range, while a mismatch need not mean poor calculations (Celik et al. 2008).
The mismatch does matter for the uncertainty. An observed order above the nominal one makes rp − 1 in (7) large and the band small, which Xing and Stern (2010) call unreasonably small uncertainty, and Roache later capped p at the nominal order for the same reason. This is why we use Oberkampf and Roy by default.
What to report
Report the fine-mesh value with its band, φ1 ± U, and state the rule and the factor of safety. How small the band has to be depends on what the result is used for, and we have not found a source that gives one number. If you compare two design variants, the difference between them should be clearly larger than their uncertainties. If it is not, a finer mesh could change which variant comes out ahead. The calculator marks U above 10 % of φ1 in amber. That is our prompt to look at the next mesh, not a limit from the literature.
The asymptotic range
The meshes are in the asymptotic range when the error falls as hp, so that the next refinement changes the result by the amount the observed order predicts.
Why there is no asymptotic-range ratio
The NASA NPARC tutorial checks the asymptotic range with the ratio GCI32/(r21p GCI21), which should be close to 1. When p is the observed order from the same three meshes, the ratio is always |φ1/φ2|, whatever the meshes. NPARC’s own worked example gives 1.002019, and 0.97050/0.96854 is 1.002024. A 1.5 % change between the two finest meshes always gives about 1.015, so the check cannot fail, and we have left it out.
What the calculator does check is P = p/pf, which sets the factor of safety, an optional fourth mesh and the reported digits. A fourth mesh gives a second observed order from the three coarsest meshes, and the calculator tells you whether the two agree within 10 %. For monotonic results it also extrapolates meshes 1–2 and 2–3 with the nominal order. In the asymptotic range the two estimates agree, and if they differ by more than U (our criterion), the meshes are probably not there yet.
Illustrative model, six grids with r = 2. Dashed teal: the nominal order, pf = 2. Dashed ochre: an iteration or round-off floor.
Reported digits
Small differences between the meshes can drown in the rounding of the numbers you copy. Exact second-order data at r = √2, copied with four significant digits as 0.3012, 0.3024 and 0.3048, give an observed order of 2.00, but within the rounding anything from 1.65 to 2.37. Under Oberkampf and Roy that decides the factor of safety, because 2.00 is inside the 10 % band and both ends are outside it. The calculator then reports the largest U the typed digits allow, 1.55 %, with the 0.50 % at the observed order next to it. That is our choice, not part of the Oberkampf and Roy rule. With Xing and Stern it does the same when the digits allow a U more than 1.5 times larger, and step 4 then uses the order that gave the larger U. If a difference is no larger than the rounding, the quantity is not resolved and gets no uncertainty. The fix is to export more digits from your solver.
Worked example: Celik et al. (2008), Table 1
The example in the paper uses data from Celik and Karatekin (1997), who computed the 2D turbulent flow over a backward-facing step on non-uniform structured grids. Column 1 of Table 1 is the dimensionless reattachment length on 18 000, 8 000 and 4 500 cells. The calculator opens with it, and we check it by hand below. Load this example into the calculator.
In 2D, r21 = (18 000/8 000)1/2 = 1.5 and r32 = (8 000/4 500)1/2 = 1.333. The changes are ε21 = 5.972 − 6.063 = −0.091 and ε32 = −0.109, with the same sign, and 0.109/0.091 = 1.198 is larger than ln 1.333/ln 1.5 = 0.710, so (4) holds. Solving (3) gives p = 1.534 and r21p = 1.8626, and (5) gives φext = (1.8626 × 6.063 − 5.972)/0.8626 = 6.1685. The relative change is ea = 0.091/6.063 = 1.5009 %.
With Celik’s Fs = 1.25, (7) gives 1.25 × 1.5009 %/0.8626 = 2.17 %, which the paper rounds to 2.2 %. The observed order is 23 % below the nominal 2, outside the 10 % band, so Oberkampf and Roy use Fs = 3 with the same p. That gives 3 × 1.5009 %/0.8626 = 5.22 %, or 6.063 ± 0.316. Xing and Stern, with P = 0.77, give Fs = 1.80 and 3.13 %.
| Value | Paper | Calculator |
|---|---|---|
| Observed order p | 1.53 | 1.534 |
| φext | 6.1685 | 6.1685 |
| GCI (Celik) | 2.2 % | 2.17 % |
| U (O&R) | – | 5.22 % |
| Value | Paper | Calculator |
|---|---|---|
| Observed order p | 0.75 | 0.752 |
| φext | 10.8801 | 10.8801 |
| GCI (Celik) | 1.1 % | 1.07 % |
| U (O&R) | – | 2.56 % |
| Value | Paper | Calculator |
|---|---|---|
| Observed order p | 1.51 | 1.508 |
| φext | 6.0269 | 6.0269 |
| GCI (Celik) | 0.5 % | 0.472 % |
| U (O&R) | – | 1.13 % |
Celik et al. (2008), Table 1, against the calculator. GCI is Celik’s value (Fs = 1.25). U (O&R) is the Oberkampf and Roy (2010) uncertainty in per cent of φ1, which the paper does not give. Velocities in m/s.
The other two columns are the axial velocity at x/H = 8 on 18 000, 4 500 and 980 cells. The first point converges monotonically with an observed order of 0.75. The second oscillates, since the medium mesh gives a lower velocity than both the others. There the rule gives 1.13 % and half the range 1.07 %, and the calculator reports 1.13 %, against 0.47 % for the GCI. Load both quantities into the calculator.
About the paper’s Table 1
The paper’s heading repeats y = 0.0526 for the third column. Its velocities, around 6 m/s against 10.8 m/s at the first point, show that it is a different point on the same profile (in Celik’s Fig. 1, y runs from 0 to 0.12 m). We call it the second point. The paper also prints eext = 0.9 % for the second column, while its own numbers give 0.85 %. Apart from that, the calculator reproduces p, φext and the GCI to the digits the paper prints.
Practical tips and pitfalls
When φ is close to zero, give a reference value
A lift coefficient near zero incidence makes the relative error explode even when the absolute changes are tiny. When φ1 is zero or very small, the NPARC tutorial notes that the fractional error may be meaningless and that another normalising value should be used. In the calculator you can give a reference value per quantity, such as the matching drag coefficient or the largest value of the quantity in the domain. The band U in the quantity’s own units does not depend on that choice.
Blended and hybrid schemes count as first order
If your scheme blends or switches between first and second order, such as the hybrid, power-law and related exponential schemes, set the nominal order to 1. The Journal of Fluids Engineering policy (Freitas 1993) treats such schemes as first-order methods unless you can show that their numerical diffusion does not swamp the modelled physical diffusion. Celik’s own profile example used a hybrid scheme, and its average local order was 1.49. Our worked example on the finite volume method shows the same on a small case: first-order upwind reaches an observed order of 0.93, central differencing 1.98.
Transient simulations need two studies
The GCI covers the spatial discretisation. In a transient case we first find a small enough time step on one mesh and then run the grid study with that time step fixed. If both change at once, you cannot tell which of them moved the result. The JFE policy asks for temporal accuracy to be demonstrated, and Celik et al. ask for iterative convergence to be checked at every time step, not only the last one.
Unstructured meshes need structured refinement
For polyhedral, tetrahedral or trimmed meshes, h = (V/N)1/D is an average, and the NPARC tutorial uses the same effective ratio for unstructured grids. The refinement itself should still be structured. In Simcenter STAR-CCM+ that means changing only the Base Size and leaving the sizes that are set relative to it as they are. If we see oscillatory convergence or an unexpected p on such a mesh, the first thing we check is whether the mesh topology changes between the levels.
Inlet sensitivity is a separate uncertainty
Small changes in the inlet conditions can move the result as well. Say a 10 % change in inlet turbulence intensity moves the drag by 5 %. Celik et al. recommend studying and reporting this sensitivity, and it comes on top of the GCI. The turbulence boundary conditions calculator helps you choose defensible inlet values.
Profiles get error bars with one average order
For a velocity or temperature profile the GCI is computed point by point, but with one average order pave as the global order of accuracy (Celik et al. 2008). In the paper’s example the local order ranged from 0.012 to 8.47 with an average of 1.49, and 20 % of the 22 points oscillated. The result is plotted as error bars on the fine-mesh profile, in the same way as experimental uncertainty. The calculator works on single values, not on profiles.
Noisy convergence calls for least squares
When the convergence is noisy, the paper suggests the least-squares version of the GCI. It fits the error model to more than three meshes, at least four, and Eça and Hoekstra (2014) describe the procedure. The calculator does not do the least-squares fit, but it accepts a fourth mesh as a check on the order.
A grid convergence study in Simcenter STAR-CCM+
This is how we set up the three meshes in Simcenter STAR-CCM+ and get the numbers into the calculator.
- 1
Base Size B
Automated Mesh gives grid 1, for example at 2 mm
- 2
B × √2 and B × 2
Grids 2 and 3 (2.83 and 4 mm), with nothing else changed
- 3
Run to convergence
Each mesh, with report monitors on the quantities
- 4
Read the results
Base Size, cell count and report values of each mesh
- 5
Paste the rows
Into the calculator’s grid table, with Base Size as the input; keep the cell counts for the report
- 6
Next Base Size
The calculator’s step 4 gives it if U is above your target
One parameter. In the automated mesh operation, sizes can be set relative to the Base Size or as absolute values. Relative sizes follow the Base Size and absolute ones do not. For a fine Base Size of 2 mm, use 2.83 mm for the medium mesh and 4 mm for the coarse one, a ratio of √2. Check your custom controls before you start, so that the sizes you want refined are relative.
Prism layers. With the same absolute prism layer settings on all three meshes, the number of prism cells follows the wall area, about 1/h2 rather than 1/h3, and the cell-count ratios no longer match the Base Size ratios. In a model case with 40 % prism cells on the fine mesh, data that are exactly second order in the Base Size give an observed order of 2.51 from the cell counts. We therefore recommend the Base Sizes as the input, and the cell counts for the report.
Quantities. Create reports for the quantities you will report and plot them while the solver runs, so that you see when they stop changing. Enter the converged values with all the digits the report gives.
A sweep in Design Manager. Make the Base Size a simulation parameter and let a Design Manager sweep study change it, as Florian Vesting does for the Backward Facing Step tutorial in our post on mesh control and uncertainty analysis. How to run Design Manager projects covers running the study on a workstation or a cluster.
The next mesh. Step 4 gives the Base Size and cell count for a target uncertainty. After the run, add it as a fourth grid. The calculator uses the three finest and checks the order against the old coarse mesh. For large, expensive models our application engineers can help you plan the levels.
References and changelog
References
I. B. Celik, U. Ghia, P. J. Roache and C. J. Freitas, “Procedure for estimation and reporting of uncertainty due to discretization in CFD applications”, Journal of Fluids Engineering 130(7), 078001, 2008. doi:10.1115/1.2960953
P. J. Roache, “Perspective: a method for uniform reporting of grid refinement studies”, Journal of Fluids Engineering 116(3), 405–413, 1994. doi:10.1115/1.2910291
F. Stern, R. V. Wilson, H. W. Coleman and E. G. Paterson, “Comprehensive approach to verification and validation of CFD simulations, Part 1: methodology and procedures”, Journal of Fluids Engineering 123(4), 793–802, 2001: convergence ratio and oscillatory uncertainty. doi:10.1115/1.1412235
W. L. Oberkampf and C. J. Roy, Verification and Validation in Scientific Computing, Cambridge University Press, 2010: the default factor-of-safety rule. doi:10.1017/CBO9780511760396
More references (11)
I. Celik and O. Karatekin, “Numerical experiments on application of Richardson extrapolation with nonuniform grids”, Journal of Fluids Engineering 119(3), 584–590, 1997: the data of the worked example. doi:10.1115/1.2819284
P. J. Roache, Verification and Validation in Computational Science and Engineering, Hermosa Publishers, 1998: factors of safety 3 for two grids and 1.25 for three or more, as quoted by NASA NPARC and ITTC 7.5-03-01-01.
ITTC, “Uncertainty analysis in CFD verification and validation, methodology and procedures”, Recommended Procedures and Guidelines 7.5-03-01-01, Revision 03, 2017. ittc.info
P. J. Roache, K. N. Ghia and F. M. White, “Editorial policy statement on the control of numerical accuracy”, Journal of Fluids Engineering 108(1), 2, 1986. doi:10.1115/1.3242537
C. J. Freitas, “Journal of Fluids Engineering editorial policy statement on the control of numerical accuracy”, Journal of Fluids Engineering 115(3), 339–340, 1993: blended and hybrid schemes. doi:10.1115/1.2910144
L. F. Richardson, “The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam”, Philosophical Transactions of the Royal Society A 210, 307–357, 1911. doi:10.1098/rsta.1911.0009
L. Eça and M. Hoekstra, “A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies”, Journal of Computational Physics 262, 104–130, 2014: least-squares GCI. doi:10.1016/j.jcp.2014.01.006
T. S. Phillips, Extrapolation-based Discretization Error and Uncertainty Estimation in Computational Fluid Dynamics, MS thesis, Virginia Tech, 2012: the Oberkampf and Roy and the Xing and Stern rules as used here. vtechworks.lib.vt.edu
T. S. Phillips and C. J. Roy, “Richardson extrapolation-based discretization uncertainty estimation for computational fluid dynamics”, Journal of Fluids Engineering 136(12), 121401, 2014. doi:10.1115/1.4027353
T. Xing and F. Stern, “Factors of safety for Richardson extrapolation”, Journal of Fluids Engineering 132(6), 061403, 2010. doi:10.1115/1.4001771
NASA Glenn Research Center, NPARC Alliance, “Examining spatial (grid) convergence”, a tutorial with a worked example. grc.nasa.gov
Review
Reviewed by Christoffer Johansson, M.Sc., Senior Application Specialist at Volupe.
Changelog
October 2026: new calculator with a bracketed solve and an existence test, Oberkampf and Roy by default with Celik’s GCI beside it, and checks on a fourth mesh and the reported digits. Checked against Table 1 of Celik et al. (2008).
Before and after a grid study
All CFD Handbook tools →y+ = uτy / ν
Calculator
y+ calculator
Size the first cell for a target y+, then the prism layer stack that stays fixed across your grids.
Open calculator →
N = mk
Planner
DOE calculator
Estimate how many simulation runs a design study needs once the mesh is settled.
Open planner →
k = ³⁄₂ (U I)2
Calculator
Turbulence boundary conditions calculator
Turn turbulence intensity into inlet k, ε, ω and ν̃.
Open calculator →
Questions about the GCI
Calculator reviewed by
Christoffer Johansson
Senior Application Specialist, M.Sc.
Ask about grid studies. We reply within one day.
Ask a technical question →
What is the Grid Convergence Index (GCI)?
GCI = Fs × ea / (r^p − 1), where ea is the relative change between the two finest meshes, r the refinement ratio, p the observed order and Fs a factor of safety, 1.25 for a study on three or more meshes. Roache proposed it in 1993 and published it in 1994, and Celik et al. (2008) give the procedure the ASME Journal of Fluids Engineering recommends. It is a band around the fine-mesh result that estimates how far that result is from the value at zero cell size.
How many meshes do I need for a mesh independence study?
At least three, refined systematically. Three meshes let you compute the observed order and compare it with the nominal one. With two you have to assume the order, and the recommended factor of safety rises from 1.25 to 3. A fourth mesh checks the observed order on a second set of three, and the least-squares GCI needs at least four.
What refinement ratio should I use?
A ratio above 1.3 between neighbouring meshes, which Celik et al. (2008) recommend from experience. In 3D a ratio of √2 ≈ 1.41 multiplies the cell count by about 2.8 per level, and a ratio of 2 by 8. The two ratios do not have to be equal. The calculator solves for the observed order with any pair of ratios and first checks that a positive order exists.
What GCI value is acceptable?
No source we know of gives one number. Judge the band by what the result is used for, and report the fine-mesh value with its band, φ1 ± U, and the rule you used. If you compare two design variants, the difference between them should be clearly larger than their uncertainties. In their own example Celik et al. (2008) report 2.2 % for the reattachment length and add that this does not account for modelling errors.
How do I know the meshes are in the asymptotic range?
The meshes are in the asymptotic range when the error falls as h^p. The best signs are an observed order close to the nominal one (Celik et al. 2008) and a fourth mesh that gives the same order as the first three. The NASA NPARC tutorial uses the ratio GCI32 / (r21^p × GCI21) as a check, but when p is fitted to the same three meshes it always equals |φ1/φ2| and cannot tell you anything new.
What is Richardson extrapolation?
A way of estimating the value at zero cell size from solutions on two meshes, first used by Richardson (1911). With the observed order p and the refinement ratio r21, φext = (r21^p × φ1 − φ2) / (r21^p − 1). It is only accurate in the asymptotic range, so be careful with φext when the observed order is far from the nominal one or the solution oscillates.
Can Volupe help me set up a grid study in Simcenter STAR-CCM+?
Yes. Our application engineers help Simcenter STAR-CCM+ users plan mesh studies, set up systematic refinement and interpret the results, as part of support or in a training course. Use the form below and choose Technical question.
Talk to an engineer in your region
Every Volupe business developer holds an engineering degree. Find your regional contact below, or send us a message.
Ask us about Grid Convergence Index (GCI) calculator
Use this form for licence quotes, demos, training and technical questions.
What happens next
- A reply within one day.
- A scoping call with one of our engineers.
- A quote, demo or other next step that fits your enquiry.
After purchase: support from Volupe simulation engineers.