OpenFOAM: verification, validation and mesh independence

The CFD credibility literature (Roache; Oberkampf and Roy; the ASME V&V 20 standard) distinguishes three levels:

• Code verification: does the program solve the mathematical model correctly? It is checked with exact or manufactured solutions, measuring that the error decreases at the formal order of the schemes. It is the developers’ job; the user repeats it on a small scale with a case of known solution, such as Poiseuille.
• Solution verification: how large is the numerical error of this run? Iterative convergence, mesh and time-step studies, GCI. It is the core of this guide.
• Validation: does the model represent reality with the required accuracy? It needs reference data with their own uncertainty.

The order matters: without solution verification you do not know how much of the difference from the experiment is numerical error, and validation becomes a comparison between two numbers of unknown origin.

Integral quantities (flow rate, mean pressure drop, force on a wall) average local errors and usually converge quickly and smoothly. Local quantities, and especially derivatives (wall shear stress, gradients, local heat fluxes, the position of a reattachment point), converge more slowly and often non-monotonically, because they depend on the resolution in that specific region and on interpolation.

Practical consequence: a mesh “independent” for ΔP may not be so for wall shear stress. The study must be stated per QoI, not for the mesh in general.

With blockMesh change only the cell counts in the blocks entry of system/blockMeshDict, doubling the in-plane directions and keeping 1 cell through the thickness:

foamDictionary -entry blocks system/blockMeshDict     # read the current entry
# coarse:  hex (0 1 2 3 4 5 6 7) (20 10 1) simpleGrading (1 1 1)
# medium:  hex (0 1 2 3 4 5 6 7) (40 20 1) simpleGrading (1 1 1)
# fine:    hex (0 1 2 3 4 5 6 7) (80 40 1) simpleGrading (1 1 1)

With simpleGrading other than 1 the last-to-first cell ratio stays the same at every level: doubling the cells halves each cell approximately uniformly and the family remains geometrically similar, even though the ratio between adjacent cells changes. Grading is explained in Part 2: blockMeshDict.

Alternatively refineMesh splits every cell in two along the chosen directions, giving exactly r = 2; in 2D cases a system/refineMeshDict excluding the thickness direction is needed. To shorten runs on fine meshes you can initialise the field from the previous mesh’s solution:

cd ~/cases/case_M2_medium
mapFields ../case_M3_coarse -consistent -sourceTime latestTime

Initialisation only changes the starting point of the iterations: the final solution, if iterative convergence is complete, does not depend on it.

For a 2D domain of 0.100 × 0.010 m (area 1.0×10⁻³ m²): with 200 cells h = (1.0×10⁻³/200)^½ = 2.24 mm; with 800 cells 1.12 mm; with 3200 cells 0.559 mm. The ratio between successive values is 2.00, consistent with doubling cells in each direction. For the script the relative h (4, 2, 1) is enough: only the ratio matters.

for c in case_M1_fine case_M2_medium case_M3_coarse; do
  echo "== $c"
  grep -E 'cells:|Max non-orthogonality|Max skewness|Max aspect ratio|Mesh OK|\*\*\*' $c/log.checkMesh
done

In a well-built family the cell count scales as r^d and quality indicators stay similar or improve with refinement. If the fine mesh has much worse maximum non-orthogonality than the others, the family is not consistent: the study would compare different mesh qualities as well as different sizes.

A recurring practical rule in the verification literature (for example Roy, Oberkampf and the NASA tutorial on iterative convergence) is that the iterative error on the QoI should be at least an order of magnitude smaller than the discretisation error being estimated. With the differences above: if the fine and medium meshes differ by 0.5%, the fine-mesh QoI must be stable to about 0.05%.

A concrete way to estimate it: take ΔP over the last N steps and compute its relative variation and slope, as the monitoring script of guide 2 does. If residualControl tolerances stop the run before this level, tighten them.

The most repeatable way to measure the same quantities on three meshes is to let OpenFOAM do it during the run: the surfaceFieldValue function objects on inlet and outlet (area-weighted averages of p and sums of phi) are described in guide 2, with a script computing ΔP, ε_Q and stability over the final window. Copying the same system/functions into the three cases makes the measurement definition identical by construction.

For internal sections not coinciding with patches, the ParaView procedure of guide 3 (Slice, Integrate Variables) should be saved as a state or script and reapplied unchanged.

This guide measures y from the lower wall (0 ≤ y ≤ H); guide 3 measures it from the axis (−H/2 ≤ y′ ≤ H/2). With y = y′ + H/2 the two forms coincide:

6 U (y/H)(1 − y/H) = 6 U (½ + y′/H)(½ − y′/H) = 6 U (¼ − y′²/H²) = ³/₂ U [1 − (2y′/H)²]

When comparing ParaView CSVs, check where the mesh puts the zero of the transverse coordinate: an H/2 offset shifts the whole parabola and produces a huge RMSE that has nothing to do with the numerical solution.

For an order-of-magnitude check of the exercise: Re = ρ U_mean D_h / μ with D_h = 2H is 400 with ρ = 1000 kg/m³, and the estimate L_e ≈ 0.05 Re D_h gives about 0.4 m: with these data the 0.100 m stretch is developed only if the inlet already imposes the parabolic profile, or if the channel is preceded by a sufficient upstream length.

• Differences: f2 − f1 = 0.6, f3 − f2 = 2.4; same sign, so monotonic convergence.
• Ratio: 2.4/0.6 = 4; p = ln 4 / ln 2 = 2, equal to the formal order of a second-order scheme.
• r^p − 1 = 3; f_ext = 100.2 + (100.2 − 100.8)/3 = 100.0.
• GCI_12 = 1.25 × (0.6/100.2) / 3 × 100 = 0.2495%; GCI_23 = 1.25 × (2.4/100.8) / 3 × 100 = 0.9921%.
• Check: 0.9921 / (4 × 0.2495) = 0.994.

The same formulas in a spreadsheet, with f1, f2, f3 in B2:B4 and r in B5:

p      =LN(ABS((B4-B3)/(B3-B2)))/LN(B5)
f_ext  =B2+(B2-B3)/(B5^p-1)
GCI12  =1.25*ABS((B3-B2)/B2)/(B5^p-1)*100
GCI23  =1.25*ABS((B4-B3)/B3)/(B5^p-1)*100
check  =GCI23/(B5^p*GCI12)

The safety factor Fs proposed by Roache is 3 for two-grid studies and 1.25 with three or more grids and a verified observed order: the smaller value is the “reward” for having demonstrated the asymptotic range.

With different ratios r21 = h2/h1 and r32 = h3/h2, the procedure by Celik and co-authors (Journal of Fluids Engineering, 2008) obtains p by iterating:

p = | ln|ε₃₂/ε₂₁| + q(p) | / ln r₂₁,   q(p) = ln[ (r₂₁^p − s) / (r₃₂^p − s) ],   s = sign(ε₃₂/ε₂₁)

with ε21 = f2 − f1 and ε32 = f3 − f2. With r21 = r32 the q term vanishes and the guide’s formula is recovered. This is one reason why a constant-ratio family simplifies everything else.

With oscillatory convergence (s < 0) there is no meaningful observed order. A common conservative practice is to state the oscillation amplitude across the three meshes (maximum minus minimum) as the discretisation uncertainty, possibly adding a fourth mesh to see whether the asymptotic range is reached.

Complete script, Python standard library only (3.9 or later). Input: a CSV with columns mesh,cells,h,qoi and rows fine, medium, coarse. It writes results_gci.txt next to the input file.

#!/usr/bin/env python3
"""Guide 04: deliberately conservative equal-ratio, monotonic 3-grid GCI.

Usage: python gci_three_grids.py results_three_grids.csv
Input columns: mesh,cells,h,qoi ; h must be a characteristic size for
geometrically equivalent meshes. The 'mesh' labels must include fine,
medium,coarse (case-insensitive); 'cells' is recorded, not used to infer h.
Generates results_gci.txt next to input. Standard library only, Python >=3.9.
GPL-3.0-or-later example, 2026.
"""
from __future__ import annotations
import csv, math, sys
from pathlib import Path


def calculate(rows):
    required = {"fine", "medium", "coarse"}
    if set(rows) != required:
        raise ValueError("CSV must contain exactly fine, medium, coarse labels")
    x = [rows[key] for key in ("fine", "medium", "coarse")]
    for row in x:
        if row['cells'] <= 0 or row['h'] <= 0:
            raise ValueError("cells and h must be positive")
        if not all(math.isfinite(row[z]) for z in ('cells', 'h', 'qoi')):
            raise ValueError("All numerical values must be finite")
    (f1, f2, f3) = (p['qoi'] for p in x)
    (h1, h2, h3) = (p['h'] for p in x)
    if not h1 < h2 < h3:
        raise ValueError("Must have h(fine) < h(medium) < h(coarse)")
    r12, r23 = h2/h1, h3/h2
    if not math.isclose(r12, r23, rel_tol=0.002, abs_tol=0):
        raise ValueError("Unequal refinement ratios: this simple GCI model does not apply")
    r = (r12*r23)**0.5
    d21, d32 = f2-f1, f3-f2
    if d21 == 0 or d32 == 0:
        raise ValueError("Zero QoI difference: observed order is indeterminate")
    if d21*d32 <= 0:
        raise ValueError("Oscillatory or non-monotonic convergence: no simple GCI")
    if f1 == 0 or f2 == 0:
        raise ValueError("Zero QoI denominator: relative GCI is undefined")
    ratio = abs(d32/d21)
    if ratio <= 1:
        raise ValueError("No positive observed order: no defensible extrapolation")
    p = math.log(ratio)/math.log(r)
    denom = r**p - 1.0
    if not (0 < p < 10) or denom <= 0 or not math.isfinite(denom):
        raise ValueError("Observed order implausible / non-finite; inspect data")
    fext = f1 + (f1-f2)/denom
    fs = 1.25
    gci12 = fs * abs((f2-f1)/f1) / denom * 100
    gci23 = fs * abs((f3-f2)/f2) / denom * 100
    asym = gci23/(r**p*gci12)
    return dict(r=r,p=p,fext=fext,gci12=gci12,gci23=gci23,asym=asym,fs=fs,
                f1=f1,f2=f2,f3=f3, h1=h1,h2=h2,h3=h3)


def read_csv(path):
    out={}
    with path.open(newline='',encoding='utf-8-sig') as f:
        reader=csv.DictReader(f)
        cols=set(reader.fieldnames or ())
        if not {'mesh','cells','h','qoi'} <= cols:
            raise ValueError("Required CSV columns: mesh,cells,h,qoi")
        for row in reader:
            label=row['mesh'].strip().lower()
            if label in out:
                raise ValueError(f"Duplicate mesh label: {label}")
            out[label]={'cells':float(row['cells']), 'h':float(row['h']), 'qoi':float(row['qoi'])}
    return out


def report(a, path):
    return f'''OpenFOAM Guide 04 -- 3-grid GCI (training utility)
Input: {path.name}
Assumptions: identical geometry/physics/measurement, uniform refinement,
             monotonic asymptotic convergence, negligible iterative noise.
Fine/medium/coarse QoI: {a['f1']:.10g} / {a['f2']:.10g} / {a['f3']:.10g}
Characteristic sizes h: {a['h1']:.10g} / {a['h2']:.10g} / {a['h3']:.10g}
Refinement ratio r:          {a['r']:.8f}
Observed order p:            {a['p']:.8f}
Richardson estimate f(h=0): {a['fext']:.10f} [input QoI unit]
Fine-grid GCI12:             {a['gci12']:.6f} %
Medium-grid GCI23:           {a['gci23']:.6f} %
Asymptotic indicator:        {a['asym']:.6f} (near 1 expected)
Safety factor Fs:            {a['fs']:.2f}

CAUTION: GCI is not a total uncertainty budget, does not validate physics,
and is not a confidence interval. Inspect logs and all assumptions.
'''


def main(argv):
    if len(argv)!=2:
        print(__doc__,file=sys.stderr)
        return 2
    path=Path(argv[1])
    try:
        a=calculate(read_csv(path))
        text=report(a,path)
        out=path.with_name('results_gci.txt')
        out.write_text(text,encoding='utf-8')
        print(text)
        print(f'Written: {out}')
        return 0
    except (ValueError, OSError, KeyError, OverflowError) as ex:
        print(f'GCI NOT COMPUTED: {ex}',file=sys.stderr)
        return 1

if __name__=='__main__':
    sys.exit(main(sys.argv))

The test file results_three_grids.csv contains the data of the table above; the expected output is:

Refinement ratio r:          2.00000000
Observed order p:            2.00000000
Richardson estimate f(h=0): 100.0000000000 [input QoI unit]
Fine-grid GCI12:             0.249501 %
Medium-grid GCI23:           0.992063 %
Asymptotic indicator:        0.994048 (near 1 expected)
Safety factor Fs:            1.25

• With a fixed step (adjustTimeStep no) change only deltaT in controlDict and keep writeControl adjustableRunTime or runTime, so written times coincide across the three cases.
• With an adaptive step the parameter to halve is maxCo: the actual step varies in time and the ratio between steps is not exactly 2; read it from the log.
• With PIMPLE also fix nOuterCorrectors and the Final tolerances: each step must be iteratively converged, otherwise the study measures the splitting error rather than the temporal one.
• The same observed-order and GCI formulas apply with Δt in place of h, provided the compared signal is the same physical instant or the same statistic.

Time-control entries are explained in Part 4: controlDict, PIMPLE algorithms in Part 6.

The ASME V&V 20 standard (verification and validation in fluid dynamics and heat transfer) formalises the comparison as: comparison error E = S − D, the difference between simulation and data; validation uncertainty u_val, which combines numerical uncertainty, simulation-input uncertainty and experimental uncertainty. The quadrature combination u_val = √(u_num² + u_input² + u_D²) assumes the three contributions are independent: exactly the condition the guide asks you to check before using it.

If |E| is of the same order as u_val, model error cannot be distinguished from the “noise” of the uncertainties; if |E| is much larger, there is a model error of about E. In the fictitious example E = −1.0 Pa: without knowing u_val nothing more can be said.

For Poiseuille flow ΔP = 12 μ L U / H² is linear in μ and U: ±5% viscosity or mean velocity gives ±5% ΔP; ±5% of H gives about ∓10%, because ΔP depends on 1/H² (at fixed mean velocity). If the simulation returns sensitivities very different from these, the measured QoI includes other effects (development region, misplaced measurement sections) to understand before going on.

The normalised sensitivity coefficient S = (ΔQoI/QoI)/(Δx/x) sums up the response in one number: 1 for μ and U, about −2 for H here.

Bilingual template to fill in for every study (Markdown, opens in Notepad++):

# CFD Verification and Validation Record / Scheda di verifica e validazione

## Identification / Identificazione
- Date / Data:
- Operator / Operatore:
- OpenFOAM edition + version / Edizione e versione:
- Case path(s) on Linux server / Percorsi dei casi Linux:
- Windows review directory / Cartella di revisione Windows:
- Input revision or hash / Revisione o hash degli input:

## Physical problem / Problema fisico
- Geometry definition and source / Geometria e origine:
- Boundary and initial conditions / Condizioni al contorno e iniziali:
- Models, properties, units / Modelli, proprietà, unità:
- QoI(s) and measurement protocol / QoI e procedura di misura:
- Acceptance threshold and reason / Soglia accettazione e motivazione:

## Mesh family / Famiglia di mesh
| Grid | Cells | h | r | Worst relevant quality indicator | QoI |
|---|---:|---:|---:|---|---:|
| Fine | | | | | |
| Medium | | | | | |
| Coarse | | | | | |

## Iterative and conservation checks / Convergenza iterativa e bilanci
- Residual and QoI trends / Residui e trend delle QoI:
- Mass imbalance definition + value / Definizione e valore sbilancio massa:
- Time-step verification (if transient) / Studio deltaT (se transitorio):

## Mesh result / Esito di raffinamento
- Observed p / Ordine osservato:
- Richardson extrapolation / Estrapolazione:
- Fine-grid GCI [%]:
- Asymptotic consistency indicator / Indicatore asintotico:
- Non-monotonic / non-asymptotic issues / Eventuali problemi:

## Physical comparisons / Confronti fisici
- Benchmark or experiment and uncertainty / Benchmark o esperimento e incertezza:
- Input and geometry sensitivities / Sensibilità di input e geometria:
- Limitations / Limiti di validità:

## Conclusions / Conclusioni
- Suitability for intended use / Adeguatezza per lo scopo:
- Recommended additional work / Ulteriori verifiche necessarie: