OpenFOAM: advanced numerical methods, convergence and performance

Before comparing schemes and timings it pays to know what is being computed. extra/case_estimates.py reads the resources/blockMesh/pitzDaily dictionary (the one used by run_case.sh), applies the blockMesh grading rules, including multi-section gradings such as $posY, and reads the inlet conditions of the case. The main numbers:

The script and its run on the prepared case are in the section 7 box (mesh quality); the tutorial inlet value of ω, 440.15 s⁻¹, equals ε/(C_μk): the inlet conditions of the different turbulence models are consistent with each other.

File contents (without headers), to be compared with the copy on the server. The complete case is in the version 14 repository. system/blockMeshDict is a copy of tutorials/resources/blockMesh/pitzDaily (as the comment in the file says): Allrun and run_case.sh use the resource version with -dict.

Allrun

#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

runApplication blockMesh -dict $FOAM_TUTORIALS/resources/blockMesh/pitzDaily
runApplication foamRun

#------------------------------------------------------------------------------

0/U

dimensions      [velocity];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (10 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    upperWall
    {
        type            noSlip;
    }

    lowerWall
    {
        type            noSlip;
    }

    frontAndBack
    {
        type            empty;
    }
}

0/p

dimensions      [kinematicPressure];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    upperWall
    {
        type            zeroGradient;
    }

    lowerWall
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

0/k

dimensions      [turbulentKineticEnergy];

internalField   uniform 0.375;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 0.375;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            kqRWallFunction;
        value           uniform 0.375;
    }
    lowerWall
    {
        type            kqRWallFunction;
        value           uniform 0.375;
    }
    frontAndBack
    {
        type            empty;
    }
}

0/epsilon

dimensions      [turbulentEpsilon];

internalField   uniform 14.855;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 14.855;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            epsilonWallFunction;
        value           uniform 14.855;
    }
    lowerWall
    {
        type            epsilonWallFunction;
        value           uniform 14.855;
    }
    frontAndBack
    {
        type            empty;
    }
}

0/nut

dimensions      [turbulentViscosity];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            calculated;
        value           uniform 0;
    }
    outlet
    {
        type            calculated;
        value           uniform 0;
    }
    upperWall
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    lowerWall
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    frontAndBack
    {
        type            empty;
    }
}

0/omega, 0/nuTilda, 0/v2, 0/f (used only by other turbulence models)

dimensions      [turbulentOmega];

internalField   uniform 440.15;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           $internalField;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            omegaWallFunction;
        value           $internalField;
    }
    lowerWall
    {
        type            omegaWallFunction;
        value           $internalField;
    }
    frontAndBack
    {
        type            empty;
    }
}

dimensions      [turbulentViscosity];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    outlet
    {
        type            zeroGradient;
    }

    upperWall
    {
        type            zeroGradient;
    }

    lowerWall
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

dimensions      [turbulentKineticEnergy];

internalField   uniform 0.25;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           $internalField;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            v2WallFunction;
        value           $internalField;
    }
    lowerWall
    {
        type            v2WallFunction;
        value           $internalField;
    }
    frontAndBack
    {
        type            empty;
    }
}

dimensions      [rate];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            fWallFunction;
        value           uniform 0;
    }
    lowerWall
    {
        type            fWallFunction;
        value           uniform 0;
    }
    frontAndBack
    {
        type            empty;
    }
}

constant/physicalProperties, constant/momentumTransport

viscosityModel  constant;

nu              1e-05;

simulationType RAS;

RAS
{
    // Tested with kEpsilon, realizableKE, kOmega, kOmega2006, kOmegaSST, v2f,
    // ShihQuadraticKE, LienCubicKE.
    model           kEpsilon;

    turbulence      on;


    viscosityModel  Newtonian;
}

system/controlDict, system/functions

solver          incompressibleFluid;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         2000;

deltaT          1;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

#includeFunc streamlinesLine
(
    name=streamlines,
    start=(-0.0205 0.001 0.00001),
    end=(-0.0205 0.0251 0.00001),
    nPoints=10,
    fields=(p k U)
)

system/fvSchemes

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,k)      bounded Gauss limitedLinear 1;
    div(phi,epsilon) bounded Gauss limitedLinear 1;
    div(phi,omega)  bounded Gauss limitedLinear 1;
    div(phi,v2)     bounded Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
    div(nonlinearStress) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

system/fvSolution

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.1;
        smoother        GaussSeidel;
    }

    pcorr
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0;
        smoother        GaussSeidel;
    }

    "(U|k|epsilon|omega|f|v2)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent      yes;

    residualControl
    {
        p               1e-2;
        U               1e-3;
        "(k|epsilon|omega|f|v2)" 1e-3;
    }
}

relaxationFactors
{
    equations
    {
        U               0.9; // 0.9 is more stable but 0.95 more convergent
        ".*"            0.9; // 0.9 is more stable but 0.95 more convergent
    }
}

system/blockMeshDict

// Note: this file is a Copy of $FOAM_TUTORIALS/resources/blockMesh/pitzDaily

units [mm];

vertices
(
    (-20.6 0 -0.5)
    (-20.6 25.4 -0.5)
    (0 -25.4 -0.5)
    (0 0 -0.5)
    (0 25.4 -0.5)
    (206 -25.4 -0.5)
    (206 0 -0.5)
    (206 25.4 -0.5)
    (290 -16.6 -0.5)
    (290 0 -0.5)
    (290 16.6 -0.5)

    (-20.6 0 0.5)
    (-20.6 25.4 0.5)
    (0 -25.4 0.5)
    (0 0 0.5)
    (0 25.4 0.5)
    (206 -25.4 0.5)
    (206 0 0.5)
    (206 25.4 0.5)
    (290 -16.6 0.5)
    (290 0 0.5)
    (290 16.6 0.5)
);

negY
(
    (2 4 1)
    (1 3 0.3)
);

posY
(
    (1 4 2)
    (2 3 4)
    (2 4 0.25)
);

posYR
(
    (2 1 1)
    (1 1 0.25)
);


blocks
(
    hex (0 3 4 1 11 14 15 12)
    (18 30 1)
    simpleGrading (0.5 $posY 1)

    hex (2 5 6 3 13 16 17 14)
    (180 27 1)
    edgeGrading (4 4 4 4 $negY 1 1 $negY 1 1 1 1)

    hex (3 6 7 4 14 17 18 15)
    (180 30 1)
    edgeGrading (4 4 4 4 $posY $posYR $posYR $posY 1 1 1 1)

    hex (5 8 9 6 16 19 20 17)
    (25 27 1)
    simpleGrading (2.5 1 1)

    hex (6 9 10 7 17 20 21 18)
    (25 30 1)
    simpleGrading (2.5 $posYR 1)
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 1 12 11)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (8 9 20 19)
            (9 10 21 20)
        );
    }
    upperWall
    {
        type wall;
        faces
        (
            (1 4 15 12)
            (4 7 18 15)
            (7 10 21 18)
        );
    }
    lowerWall
    {
        type wall;
        faces
        (
            (0 3 14 11)
            (3 2 13 14)
            (2 5 16 13)
            (5 8 19 16)
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (
            (0 3 4 1)
            (2 5 6 3)
            (3 6 7 4)
            (5 8 9 6)
            (6 9 10 7)
            (11 14 15 12)
            (13 16 17 14)
            (14 17 18 15)
            (16 19 20 17)
            (17 20 21 18)
        );
    }
);

Useful details: kEpsilon is the active model (the comment lists the models tested with the same case); residualControl would stop the run before the 2000 iterations of endTime; consistent yes activates the SIMPLEC variant (next box); wallDist { method meshWave; } serves models that use the wall distance; the streamlinesLine function writes streamlines at every write time.

The incompressibleFluid module builds the momentum equation in momentumPredictor.C, relaxes it and solves it with the pressure gradient of the previous iteration:

tUEqn =
(
    fvm::ddt(U) + fvm::div(phi, U)
  + MRF.DDt(U)
  + momentumTransport->divDevSigma(U)
 ==
    fvModels().source(U)
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();
...
if (pimple.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p));
}

In correctPressure.C, with consistent yes, the pressure-equation coefficient is not 1/A but the “consistent” one of J. P. Van Doormaal, G. D. Raithby, “Enhancements of the SIMPLE method for predicting incompressible fluid flows”, Numerical Heat Transfer 7, 147–163, 1984, doi:10.1080/01495728408961817:

tmp<volScalarField> rAtU(rAU);

if (pimple.consistent())
{
    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
    phiHbyA +=
        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
    HbyA -= (rAU - rAtU())*fvc::grad(p);
}
...
fvScalarMatrix pEqn
(
    fvm::laplacian(rAtU(), p)
 ...
);
...
// Explicitly relax pressure for momentum corrector
p.relax();

U = HbyA - rAtU*fvc::grad(p);

• With SIMPLEC the velocity correction accounts for the coefficients of neighbouring cells (H1): this is why the tutorial does not under-relax pressure (no fields { p … } entry, so p.relax() has no effect) and uses 0.9 for all equations, whereas classic SIMPLE by S. V. Patankar, D. B. Spalding (International Journal of Heat and Mass Transfer 15, 1787–1806, 1972, doi:10.1016/0017-9310(72)90054-3) typically needs 0.3 on p and 0.7 on U.
• Interpolating the fluxes phiHbyA from cell values follows the idea of C. M. Rhie, W. L. Chow (AIAA Journal 21, 1525–1532, 1983, doi:10.2514/3.8284), which avoids checkerboard pressure oscillations on co-located grids.
• The turbulence model (kEpsilon, B. E. Launder, D. B. Spalding, Computer Methods in Applied Mechanics and Engineering 3, 269–289, 1974, doi:10.1016/0045-7825(74)90029-2) is solved after pressure: for every iteration the log shows Ux, Uy, p, the continuity line, epsilon and k.

• Iterative error. It is the difference between the solution at the current iteration and the one the discrete system would have at convergence. It is not measured by the residual of a single linear solver but by watching the quantities of interest: if the pressure drop and the reattachment point still change by 0.1 % every 100 iterations, the iterative error is of that order (section 15, observables.py).
• Discretisation error. It depends on mesh and schemes; once iteratively converged it is estimated with at least three meshes in a constant ratio, the observed order and the GCI: P. J. Roache, Journal of Fluids Engineering 116, 405–413, 1994, doi:10.1115/1.2910291; I. B. Celik et al., Journal of Fluids Engineering 130, 078001, 2008, doi:10.1115/1.2960953 (Guide 4). Comparing schemes on a single mesh shows how much the discretisation error depends on the scheme, not which scheme is “right”.
• Rule of thumb. The iterative error must be at least an order of magnitude smaller than the differences one wants to attribute to the schemes: otherwise the comparison only measures how far each case is from convergence after 600 iterations.
• Reference texts: J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, 4th ed., Springer, 2020, doi:10.1007/978-3-319-99693-6; F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, Springer, 2016, doi:10.1007/978-3-319-16874-6 (with OpenFOAM examples).

bounded. In boundedConvectionScheme.C the bounded variant subtracts from the convection term the product of the flux divergence and the field:

return
    scheme_().fvmDiv(faceFlux, vf)
  - fvm::Sp(fvc::surfaceIntegrate(faceFlux), vf);

That is, it discretises ∇·(φU) − (∇·φ)U. At convergence ∇·φ = 0 and the term vanishes; during the iterations of a steady case, when continuity is not yet satisfied, it prevents the mass error from acting as a spurious momentum source. It is a stability device for the SIMPLE loop, independent of the interpolation scheme that follows.

linearUpwind. In linearUpwind.C the face value is that of the upstream cell plus an explicit gradient correction (here grad(U), computed with Gauss linear):

forAll(faceFlux, facei)
{
    const label celli =
        (faceFlux[facei] > 0) ? owner[facei] : neighbour[facei];

    sfCorr[facei] = (Cf[facei] - C[celli]) & gradVf[celli];
}

• It is the second-order extrapolation of R. F. Warming, R. M. Beam, “Upwind second-order difference schemes and applications in aerodynamic flows”, AIAA Journal 14, 1241–1249, 1976, doi:10.2514/3.61457. The upwind part is implicit, the correction explicit: it is unlimited, so it can produce small over/undershoots where gradients change abruptly, such as in the shear layer after the step.
• upwind. The upstream value only: bounded and robust, but first-order, with a truncation error equivalent to a diffusivity |U|Δx/2 along the flow; in this case it is 1 to 30 times the inlet eddy viscosity (section 1). It typically shortens recirculation zones.
• limitedLinearV 1. Linear interpolation limited with the limitedLinear TVD limiter (coefficient 1 = strongest limiting, towards upwind); the “V” denotes the vector version, which computes the limiter once in the direction of maximum change of the vector instead of component by component, preserving the velocity direction. Limiter theory: P. K. Sweby, “High resolution schemes using flux limiters for hyperbolic conservation laws”, SIAM Journal on Numerical Analysis 21, 995–1011, 1984, doi:10.1137/0721062.

• The reattachment point on the lower wall is the observable most sensitive to numerical diffusion: observables.py derives it from the wall shear stress at every write time (section 15). With upwind a shorter recirculation is expected.
• The pressure drop between inlet and outlet is an integral quantity, less sensitive; Ux profiles at sections x = 2h, 4h, 6h (Plot Over Line in ParaView) show where the schemes differ.
• Only the U scheme changes: k and ε remain bounded Gauss limitedLinear 1 in all cases, as in the tutorial.

The residual printed in the log is normalised in lduMatrixSolver.C: the sum of the magnitudes of b − Ax is divided by a factor depending on the solution itself, so the value does not change with the scale of the quantities:

// --- Calculate A dot reference value of psi
matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);

tmpField *= gAverage(psi, matrix_.lduMesh_.comm());

return
    gSum
    (
        (mag(Apsi - tmpField) + mag(source - tmpField))(),
        matrix_.lduMesh_.comm()
    )
  + solverPerformance::small_;

The linear solver stopping test, in SolverPerformance.C, is “final residual < tolerance” or “final residual < relTol × initial residual”:

if
(
    finalResidual_ < Tolerance
 || (
        RelTolerance
      > small_*pTraits<Type>::one
     && finalResidual_ < cmptMultiply(RelTolerance, initialResidual_)
    )
)
{
    converged_ = true;
}

• With relTol 0.1 every p solution reduces the residual by a factor of 10 and stops: in the first iterations a few GAMG cycles suffice, towards convergence the initial residual is already small and the work per iteration changes. This is why the number of linear iterations (No Iterations) must be recorded: residual_history.py sums them per case.
• pcorr (relTol 0) is used only for the initial flux correction; GAMG is an algebraic multigrid (K. Stüben, “A review of algebraic multigrid”, Journal of Computational and Applied Mathematics 128, 281–309, 2001, doi:10.1016/S0377-0427(00)00516-1) with a GaussSeidel smoother; for U, k and ε the tutorial uses smoothSolver with symGaussSeidel, suitable for asymmetric, diagonally dominant matrices.

residualControl. In singleRegionConvergenceControl.C the loop stops when, for all listed fields, the initial residual of the first solution in the iteration is below the threshold; for vectors the largest component counts (cmptMax in convergenceControlTemplates.C):

getInitialResiduals(mesh_, fieldName, 0, residual, residual);
...
const bool absCheck = residual < residualControl_[fieldi].absTol;
cs.satisfied = cs.satisfied && absCheck;

// convergenceControlTemplates.C
r0 = cmptMax(sp[0].initialResidual());

The script reads log.foamRun of each case iteration by iteration and computes: the initial residual of the first solution of each field (the one used by residualControl), the iteration at which the tutorial criteria (p 10⁻², U 10⁻³, k and ε 10⁻³, changeable with --criteria) would have been met, total and mean linear iterations of p, the final continuity error and the mean cost per iteration from successive ExecutionTime values (excluding the first 10). It writes residual_history.csv and, with --plot, the plot of the Ux and p residuals.

#!/usr/bin/env python3
"""Residual and cost histories of the Guide 19 campaign from log.foamRun.

For every case listed in campaign_manifest.csv (or given on the command
line) the script parses log.foamRun step by step and reports:

  * the initial residual of the FIRST solution of each field in each step
    (the quantity used by SIMPLE/residualControl; for U the largest of the
    solved components, as cmptMax does in OpenFOAM 14);
  * the step at which the tutorial residualControl criteria
    (p 1e-2, U 1e-3, k/epsilon/omega/f/v2 1e-3, or --criteria) would have
    stopped the original tutorial;
  * linear-solver work: total and mean "No Iterations" of p and U;
  * the local continuity error at the last step;
  * the mean cost per step from consecutive ExecutionTime values,
    excluding the first 10 steps (start-up).

Outputs: a summary on screen, residual_history.csv (one row per case and
step) and, with --plot, residual_history.png (needs matplotlib). Existing
output files are never overwritten.

Usage:
  python3 residual_history.py CAMPAIGN_DIR [--cases reference upwind ...]
                              [--criteria "p=1e-2,U=1e-3,k=1e-3,epsilon=1e-3"] [--plot]
"""
import argparse, csv, math, os, re, sys

STEP = re.compile(r'^Time = ([0-9.eE+-]+)s?\s*$')
SOLVE = re.compile(r'^(\w+):\s+Solving for (\w+), Initial residual = ([0-9.eE+-]+), '
                   r'Final residual = ([0-9.eE+-]+), No Iterations (\d+)')
CONT = re.compile(r'time step continuity errors : sum local = ([0-9.eE+-]+), global = ([0-9.eE+-]+)')
TIMES = re.compile(r'^ExecutionTime = ([0-9.eE+-]+) s\s+ClockTime = ([0-9.eE+-]+) s')


def parse(log):
    steps = []
    cur = None
    for line in open(log, encoding='utf-8', errors='replace'):
        line = line.rstrip('\n')
        m = STEP.match(line)
        if m:
            cur = dict(step=float(m.group(1)), first={}, iters={}, cont=None, cpu=None, clock=None)
            steps.append(cur)
            continue
        if cur is None:
            continue
        m = SOLVE.match(line)
        if m:
            f = m.group(2)
            cur['first'].setdefault(f, float(m.group(3)))
            cur['iters'][f] = cur['iters'].get(f, 0) + int(m.group(5))
            continue
        m = CONT.search(line)
        if m:
            cur['cont'] = float(m.group(1))
            continue
        m = TIMES.match(line)
        if m:
            cur['cpu'], cur['clock'] = float(m.group(1)), float(m.group(2))
    return [s for s in steps if s['cpu'] is not None]


def field_residual(step, name):
    """Residual of a field as residualControl sees it (vector -> max component)."""
    if name in step['first']:
        return step['first'][name]
    comps = [v for k, v in step['first'].items() if k.startswith(name) and k[len(name):] in ('x', 'y', 'z')]
    return max(comps) if comps else None


def satisfied(step, crit):
    """residualControl test: every controlled field that is solved is below its tolerance."""
    checked = [(field_residual(step, f), tol) for f, tol in crit.items()]
    checked = [(r, tol) for r, tol in checked if r is not None]     # fields not solved (e.g. omega with kEpsilon) are ignored
    return bool(checked) and all(r < tol for r, tol in checked)


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('campaign')
    ap.add_argument('--cases', nargs='*')
    ap.add_argument('--criteria', default='p=1e-2,U=1e-3,k=1e-3,epsilon=1e-3,omega=1e-3,f=1e-3,v2=1e-3')
    ap.add_argument('--plot', action='store_true')
    a = ap.parse_args()
    base = os.path.abspath(os.path.expanduser(a.campaign))
    if a.cases:
        cases = a.cases
    else:
        with open(os.path.join(base, 'campaign_manifest.csv'), encoding='utf-8') as f:
            cases = [r['case'] for r in csv.DictReader(f)]
    crit = {k: float(v) for k, v in (c.split('=') for c in a.criteria.split(','))}
    out_csv = os.path.join(base, 'residual_history.csv')
    out_png = os.path.join(base, 'residual_history.png')
    for p in (out_csv,) + ((out_png,) if a.plot else ()):
        if os.path.exists(p):
            sys.exit('refusing to overwrite ' + p)
    data = {}
    for c in cases:
        log = os.path.join(base, c, 'log.foamRun')
        if not os.path.isfile(log):
            print(f'{c:12s} no log.foamRun')
            continue
        st = parse(log)
        if not st:
            print(f'{c:12s} log without complete steps')
            continue
        data[c] = st
    fields = ['Ux', 'Uy', 'p', 'k', 'epsilon']
    print(f'{"case":12s} {"steps":>5s} ' + ' '.join(f'{f:>9s}' for f in fields) +
          f' {"stop@":>6s} {"p iters":>8s} {"p it/step":>9s} {"cont.loc":>9s} {"ms/step":>8s}')
    rows = []
    for c, st in data.items():
        last = st[-1]
        stop = next((int(s['step']) for s in st if satisfied(s, crit)), None)
        pit = sum(s['iters'].get('p', 0) for s in st)
        cpu = [s['cpu'] for s in st]
        k0 = min(10, len(cpu) - 2)
        ms = 1e3 * (cpu[-1] - cpu[k0]) / (len(cpu) - 1 - k0) if len(cpu) > k0 + 1 else float('nan')
        print(f'{c:12s} {int(last["step"]):5d} ' + ' '.join(f'{last["first"].get(f, float("nan")):9.2e}' for f in fields) +
              f' {str(stop) if stop else "-":>6s} {pit:8d} {pit / len(st):9.1f} {last["cont"] if last["cont"] is not None else float("nan"):9.2e} {ms:8.2f}')
        for s in st:
            rows.append([c, int(s['step'])] + [s['first'].get(f, '') for f in fields] +
                        [s['iters'].get('p', ''), s['iters'].get('Ux', ''), s['cont'], s['cpu'], s['clock']])
    print('stop@ = first step at which all residualControl criteria (' + a.criteria + ') hold, i.e. where the tutorial would stop;')
    print('        "-" = never within the run. Residuals are initial residuals of the first solution in each step.')
    with open(out_csv, 'x', newline='', encoding='utf-8') as f:
        w = csv.writer(f)
        w.writerow(['case', 'step'] + ['r0_' + x for x in fields] + ['p_iterations', 'Ux_iterations', 'continuity_local', 'ExecutionTime_s', 'ClockTime_s'])
        w.writerows(rows)
    print('Wrote', out_csv)
    if a.plot and data:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        fig, axs = plt.subplots(1, 2, figsize=(11, 4.2), dpi=150, sharey=True)
        for ax, f in zip(axs, ('Ux', 'p')):
            for c, st in data.items():
                if c.startswith('parallel'):
                    continue
                ax.semilogy([s['step'] for s in st], [s['first'].get(f, math.nan) for s in st], lw=1.2, label=c)
            ax.axhline(crit['U' if f.startswith('U') else f], color='gray', ls=':', lw=1)
            ax.set_title(f'initial residual of {f}', fontsize=10)
            ax.set_xlabel('SIMPLE iteration')
            ax.grid(alpha=0.3, which='both')
        axs[0].legend(fontsize=8)
        fig.tight_layout()
        fig.savefig(out_png)
        print('Wrote', out_png)


if __name__ == '__main__':
    main()

Test on the synthetic logs written by the test generator in the section 15 box (exponential models, not OpenFOAM results): the line format is that of the OpenFOAM 14 foamRun log (Time = 600s, Solving for Ux, Initial residual = …, ExecutionTime = … s ClockTime = … s).

$ python3 residual_history.py ~/OpenFOAM/guide19_run01 --plot
case         steps        Ux        Uy         p         k   epsilon  stop@  p iters p it/step  cont.loc  ms/step
reference      600  3.81e-05  4.31e-05  3.27e-04  2.82e-05  1.46e-05    417     5575       9.3  5.48e-04    19.40
upwind         600  2.50e-06  3.55e-06  1.69e-05  1.37e-06  1.19e-06    289     5575       9.3  2.03e-05    19.40
limited        600  9.75e-06  1.16e-05  1.04e-04  6.85e-06  3.92e-06    366     5575       9.3  1.68e-04    19.40
relax_06       600  3.84e-03  4.27e-03  1.17e-02  2.90e-03  1.45e-03      -     5575       9.3  1.98e-02    19.40
parallel_2     600  3.81e-05  4.31e-05  3.27e-04  2.82e-05  1.46e-05    417     5575       9.3  5.48e-04    11.45
parallel_4     600  3.81e-05  4.31e-05  3.27e-04  2.82e-05  1.46e-05    417     5575       9.3  5.48e-04     7.66
stop@ = first step at which all residualControl criteria (p=1e-2,U=1e-3,k=1e-3,epsilon=1e-3,omega=1e-3,f=1e-3,v2=1e-3) hold, i.e. where the tutorial would stop;
        "-" = never within the run. Residuals are initial residuals of the first solution in each step.
Wrote /home/user/OpenFOAM/guide19_run01/residual_history.csv
Wrote /home/user/OpenFOAM/guide19_run01/residual_history.png

For equations (relaxationFactors/equations) OpenFOAM applies implicit relaxation in fvMatrix.C: it makes the matrix diagonally dominant, divides the diagonal by α and adds the difference to the source, so that the relaxed solution remains consistent with the original equation at convergence:

// Ensure the matrix is diagonally dominant...
forAll(D, celli)
{
    D[celli] = max(mag(D[celli]), sumOff[celli]);
}

// ... then relax
D /= alpha;
...
// Finally add the relaxation contribution to the source.
S += (D - D0)*psi_.primitiveField();

• At convergence U no longer changes between two iterations and the added term vanishes: implicit relaxation changes the path, not the final solution. A difference between reference and relax_06 after 600 iterations therefore measures only the different distance from convergence.
• With SIMPLEC the consistent coefficient rAtU is built from the relaxed matrix: lowering α from 0.9 to 0.6 also slows the pressure correction. A clearly higher number of iterations to reach the same residuals is expected; residual_history.py measures it with the “stop@” column.
• The tutorial comment (“0.9 is more stable but 0.95 more convergent”) points the other way: with SIMPLEC one can go up towards 0.95 if the case is stable.

The script rebuilds the cell sizes with the blockMesh rules: for each block edge it reads the grading (simpleGrading or the 12 entries of edgeGrading, including multi-section “length fraction, cell fraction, ratio” gradings), computes the geometric progression of the cells and reports minima, maxima and the largest aspect ratio per block; then the physical inlet quantities and the upwind numerical diffusion. Python standard library only.

#!/usr/bin/env python3
"""A-priori numbers for the pitzDailySteady campaign (OpenFOAM 14).

Reads the blockMesh dictionary (vertices, macros, simpleGrading/edgeGrading
with multi-grading sections), 0/U, 0/k, 0/epsilon and
constant/physicalProperties of a case and prints:

  * cells per block and in total, cell sizes along every block edge
    (blockMesh grading rules), smallest/largest cells and aspect ratios;
  * step height, Reynolds number, inlet turbulence intensity, length scale
    and eddy viscosity from the inlet k and epsilon;
  * the numerical diffusivity of first-order upwind, |U| dx / 2, compared
    with the laminar and inlet eddy viscosities;
  * cells per MPI process for 1, 2, 4 and 8 processes.

Usage:
  python3 case_estimates.py CASE [--blockmesh FILE]
  (default FILE: $FOAM_TUTORIALS/resources/blockMesh/pitzDaily if it exists,
   otherwise CASE/system/blockMeshDict)
Only the Python standard library is needed.
"""
import argparse, math, os, re, sys

UNITS = {'m': 1.0, 'cm': 1e-2, 'mm': 1e-3, 'um': 1e-6}


def strip(t):
    return re.sub(r'//[^\n]*', ' ', re.sub(r'/\*.*?\*/', ' ', t, flags=re.S))


def balanced(text, start, o='(', c=')'):
    i = text.index(o, start)
    depth = 0
    for j in range(i, len(text)):
        if text[j] == o:
            depth += 1
        elif text[j] == c:
            depth -= 1
            if depth == 0:
                return text[i + 1:j], j + 1
    raise ValueError('unbalanced ' + o)


def parse_list(s):
    """Nested OpenFOAM list of numbers -> nested Python lists."""
    toks = re.findall(r'\(|\)|[-+0-9.eE]+', s)
    stack = [[]]
    for t in toks:
        if t == '(':
            stack.append([])
        elif t == ')':
            v = stack.pop()
            stack[-1].append(v)
        else:
            stack[-1].append(float(t))
    return stack[0]


def read_blockmesh(path):
    t = strip(open(path, encoding='utf-8').read())
    m = re.search(r'units\s+\[(\w+)\]\s*;', t)
    scale = UNITS[m.group(1)] if m else float((re.search(r'convertToMeters\s+([-+0-9.eE]+)', t) or [0, 1])[1])
    macros = {}
    for mm in re.finditer(r'(?m)^\s*(\w+)\s*\n?\s*\(', t):
        name = mm.group(1)
        if name in ('vertices', 'blocks', 'boundary', 'edges', 'faces', 'defaultPatch', 'mergePatchPairs'):
            continue
        body, _ = balanced(t, mm.end() - 1)
        macros[name] = '(' + body + ')'
    vb, _ = balanced(t, t.index('vertices'))
    V = [[scale * float(c) for c in v] for v in parse_list(vb)]
    bb, _ = balanced(t, t.index('blocks'))
    for k, v in macros.items():
        bb = bb.replace('$' + k, v)
    blocks = []
    for hm in re.finditer(r'hex\s*\(([^)]*)\)\s*\(([^)]*)\)\s*(simpleGrading|edgeGrading)', bb):
        verts = [int(x) for x in hm.group(1).split()]
        n = [int(float(x)) for x in hm.group(2).split()]
        g, _ = balanced(bb, hm.end())
        blocks.append(dict(verts=verts, n=n, kind=hm.group(3), grading=parse_list(g)))
    return V, blocks


def sections(spec):
    """Grading spec of one edge: number or list of (lengthFrac, cellFrac, ratio)."""
    if isinstance(spec, (int, float)):
        return [(1.0, 1.0, float(spec))]
    tl = sum(s[0] for s in spec)
    tc = sum(s[1] for s in spec)
    return [(s[0] / tl, s[1] / tc, s[2]) for s in spec]


def edge_cells(L, n, spec):
    secs = sections(spec)
    counts = [max(1, int(round(cf * n))) for _, cf, _ in secs]
    counts[-1] = n - sum(counts[:-1])
    out = []
    for (lf, _, r), nc in zip(secs, counts):
        ls = lf * L
        if nc == 1 or abs(r - 1) < 1e-12:
            out += [ls / nc] * nc
        else:
            q = r ** (1.0 / (nc - 1))
            d1 = ls * (1 - q) / (1 - q ** nc)
            out += [d1 * q ** i for i in range(nc)]
    return out


EDGES = [(0, 1), (3, 2), (7, 6), (4, 5), (0, 3), (1, 2), (5, 6), (4, 7), (0, 4), (1, 5), (2, 6), (3, 7)]


def block_sizes(V, b):
    res = []
    for d in range(3):
        sizes = []
        for e in range(4):
            i, j = EDGES[4 * d + e]
            p, q = V[b['verts'][i]], V[b['verts'][j]]
            L = math.dist(p, q)
            spec = b['grading'][d] if b['kind'] == 'simpleGrading' else b['grading'][4 * d + e]
            sizes.append(edge_cells(L, b['n'][d], spec))
        res.append(sizes)
    return res


def field_value(case, name, patch):
    t = strip(open(os.path.join(case, '0', name), encoding='utf-8').read())
    m = re.search(re.escape(patch) + r'\s*\{[^}]*?value\s+uniform\s+(\([^)]*\)|[-+0-9.eE]+)', t, re.S)
    v = m.group(1)
    if v.startswith('$'):
        v = re.search(r'internalField\s+uniform\s+([-+0-9.eE]+)', t).group(1)
    return [float(x) for x in v.strip('()').split()] if v.startswith('(') else float(v)


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('case')
    ap.add_argument('--blockmesh')
    a = ap.parse_args()
    bm = a.blockmesh
    if not bm:
        tut = os.path.join(os.environ.get('FOAM_TUTORIALS', ''), 'resources', 'blockMesh', 'pitzDaily')
        bm = tut if os.environ.get('FOAM_TUTORIALS') and os.path.isfile(tut) else os.path.join(a.case, 'system', 'blockMeshDict')
    V, blocks = read_blockmesh(bm)
    print(f'blockMesh dictionary: {os.path.basename(bm)}  ({len(V)} vertices, {len(blocks)} blocks)')
    tot = 0
    dmin, dmax, armax = 1e9, 0, 0
    print('block  cells (nx x ny)   dx min-max [mm]     dy min-max [mm]     max aspect dx/dy or dy/dx')
    for k, b in enumerate(blocks):
        nc = b['n'][0] * b['n'][1] * b['n'][2]
        tot += nc
        sx, sy, _ = block_sizes(V, b)
        allx = [d for s in sx for d in s]
        ally = [d for s in sy for d in s]
        # aspect ratio estimated per edge pair: largest dx / smallest dy and vice versa
        ar = max(max(allx) / min(ally), max(ally) / min(allx))
        armax = max(armax, ar)
        dmin = min(dmin, min(allx), min(ally))
        dmax = max(dmax, max(allx), max(ally))
        print(f'{k:5d}  {nc:6d} ({b["n"][0]:3d} x {b["n"][1]:3d})   {1e3 * min(allx):6.3f} - {1e3 * max(allx):6.3f}     '
              f'{1e3 * min(ally):6.3f} - {1e3 * max(ally):6.3f}     {ar:8.1f}')
    print(f'total cells: {tot}   smallest cell edge {1e3 * dmin:.3f} mm, largest {1e3 * dmax:.3f} mm, '
          f'bounding aspect ratio <= {armax:.1f}')
    ys = [v[1] for v in V]
    xs = [v[0] for v in V]
    h = -min(ys)                      # step height: lower wall below y = 0
    Hin = max(ys)
    U = field_value(a.case, 'U', 'inlet')[0]
    k = field_value(a.case, 'k', 'inlet')
    eps = field_value(a.case, 'epsilon', 'inlet')
    t = strip(open(os.path.join(a.case, 'constant', 'physicalProperties'), encoding='utf-8').read())
    nu = float(re.search(r'nu\s+(?:\[[^\]]*\]\s*)?([-+0-9.eE]+)', t).group(1))
    Cmu = 0.09
    print(f'\ngeometry: inlet channel height {1e3 * Hin:.1f} mm, step height h = {1e3 * h:.1f} mm, '
          f'domain x from {1e3 * min(xs):.1f} to {1e3 * max(xs):.1f} mm ({(max(xs)) / h:.1f} h downstream of the step)')
    print(f'inlet velocity {U:g} m/s, nu = {nu:g} m2/s -> Re_h = U h / nu = {U * h / nu:.0f}, Re_H(inlet) = {U * Hin / nu:.0f}')
    I = math.sqrt(2 * k / 3) / U
    nut = Cmu * k * k / eps
    lm = Cmu ** 0.75 * k ** 1.5 / eps
    print(f'inlet turbulence: k = {k:g} m2/s2 -> intensity {100 * I:.1f} %; epsilon = {eps:g} m2/s3 -> '
          f'length scale Cmu^0.75 k^1.5/eps = {1e3 * lm:.2f} mm ({lm / Hin:.2f} H), nut = Cmu k^2/eps = {nut:.3e} m2/s ({nut / nu:.0f} nu)')
    print(f'consistency: omega = eps/(Cmu k) = {eps / (Cmu * k):.2f} 1/s (tutorial 0/omega uses 440.15)')
    print('\nfirst-order upwind numerical diffusivity |U| dx / 2 (flow aligned with the grid):')
    for dx in sorted({dmin, 0.5e-3, 1e-3, 2e-3, dmax}):
        nn = U * dx / 2
        print(f'  dx = {1e3 * dx:6.3f} mm  ->  {nn:.3e} m2/s = {nn / nu:7.0f} nu = {nn / nut:5.2f} x inlet nut')
    print('\ncells per MPI process (decomposition balance ignored):')
    for n in (1, 2, 4, 8):
        print(f'  {n} process(es): {tot / n:8.0f} cells each')


if __name__ == '__main__':
    main()

$ python3 case_estimates.py ~/OpenFOAM/guide19_run01/reference
blockMesh dictionary: pitzDaily  (22 vertices, 5 blocks)
block  cells (nx x ny)   dx min-max [mm]     dy min-max [mm]     max aspect dx/dy or dy/dx
    0     540 ( 18 x  30)    0.791 -  1.583      0.319 -  2.296          5.0
    1    4860 (180 x  27)    0.528 -  2.114      0.360 -  1.201          5.9
    2    5400 (180 x  30)    0.528 -  2.114      0.319 -  2.296          6.6
    3     675 ( 25 x  27)    2.047 -  5.145      0.615 -  0.941          8.4
    4     750 ( 25 x  30)    2.047 -  5.145      0.208 -  2.296         24.7
total cells: 12225   smallest cell edge 0.208 mm, largest 5.145 mm, bounding aspect ratio <= 24.7

geometry: inlet channel height 25.4 mm, step height h = 25.4 mm, domain x from -20.6 to 290.0 mm (11.4 h downstream of the step)
inlet velocity 10 m/s, nu = 1e-05 m2/s -> Re_h = U h / nu = 25400, Re_H(inlet) = 25400
inlet turbulence: k = 0.375 m2/s2 -> intensity 5.0 %; epsilon = 14.855 m2/s3 -> length scale Cmu^0.75 k^1.5/eps = 2.54 mm (0.10 H), nut = Cmu k^2/eps = 8.520e-04 m2/s (85 nu)
consistency: omega = eps/(Cmu k) = 440.15 1/s (tutorial 0/omega uses 440.15)

first-order upwind numerical diffusivity |U| dx / 2 (flow aligned with the grid):
  dx =  0.208 mm  ->  1.042e-03 m2/s =     104 nu =  1.22 x inlet nut
  dx =  0.500 mm  ->  2.500e-03 m2/s =     250 nu =  2.93 x inlet nut
  dx =  1.000 mm  ->  5.000e-03 m2/s =     500 nu =  5.87 x inlet nut
  dx =  2.000 mm  ->  1.000e-02 m2/s =    1000 nu = 11.74 x inlet nut
  dx =  5.145 mm  ->  2.572e-02 m2/s =    2572 nu = 30.19 x inlet nut

cells per MPI process (decomposition balance ignored):
  1 process(es):    12225 cells each
  2 process(es):     6112 cells each
  4 process(es):     3056 cells each
  8 process(es):     1528 cells each

• The total, 12,225 cells, matches the sum of the dictionary blocks (540 + 4860 + 5400 + 675 + 750): it is the number to find in log.checkMesh (“cells: 12225”).
• The reported aspect ratio is an upper bound (longest cell over thinnest cell of the block); the actual one, cell by cell, is in the “Max aspect ratio” line of checkMesh. Block 4, tapered and with the $posYR grading, has the most elongated cells.
• From log.checkMesh record in the protocol: number of cells, “Max non-orthogonality”, “Max skewness”, “Max aspect ratio” and any “Failed … mesh checks” line. The faces of the tapered blocks are not orthogonal: the Gauss linear corrected correction of the Laplacians compensates non-orthogonality, with nNonOrthogonalCorrectors 0.

run_case.sh always builds the mesh from the resources dictionary (blockMesh -dict $FOAM_TUTORIALS/resources/blockMesh/pitzDaily), so a refined mesh must be prepared in a new copy and launched by hand, using the case system/blockMeshDict:

cp -r ~/OpenFOAM/guide19_run01/reference ~/OpenFOAM/guide19_mesh2x     # copy BEFORE running reference
cd ~/OpenFOAM/guide19_mesh2x
sed -i 's/(18 30 1)/(36 60 1)/; s/(180 27 1)/(360 54 1)/; s/(180 30 1)/(360 60 1)/; s/(25 27 1)/(50 54 1)/; s/(25 30 1)/(50 60 1)/' system/blockMeshDict
grep -n "hex" -A1 system/blockMeshDict        # check the five blocks
blockMesh > log.blockMesh 2>&1
checkMesh > log.checkMesh 2>&1                 # expect cells: 48900
foamRun > log.foamRun 2>&1

The grading ratio is defined between the last and the first cell of the section: doubling the cells with the same ratio halves the end cells and keeps the progression geometrically similar. With 4 times the cells more SIMPLE iterations are needed to converge: 600 iterations are not enough in general, and the comparison with the base mesh must be made at iterative convergence (section 3).

• In OpenFOAM 14 foamRun prints Time = followed by the time name with its unit (foamRun.C, runTime.userTimeName()): the log reads “Time = 600s”. The scripts’ regular expressions accept both “600” and “600s”.
• With writeInterval 100 and 600 iterations the folders 100…600 are written: they are iterations, not seconds. foamListTimes lists them; in parallel cases they are inside processorN/ until reconstructPar -latestTime rebuilds the last one.
• The original tutorial, with residualControl active, stops by itself when the residuals fall below the thresholds and writes that last iteration; in the campaign copies the control is removed and all variants run exactly the same number of iterations.

• scotch partitions the cell graph minimising the boundary faces between processors, with no geometric coefficients (F. Pellegrini, J. Roman, “Scotch: a software package for static mapping by dual recursive bipartitioning of process and architecture graphs”, HPCN Europe 1996, LNCS 1067, 493–498, doi:10.1007/3-540-61142-8_588). The OpenFOAM 14 manual (section 3.4) also describes the geometric simple and hierarchical methods.
• In log.decomposePar read, for each processor, the number of cells and of processor faces: their sum relative to the internal faces indicates how much communication each iteration needs. With about 3000 cells per process (4 processes) a non-negligible fraction of the cells borders another subdomain.
• decomposePar -cellProc writes the cellProc field (processor number of each cell), to be viewed in ParaView to see the subdomains.
• Use physical cores only: lscpu shows “Core(s) per socket” and “Thread(s) per core”. With Open MPI, --bind-to core pins each process to a core and reduces the variability between repetitions; “There are not enough slots available in the system” means more processes are requested than available cores (--oversubscribe forces it, but the measured time is no longer a meaningful speedup).
• The parallel result must match the serial one within the solver tolerance: the initial residuals in the log may differ slightly because GAMG and the smoothers work per subdomain; pressure drop and reattachment point must stay equal to the chosen significant digits.

ParaView reads decomposed cases directly: in the OpenFOAM reader choose “Decomposed Case” in Case Type. This is useful to check the subdomains (with cellProc) and the intermediate times, which reconstructPar -latestTime does not rebuild. To rebuild everything: reconstructPar with no options, in a copy after the benchmark, not in the folder used for timing.

The line printed at the end of every iteration comes from foamRun.C:

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
    << "  ClockTime = " << runTime.elapsedClockTime() << " s"
    << nl << endl;

ExecutionTime is the CPU time of the process: cpuTime.C uses the C library clock(). Info writes only from the master process, so in a parallel run it is the CPU time of rank 0, not the sum over processes (the guide has been corrected on this point). With MPI it can approach the elapsed time, because busy-waiting for messages consumes CPU.

// src/OSspecific/POSIX/cpuTime/cpuTime.C
double Foam::cpuTime::elapsedCpuTime() const
{
    curTime_ = clock();
    return timeDifference(startTime_, curTime_);   // (cur - prev)/CLOCKS_PER_SEC
}

// src/OpenFOAM/global/clock/clock.C
time_t Foam::clock::elapsedClockTime() const
{
    newTime_ = getTime();
    return newTime_ - startTime_;
}

ClockTime is a time_t: whole seconds since foamRun started (clock.C). On a 10 s run the rounding error reaches 10 %, and a speedup computed from two integer values can vary a lot from one test to the next. analyze_campaign.py uses exactly this value; for accurate measurements extra/benchmark_repeats.py measures the elapsed time of each phase with time.perf_counter() and repeats the runs.

The laws for reading the numbers. Amdahl’s law, S(n) = 1/(f + (1 − f)/n) with f the non-parallelisable fraction (G. M. Amdahl, “Validity of the single processor approach to achieving large scale computing capabilities”, AFIPS Spring Joint Computer Conference, 483–485, 1967, doi:10.1145/1465482.1465560); scaling with growing workload (J. L. Gustafson, “Reevaluating Amdahl’s law”, Communications of the ACM 31, 532–533, 1988, doi:10.1145/42411.42415): for a fixed-size problem Amdahl applies, for meshes that grow with the processes the “weak scaling” efficiency matters; Karp–Flatt metric, e(n) = (1/S − 1/n)/(1 − 1/n), serial fraction derived from measurements (A. H. Karp, H. P. Flatt, “Measuring parallel processor performance”, Communications of the ACM 33, 539–543, 1990, doi:10.1145/78607.78614): if e grows with n the limit is not serial code but communication cost.

benchmark_repeats.py takes a prepared, never-run case (for example parallel_4), copies it N times into a new folder and repeats the same sequence as run_case.sh, measuring the elapsed time of each phase (blockMesh, checkMesh, decomposePar, solver, reconstructPar). In timing.csv it also records the last ExecutionTime/ClockTime of the log, the number of iterations, the server load average before each run, the number of cores and the CPU model. The source case is not touched; an error stops the repetitions.

#!/usr/bin/env python3
"""Repeated, timed runs of one prepared Guide 19 case (OpenFOAM 14).

The ClockTime printed by foamRun is an integer number of seconds (time_t in
src/OpenFOAM/global/clock), too coarse for short runs, and a single run
does not give the median the guide asks for. This script takes a prepared,
never-run case (for example parallel_4 made by prepare_campaign.py), copies
it N times into a NEW directory and runs each copy with the same sequence
as run_case.sh, timing every phase with a sub-millisecond wall clock:

  blockMesh -dict $FOAM_TUTORIALS/resources/blockMesh/pitzDaily
  checkMesh
  decomposePar                      (only with more than one process)
  foamRun | mpirun -np N foamRun -parallel
  reconstructPar -latestTime        (only with more than one process)

Each copy keeps its own logs. timing.csv records, per repetition, the wall
time of every phase, the last ExecutionTime/ClockTime of the log, the number
of steps, the 1-minute load average before the run and the host data
(cores, CPU model). The source case is never modified.

Usage:
  python3 benchmark_repeats.py CASE --repeats 3 [--out CASE_bench] [--mpirun-args "--bind-to core"]
"""
import argparse, csv, glob, os, re, shlex, shutil, subprocess, sys, time


def host():
    model = ''
    try:
        for line in open('/proc/cpuinfo'):
            if line.startswith('model name'):
                model = line.split(':', 1)[1].strip()
                break
    except OSError:
        pass
    return os.cpu_count(), model


def pristine(case):
    bad = [f for f in ('log.blockMesh', 'log.foamRun', 'RUN_STARTED', 'postProcessing') if os.path.exists(os.path.join(case, f))]
    bad += glob.glob(os.path.join(case, 'processor*'))
    if os.path.isdir(os.path.join(case, 'constant', 'polyMesh')):
        bad.append('constant/polyMesh')
    return bad


def run(cmd, cwd, log):
    t0 = time.perf_counter()
    with open(os.path.join(cwd, log), 'w') as f:
        rc = subprocess.call(cmd, cwd=cwd, stdout=f, stderr=subprocess.STDOUT)
    return rc, time.perf_counter() - t0


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('case')
    ap.add_argument('--repeats', type=int, default=3)
    ap.add_argument('--out')
    ap.add_argument('--mpirun-args', default='', help='extra mpirun options, e.g. "--bind-to core"')
    a = ap.parse_args()
    case = os.path.abspath(os.path.expanduser(a.case))
    out = os.path.abspath(os.path.expanduser(a.out or case.rstrip('/') + '_bench'))
    tut = os.environ.get('FOAM_TUTORIALS')
    if not tut:
        sys.exit('source the OpenFOAM 14 environment first (FOAM_TUTORIALS is not set)')
    for c in ('blockMesh', 'checkMesh', 'foamRun'):
        if not shutil.which(c):
            sys.exit('missing command: ' + c)
    if not os.path.isdir(os.path.join(case, 'system')):
        sys.exit('not an OpenFOAM case: ' + case)
    bad = pristine(case)
    if bad:
        sys.exit('the source case has already been run or meshed (' + ', '.join(os.path.basename(b) for b in bad) + '): use a fresh prepared case')
    if os.path.exists(out):
        sys.exit('refusing to overwrite ' + out)
    if not 1 <= a.repeats <= 50:
        sys.exit('--repeats must be between 1 and 50')
    dp = os.path.join(case, 'system', 'decomposeParDict')
    n = 1
    if os.path.isfile(dp):
        m = re.search(r'numberOfSubdomains\s+(\d+)\s*;', open(dp).read())
        n = int(m.group(1)) if m else 1
    if n > 1:
        for c in ('decomposePar', 'reconstructPar', 'mpirun'):
            if not shutil.which(c):
                sys.exit('missing command: ' + c)
    cores, model = host()
    os.makedirs(out)
    fields = ['case', 'repeat', 'processes', 'status', 'load1_before', 't_blockMesh_s', 't_checkMesh_s', 't_decomposePar_s',
              't_solver_wall_s', 't_reconstructPar_s', 'log_ExecutionTime_s', 'log_ClockTime_s', 'steps', 'host_cores', 'cpu_model']
    with open(os.path.join(out, 'timing.csv'), 'x', newline='', encoding='utf-8') as f:
        w = csv.DictWriter(f, fieldnames=fields)
        w.writeheader()
        for r in range(1, a.repeats + 1):
            rep = os.path.join(out, f'rep{r:02d}')
            shutil.copytree(case, rep)
            row = dict(case=os.path.basename(case), repeat=r, processes=n, load1_before=round(os.getloadavg()[0], 2),
                       host_cores=cores, cpu_model=model)
            steps = [('t_blockMesh_s', ['blockMesh', '-dict', os.path.join(tut, 'resources', 'blockMesh', 'pitzDaily')], 'log.blockMesh'),
                     ('t_checkMesh_s', ['checkMesh'], 'log.checkMesh')]
            if n > 1:
                steps += [('t_decomposePar_s', ['decomposePar'], 'log.decomposePar'),
                          ('t_solver_wall_s', ['mpirun', '-np', str(n)] + shlex.split(a.mpirun_args) + ['foamRun', '-parallel'], 'log.foamRun'),
                          ('t_reconstructPar_s', ['reconstructPar', '-latestTime'], 'log.reconstructPar')]
            else:
                steps += [('t_solver_wall_s', ['foamRun'], 'log.foamRun')]
            status = 'ok'
            for key, cmd, log in steps:
                rc, dt = run(cmd, rep, log)
                row[key] = f'{dt:.4f}'
                if rc != 0:
                    status = f'failed in {cmd[0] if cmd[0] != "mpirun" else "foamRun"} (exit {rc})'
                    break
            lt = open(os.path.join(rep, 'log.foamRun'), errors='replace').read() if os.path.isfile(os.path.join(rep, 'log.foamRun')) else ''
            tm = re.findall(r'ExecutionTime = ([0-9.eE+-]+) s\s+ClockTime = ([0-9.eE+-]+) s', lt)
            st = re.findall(r'(?m)^Time = ([0-9.eE+-]+)', lt)
            if tm:
                row['log_ExecutionTime_s'], row['log_ClockTime_s'] = tm[-1]
            row['steps'] = st[-1] if st else ''
            if status == 'ok' and not re.search(r'(?m)^End\s*$', lt):
                status = 'foamRun log without End'
            row['status'] = status
            w.writerow(row)
            f.flush()
            print(f'repeat {r}: {status}, solver wall time {row.get("t_solver_wall_s", "-")} s, ClockTime {row.get("log_ClockTime_s", "-")} s')
            if status != 'ok':
                sys.exit(2)
    print('Wrote', os.path.join(out, 'timing.csv'))


if __name__ == '__main__':
    main()

benchmark_stats.py merges the timing.csv files, checks that all runs have the same number of iterations and, for each number of processes, reports median, min–max, spread, milliseconds per iteration, the log ClockTime, total time with all phases, speedup, efficiency, the Karp–Flatt metric and the Amdahl serial fraction fitted to the data, with predictions for 8 and 16 processes.

#!/usr/bin/env python3
"""Speedup statistics from the timing.csv files of benchmark_repeats.py.

For every number of processes n it takes the successful repetitions and
reports the median, minimum and maximum solver wall time, the time per
iteration, the median "time to solution" (all phases, including blockMesh,
checkMesh, decomposePar and reconstructPar), and, relative to the serial
median T1:

  speedup        S(n) = T1 / Tn
  efficiency     E(n) = S(n) / n
  Karp-Flatt     e(n) = (1/S - 1/n) / (1 - 1/n)   experimentally determined
                 serial fraction (Karp & Flatt, 1990); e growing with n
                 points to parallel overhead rather than serial code
  Amdahl         serial fraction f fitted by least squares on
                 1/S = f + (1 - f)/n, and the speedups it predicts for 8 and
                 16 processes (an upper bound: overhead is not modelled).

It also reports the spread of the repetitions, (max - min)/median, and
compares the medians with the integer ClockTime of the logs.

Usage: python3 benchmark_stats.py reference_bench/timing.csv parallel_2_bench/timing.csv parallel_4_bench/timing.csv [--plot speedup_repeats.png]
"""
import argparse, csv, os, statistics, sys


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('timing', nargs='+')
    ap.add_argument('--plot')
    a = ap.parse_args()
    runs = {}
    for t in a.timing:
        for r in csv.DictReader(open(t, encoding='utf-8')):
            if r['status'] != 'ok':
                continue
            runs.setdefault(int(r['processes']), []).append(r)
    if 1 not in runs:
        sys.exit('a successful serial run (1 process) is required')
    steps = {float(r['steps']) for rr in runs.values() for r in rr}
    if len(steps) != 1:
        sys.exit(f'runs with different numbers of iterations: {sorted(steps)}')
    nsteps = steps.pop()
    stat = {}
    for n, rr in sorted(runs.items()):
        ts = [float(r['t_solver_wall_s']) for r in rr]
        tot = [sum(float(r[k] or 0) for k in ('t_blockMesh_s', 't_checkMesh_s', 't_decomposePar_s', 't_solver_wall_s', 't_reconstructPar_s')) for r in rr]
        ck = [float(r['log_ClockTime_s']) for r in rr if r['log_ClockTime_s']]
        stat[n] = dict(k=len(ts), med=statistics.median(ts), lo=min(ts), hi=max(ts), tot=statistics.median(tot),
                       clock=statistics.median(ck) if ck else float('nan'), load=max(float(r['load1_before']) for r in rr))
    T1 = stat[1]['med']
    print(f'iterations per run: {nsteps:g}; host: {runs[1][0]["host_cores"]} cores, {runs[1][0]["cpu_model"] or "CPU model n/a"}')
    print(f'{"n":>3s} {"runs":>4s} {"median s":>9s} {"min-max s":>15s} {"spread":>7s} {"ms/iter":>8s} {"ClockTime":>9s} {"total s":>8s} '
          f'{"S(n)":>6s} {"E(n)":>6s} {"Karp-Flatt":>10s} {"load1":>6s}')
    xs, ys = [], []
    for n, s in stat.items():
        S = T1 / s['med']
        E = S / n
        kf = (1 / S - 1 / n) / (1 - 1 / n) if n > 1 else float('nan')
        if n > 1:
            xs.append(1 / n)
            ys.append(1 / S)
        s.update(S=S, E=E)
        print(f'{n:3d} {s["k"]:4d} {s["med"]:9.3f} {s["lo"]:7.3f}-{s["hi"]:<7.3f} {(s["hi"] - s["lo"]) / s["med"]:7.1%} '
              f'{1e3 * s["med"] / nsteps:8.3f} {s["clock"]:9.0f} {s["tot"]:8.3f} {S:6.2f} {E:6.1%} '
              + (f'{kf:10.3f}' if n > 1 else f'{"-":>10s}') + f' {s["load"]:6.2f}')
    if xs:
        # 1/S = f + (1 - f)/n  ->  1/S - 1/n = f (1 - 1/n): least squares through the origin
        num = sum((y - x) * (1 - x) for x, y in zip(xs, ys))
        den = sum((1 - x) ** 2 for x in xs)
        f = num / den
        print(f'Amdahl serial fraction fitted on n > 1: f = {f:.3f}; predicted S(8) = {1 / (f + (1 - f) / 8):.2f}, '
              f'S(16) = {1 / (f + (1 - f) / 16):.2f}, limit 1/f = {1 / f if f > 0 else float("inf"):.1f}')
    print('Speedup and efficiency use the median solver wall time; "total" includes mesh, checks, decomposition and reconstruction.')
    if a.plot:
        if os.path.exists(a.plot):
            sys.exit('refusing to overwrite ' + a.plot)
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        ns = sorted(stat)
        fig, ax = plt.subplots(figsize=(6.8, 4.2), dpi=150)
        ax.plot(ns, [stat[n]['S'] for n in ns], 'o-', label='median speedup')
        ax.vlines(ns, [T1 / stat[n]['hi'] for n in ns], [T1 / stat[n]['lo'] for n in ns], color='C0', alpha=0.5, label='range of repetitions')
        ax.plot(ns, ns, '--', color='gray', label='ideal')
        ax.set(xlabel='MPI processes', ylabel='S(n) = T1 / Tn', title='Repeated runs, same mesh and iterations')
        ax.grid(alpha=0.3)
        ax.legend(fontsize=8)
        fig.tight_layout()
        fig.savefig(a.plot)
        print('Wrote', a.plot)


if __name__ == '__main__':
    main()

Test with the OpenFOAM commands replaced by placeholders that write synthetic logs and wait for a time proportional to an imposed cost per iteration (section 15): the numbers show how the scripts work and the difference between the integer log ClockTime and the measured time, not OpenFOAM performance.

$ python3 prepare_campaign.py --out ~/OpenFOAM/guide19_bench --steps 600
Created six independent cases: /home/user/OpenFOAM/guide19_bench
[exit 0]
$ python3 benchmark_repeats.py ~/OpenFOAM/guide19_run01/reference --repeats 3
the source case has already been run or meshed (log.blockMesh, log.foamRun, RUN_STARTED, polyMesh): use a fresh prepared case
[exit 1]
$ python3 benchmark_repeats.py ~/OpenFOAM/guide19_bench/reference --repeats 3
repeat 1: ok, solver wall time 1.2339 s, ClockTime 12 s
repeat 2: ok, solver wall time 1.2350 s, ClockTime 12 s
repeat 3: ok, solver wall time 1.2409 s, ClockTime 12 s
Wrote /home/user/OpenFOAM/guide19_bench/reference_bench/timing.csv
[exit 0]
$ python3 benchmark_repeats.py ~/OpenFOAM/guide19_bench/parallel_2 --repeats 3
repeat 1: ok, solver wall time 0.7437 s, ClockTime 7 s
repeat 2: ok, solver wall time 0.7439 s, ClockTime 7 s
repeat 3: ok, solver wall time 0.7439 s, ClockTime 7 s
Wrote /home/user/OpenFOAM/guide19_bench/parallel_2_bench/timing.csv
[exit 0]
$ python3 benchmark_repeats.py ~/OpenFOAM/guide19_bench/parallel_4 --repeats 3 --mpirun-args "--bind-to core"
repeat 1: ok, solver wall time 0.5082 s, ClockTime 5 s
repeat 2: ok, solver wall time 0.5093 s, ClockTime 5 s
repeat 3: ok, solver wall time 0.5107 s, ClockTime 5 s
Wrote /home/user/OpenFOAM/guide19_bench/parallel_4_bench/timing.csv
[exit 0]
$ ls ~/OpenFOAM/guide19_bench ~/OpenFOAM/guide19_bench/parallel_4_bench ~/OpenFOAM/guide19_bench/parallel_4_bench/rep01
/home/user/OpenFOAM/guide19_bench:
campaign_manifest.csv
limited
parallel_2
parallel_2_bench
parallel_4
parallel_4_bench
reference
reference_bench
relax_06
upwind

/home/user/OpenFOAM/guide19_bench/parallel_4_bench:
rep01
rep02
rep03
timing.csv

/home/user/OpenFOAM/guide19_bench/parallel_4_bench/rep01:
0
EXPERIMENT.txt
constant
log.blockMesh
log.checkMesh
log.decomposePar
log.foamRun
log.reconstructPar
processor0
processor1
processor2
processor3
system
[exit 0]
$ cut -d, -f1-12 ~/OpenFOAM/guide19_bench/parallel_4_bench/timing.csv
case,repeat,processes,status,load1_before,t_blockMesh_s,t_checkMesh_s,t_decomposePar_s,t_solver_wall_s,t_reconstructPar_s,log_ExecutionTime_s,log_ClockTime_s
parallel_4,1,4,ok,0.12,0.0030,0.0015,0.0083,0.5082,0.0015,5.76251,5
parallel_4,2,4,ok,0.12,0.0032,0.0014,0.0093,0.5093,0.0016,5.76251,5
parallel_4,3,4,ok,0.12,0.0036,0.0014,0.0092,0.5107,0.0017,5.76251,5
[exit 0]
$ python3 benchmark_stats.py ~/OpenFOAM/guide19_bench/*_bench/timing.csv --plot ~/OpenFOAM/guide19_bench/speedup_repeats.png
iterations per run: 600; host: 2 cores, Intel(R) Xeon(R) Processor @ 2.10GHz
  n runs  median s       min-max s  spread  ms/iter ClockTime  total s   S(n)   E(n) Karp-Flatt  load1
  1    3     1.235   1.234-1.241      0.6%    2.058        12    1.241   1.00 100.0%          -   0.15
  2    3     0.744   0.744-0.744      0.0%    1.240         7    0.757   1.66  83.0%      0.205   0.14
  4    3     0.509   0.508-0.511      0.5%    0.849         5    0.525   2.42  60.6%      0.217   0.12
Amdahl serial fraction fitted on n > 1: f = 0.213; predicted S(8) = 3.21, S(16) = 3.82, limit 1/f = 4.7
Speedup and efficiency use the median solver wall time; "total" includes mesh, checks, decomposition and reconstruction.
Wrote /home/user/OpenFOAM/guide19_bench/speedup_repeats.png
[exit 0]

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Prepare isolated OpenFOAM 14 numerical-method and MPI experiments.

The source tutorial is never changed. Existing output directories are never replaced.
"""
import argparse
import csv
import os
from pathlib import Path
import re
import shutil
import sys

STUDIES = [
    ("reference", "bounded Gauss linearUpwind grad(U)", 1, "Reference scheme"),
    ("upwind", "bounded Gauss upwind", 1, "First-order upwind"),
    ("limited", "bounded Gauss limitedLinearV 1", 1, "Limited vector scheme"),
    ("relax_06", "bounded Gauss linearUpwind grad(U)", 1, "U relaxation 0.6"),
    ("parallel_2", "bounded Gauss linearUpwind grad(U)", 2, "Reference MPI 2"),
    ("parallel_4", "bounded Gauss linearUpwind grad(U)", 4, "Reference MPI 4"),
]


def replace_single_line(text, pattern, replacement, label):
    value, count = re.subn(pattern, replacement, text, flags=re.M)
    if count != 1:
        raise ValueError(f"Expected exactly one {label} match, found {count}")
    return value


def strip_residual_control(text):
    """Remove one balanced residualControl dictionary for fixed-step benchmarks."""
    matches = list(re.finditer(r"(?m)^\s*residualControl\s*\{", text))
    if len(matches) != 1:
        raise ValueError("Expected exactly one residualControl block")
    match = matches[0]
    open_idx = text.find("{", match.start())
    depth = 0
    end = None
    for idx in range(open_idx, len(text)):
        if text[idx] == "{":
            depth += 1
        elif text[idx] == "}":
            depth -= 1
            if depth == 0:
                end = idx + 1
                break
    if end is None:
        raise ValueError("Unbalanced residualControl block")
    return text[:match.start()] + "\n    // residualControl disabled for equal-step benchmarking.\n" + text[end:]


def prepare_case(src: Path, dst: Path, scheme: str, np: int, relaxation: bool, steps: int):
    if dst.exists():
        raise FileExistsError(f"Refusing to replace existing directory: {dst}")
    dst.mkdir(parents=True)
    try:
        for sub in ("0", "constant", "system"):
            if not (src / sub).is_dir():
                raise FileNotFoundError(f"Tutorial folder missing: {src/sub}")
            shutil.copytree(src / sub, dst / sub)

        schemes = dst / "system/fvSchemes"
        original = schemes.read_text(encoding="utf-8")
        modified = replace_single_line(
            original, r"^\s*div\(phi,U\)\s+[^;]+;",
            f"    div(phi,U)      {scheme};", "div(phi,U)"
        )
        schemes.write_text(modified, encoding="utf-8")

        control = dst / "system/controlDict"
        v = control.read_text(encoding="utf-8")
        v = replace_single_line(v, r"^\s*endTime\s+[^;]+;", f"endTime         {steps};", "endTime")
        v = replace_single_line(v, r"^\s*writeInterval\s+[^;]+;", "writeInterval   100;", "writeInterval")
        control.write_text(v, encoding="utf-8")

        solution = dst / "system/fvSolution"
        v = strip_residual_control(solution.read_text(encoding="utf-8"))
        if relaxation:
            v = replace_single_line(v, r"^\s*U\s+0\.9;[^\n]*$", "        U               0.6; // Method experiment", "relaxation U")
        solution.write_text(v, encoding="utf-8")

        if np > 1:
            (dst / "system/decomposeParDict").write_text(
                "FoamFile\n{\n    format ascii;\n    class dictionary;\n    object decomposeParDict;\n}\n"
                f"numberOfSubdomains {np};\nmethod scotch;\n", encoding="utf-8")
        (dst / "EXPERIMENT.txt").write_text(
            "OpenFOAM 14. Isolated teaching case, derived from "
            "tutorials/incompressibleFluid/pitzDailySteady.\n"
            f"Scheme: {scheme}\nMPI processes: {np}\nFixed pseudo-iterations: {steps}\n"
            "In fvSolution, SIMPLE/residualControl is removed to permit fixed-step timing.\n"
            "The numerical time index is NOT a physical time in this steady case.\n"
            "Real CFD solver execution required for scientific claims.\n", encoding="utf-8")
    except Exception:
        # Only remove a directory created by this invocation; source is untouched.
        shutil.rmtree(dst)
        raise


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--out", required=True, help="New directory that does not exist yet")
    ap.add_argument("--source", default=None, help="Override source tutorial directory")
    ap.add_argument("--steps", type=int, default=600, help="Fixed pseudo-iterations (default 600)")
    args = ap.parse_args()
    if args.steps < 10 or args.steps > 100000:
        ap.error("--steps must be between 10 and 100000")
    default_src = Path(os.environ.get("FOAM_TUTORIALS", "")) / "incompressibleFluid" / "pitzDailySteady"
    src = Path(args.source) if args.source else default_src
    if not args.source and not os.environ.get("FOAM_TUTORIALS"):
        ap.error("Set FOAM_TUTORIALS in your OpenFOAM shell, or supply --source")
    src = src.expanduser().resolve()
    out = Path(args.out).expanduser().resolve()
    if not src.is_dir():
        ap.error(f"Source tutorial not found: {src}")
    if out.exists():
        ap.error(f"Refusing to overwrite existing output: {out}")
    out.mkdir(parents=True)
    try:
        with (out / "campaign_manifest.csv").open("w", newline="", encoding="utf-8") as f:
            writer = csv.writer(f)
            writer.writerow(["case", "scheme", "mpi_processes", "steps", "notes"])
            for name, scheme, np, note in STUDIES:
                prepare_case(src, out / name, scheme, np, name == "relax_06", args.steps)
                writer.writerow([name, scheme, np, args.steps, note])
        print(f"Created six independent cases: {out}")
    except Exception as exc:
        # Existing directory was checked before creation; remove only our new tree.
        shutil.rmtree(out)
        print(f"Preparation failed safely: {exc}", file=sys.stderr)
        raise SystemExit(2)


if __name__ == "__main__":
    main()

• It copies only 0/, constant/ and system/ of the tutorial; every substitution must find exactly one line (div(phi,U), endTime, writeInterval, the U factor at 0.9) and a single residualControl block, otherwise it stops and deletes only the folder it created.
• The diff with the tutorial shows only the expected changes, plus a cosmetic side effect: the ^\s* expression with re.M also captures the preceding empty lines, which are removed in front of endTime and writeInterval. The dictionary content does not change.
• writeInterval stays 100 (as in the tutorial): with --steps not a multiple of 100 the last iteration may not be written to disk; use multiples of 100.

Test with the version 14 tutorial (FOAM_TUTORIALS set as etc/bashrc does):

$ python3 prepare_campaign.py --out ~/OpenFOAM/guide19_run01 --steps 600
Created six independent cases: /home/user/OpenFOAM/guide19_run01
[exit 0]
$ ls ~/OpenFOAM/guide19_run01 ~/OpenFOAM/guide19_run01/parallel_4/system
/home/user/OpenFOAM/guide19_run01:
campaign_manifest.csv
limited
parallel_2
parallel_4
reference
relax_06
upwind

/home/user/OpenFOAM/guide19_run01/parallel_4/system:
blockMeshDict
controlDict
decomposeParDict
functions
fvSchemes
fvSolution
[exit 0]
$ cat ~/OpenFOAM/guide19_run01/campaign_manifest.csv
case,scheme,mpi_processes,steps,notes
reference,bounded Gauss linearUpwind grad(U),1,600,Reference scheme
upwind,bounded Gauss upwind,1,600,First-order upwind
limited,bounded Gauss limitedLinearV 1,1,600,Limited vector scheme
relax_06,bounded Gauss linearUpwind grad(U),1,600,U relaxation 0.6
parallel_2,bounded Gauss linearUpwind grad(U),2,600,Reference MPI 2
parallel_4,bounded Gauss linearUpwind grad(U),4,600,Reference MPI 4
[exit 0]
$ cat ~/OpenFOAM/guide19_run01/parallel_4/EXPERIMENT.txt
OpenFOAM 14. Isolated teaching case, derived from tutorials/incompressibleFluid/pitzDailySteady.
Scheme: bounded Gauss linearUpwind grad(U)
MPI processes: 4
Fixed pseudo-iterations: 600
In fvSolution, SIMPLE/residualControl is removed to permit fixed-step timing.
The numerical time index is NOT a physical time in this steady case.
Real CFD solver execution required for scientific claims.
[exit 0]
$ diff -r "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady/system" ~/OpenFOAM/guide19_run01/reference/system
diff -r /opt/OpenFOAM-14/tutorials/incompressibleFluid/pitzDailySteady/system/controlDict /home/user/OpenFOAM/guide19_run01/reference/system/controlDict
24,25c24
< 
< endTime         2000;
---
> endTime         600;
30d28
< 
diff -r /opt/OpenFOAM-14/tutorials/incompressibleFluid/pitzDailySteady/system/fvSolution /home/user/OpenFOAM/guide19_run01/reference/system/fvSolution
49,54c49,50
<     residualControl
<     {
<         p               1e-2;
<         U               1e-3;
<         "(k|epsilon|omega|f|v2)" 1e-3;
<     }
---
>     // residualControl disabled for equal-step benchmarking.
> 
[exit 1]
$ diff ~/OpenFOAM/guide19_run01/reference/system/fvSchemes ~/OpenFOAM/guide19_run01/upwind/system/fvSchemes
30c30
<     div(phi,U)      bounded Gauss linearUpwind grad(U);
---
>     div(phi,U)      bounded Gauss upwind;
[exit 1]
$ diff ~/OpenFOAM/guide19_run01/reference/system/fvSchemes ~/OpenFOAM/guide19_run01/limited/system/fvSchemes
30c30
<     div(phi,U)      bounded Gauss linearUpwind grad(U);
---
>     div(phi,U)      bounded Gauss limitedLinearV 1;
[exit 1]
$ diff ~/OpenFOAM/guide19_run01/reference/system/fvSolution ~/OpenFOAM/guide19_run01/relax_06/system/fvSolution
57c57
<         U               0.9; // 0.9 is more stable but 0.95 more convergent
---
>         U               0.6; // Method experiment
[exit 1]
$ cat ~/OpenFOAM/guide19_run01/parallel_4/system/decomposeParDict
FoamFile
{
    format ascii;
    class dictionary;
    object decomposeParDict;
}
numberOfSubdomains 4;
method scotch;
[exit 0]
$ python3 prepare_campaign.py --out ~/OpenFOAM/guide19_run01 --steps 600
usage: prepare_campaign.py [-h] --out OUT [--source SOURCE] [--steps STEPS]
prepare_campaign.py: error: Refusing to overwrite existing output: /home/user/OpenFOAM/guide19_run01
[exit 2]

Physical quantities must be defined before the test (section 13). The script adds to the tutorial system/functions, kept as it is (copy in functions.tutorial), standard functions from etc/caseDicts/functions with the #includeFunc argument syntax: inlet and outlet flow rates, mean pressure on both sections, three probes of U and p, the wall shear stress and its profile along x on the lower wall (graphPatchCutLayerAverage). It refuses cases already started, meshed or with a functions different from the tutorial; it checks all cases before changing any.

#!/usr/bin/env python3
"""Add the Guide 19 physical observables to NEW, never-run pitzDailySteady copies.

The tutorial system/functions (streamlines) is kept verbatim and saved as
system/functions.tutorial; the following OpenFOAM 14 functions are
appended, all from etc/caseDicts/functions:

  patchFlowRate   inlet and outlet volumetric flow rate (sum of phi) [m3/s]
  patchAverage    area-averaged kinematic pressure p on inlet and outlet
                  -> pressure drop [m2/s2]; multiply by rho for Pa
  probes          U and p at three points: recirculation, reattachment
                  region, downstream channel
  wallShearStress + graphPatchCutLayerAverage on lowerWall along x
                  -> wall shear stress profile, from which observables.py
                  finds the reattachment point

Flow rates, pressures and probes are written at every iteration, the wall
shear stress graph at every write time. The physics of the case is not
changed. Cases with logs, RUN_STARTED, a mesh, processor directories or a
modified system/functions are refused.

Usage: python3 add_observables.py CASE [CASE ...]   (or CAMPAIGN_DIR/*)
       python3 add_observables.py --print            (show the appended block)
"""
import glob, os, sys

MARK = '// Guide 19 observables (added by add_observables.py)'
BLOCK = MARK + '''
// Volumetric flow rates [m3/s]: inlet negative (entering), outlet positive.
#includeFunc patchFlowRate(name=inletFlowRate, patch=inlet)
#includeFunc patchFlowRate(name=outletFlowRate, patch=outlet)

// Area-averaged kinematic pressure [m2/s2] on inlet and outlet.
#includeFunc patchAverage(name=inletPressure, patch=inlet, p)
#includeFunc patchAverage(name=outletPressure, patch=outlet, p)

// Probes: recirculation (x = 50 mm, y = -12.5 mm), reattachment region
// (x = 150 mm, y = -20 mm) and downstream channel (x = 250 mm, y = 0).
#includeFunc probes(name=probes, points=((0.05 -0.0125 0) (0.15 -0.02 0) (0.25 0 0)), U, p)

// Wall shear stress on the lower wall, averaged in 200 layers along x.
#includeFunc wallShearStress
#includeFunc graphPatchCutLayerAverage(name=lowerWallShearStress, patch=lowerWall, direction=(1 0 0), nPoints=200, wallShearStress)
'''


def check(case):
    if not all(os.path.isdir(os.path.join(case, d)) for d in ('0', 'constant', 'system')):
        return 'not an OpenFOAM case'
    bad = [f for f in ('log.blockMesh', 'log.foamRun', 'RUN_STARTED', 'RUN_SUCCESS', 'postProcessing') if os.path.exists(os.path.join(case, f))]
    bad += [os.path.basename(p) for p in glob.glob(os.path.join(case, 'processor*'))]
    if os.path.isdir(os.path.join(case, 'constant', 'polyMesh')):
        bad.append('constant/polyMesh')
    if bad:
        return 'already started or meshed (' + ', '.join(bad) + ')'
    fn = os.path.join(case, 'system', 'functions')
    if not os.path.isfile(fn):
        return 'system/functions missing'
    t = open(fn, encoding='utf-8').read()
    if MARK in t:
        return 'observables already installed'
    if 'streamlinesLine' not in t:
        return 'system/functions differs from the pitzDailySteady tutorial'
    if os.path.exists(fn + '.tutorial'):
        return 'system/functions.tutorial already exists'
    return None


def main():
    args = sys.argv[1:]
    if args == ['--print']:
        sys.stdout.write(BLOCK)
        return
    if not args:
        sys.exit(__doc__)
    cases = [c for c in args if os.path.isdir(c)]
    errors = {c: check(c) for c in cases}
    errors.update({c: 'not a directory' for c in args if c not in cases})
    for c, e in errors.items():
        if e:
            print(f'{c}: {e}')
    if any(errors.values()):
        sys.exit('nothing changed')
    for c in cases:
        fn = os.path.join(c, 'system', 'functions')
        t = open(fn, encoding='utf-8').read()
        os.rename(fn, fn + '.tutorial')
        end = '// ************************************************************************* //'
        body = t.rsplit(end, 1)[0].rstrip('\n') if end in t else t.rstrip('\n')
        with open(fn, 'w', encoding='utf-8') as f:
            f.write(body + '\n\n' + BLOCK + '\n' + end + '\n')
        print(f'{c}: observables installed (original kept as system/functions.tutorial)')


if __name__ == '__main__':
    main()

The resulting file (extra/functions_guide19 in the archive):

#includeFunc streamlinesLine
(
    name=streamlines,
    start=(-0.0205 0.001 0.00001),
    end=(-0.0205 0.0251 0.00001),
    nPoints=10,
    fields=(p k U)
)

// Guide 19 observables (added by add_observables.py)
// Volumetric flow rates [m3/s]: inlet negative (entering), outlet positive.
#includeFunc patchFlowRate(name=inletFlowRate, patch=inlet)
#includeFunc patchFlowRate(name=outletFlowRate, patch=outlet)

// Area-averaged kinematic pressure [m2/s2] on inlet and outlet.
#includeFunc patchAverage(name=inletPressure, patch=inlet, p)
#includeFunc patchAverage(name=outletPressure, patch=outlet, p)

// Probes: recirculation (x = 50 mm, y = -12.5 mm), reattachment region
// (x = 150 mm, y = -20 mm) and downstream channel (x = 250 mm, y = 0).
#includeFunc probes(name=probes, points=((0.05 -0.0125 0) (0.15 -0.02 0) (0.25 0 0)), U, p)

// Wall shear stress on the lower wall, averaged in 200 layers along x.
#includeFunc wallShearStress
#includeFunc graphPatchCutLayerAverage(name=lowerWallShearStress, patch=lowerWall, direction=(1 0 0), nPoints=200, wallShearStress)

Test on a new campaign (the first command refuses the list because it contains an already computed case and changes nothing):

$ python3 prepare_campaign.py --out ~/OpenFOAM/guide19_run02 --steps 600
Created six independent cases: /home/user/OpenFOAM/guide19_run02
[exit 0]
$ python3 add_observables.py ~/OpenFOAM/guide19_run02/reference ~/OpenFOAM/guide19_run01/upwind
/home/user/OpenFOAM/guide19_run01/upwind: already started or meshed (log.blockMesh, log.foamRun, RUN_STARTED, RUN_SUCCESS, constant/polyMesh)
nothing changed
[exit 1]
$ python3 add_observables.py ~/OpenFOAM/guide19_run02/reference ~/OpenFOAM/guide19_run02/upwind ~/OpenFOAM/guide19_run02/limited ~/OpenFOAM/guide19_run02/relax_06
/home/user/OpenFOAM/guide19_run02/reference: observables installed (original kept as system/functions.tutorial)
/home/user/OpenFOAM/guide19_run02/upwind: observables installed (original kept as system/functions.tutorial)
/home/user/OpenFOAM/guide19_run02/limited: observables installed (original kept as system/functions.tutorial)
/home/user/OpenFOAM/guide19_run02/relax_06: observables installed (original kept as system/functions.tutorial)
[exit 0]
$ python3 add_observables.py ~/OpenFOAM/guide19_run02/reference
/home/user/OpenFOAM/guide19_run02/reference: observables already installed
nothing changed
[exit 1]

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run an unexecuted prepared case. Never clean or replace earlier data.
set -Eeuo pipefail
if (( $# < 1 || $# > 2 )); then
  echo "Usage: $0 /path/to/CASE [number-of-MPI-processes]" >&2
  exit 64
fi
CASE_DIR="$(realpath "$1")"
NPROC="${2:-1}"
[[ -d "$CASE_DIR/0" && -d "$CASE_DIR/system" && -d "$CASE_DIR/constant" ]] || { echo "Not an OpenFOAM case" >&2; exit 64; }
[[ "$NPROC" =~ ^(1|2|4)$ ]] || { echo "MPI count must be 1, 2, or 4" >&2; exit 64; }
[[ -n "${FOAM_TUTORIALS:-}" ]] || { echo "Source OpenFOAM environment first" >&2; exit 64; }
for cmd in blockMesh checkMesh foamRun; do command -v "$cmd" >/dev/null || { echo "Missing: $cmd" >&2; exit 69; }; done
if (( NPROC > 1 )); then
  for cmd in decomposePar reconstructPar mpirun; do command -v "$cmd" >/dev/null || { echo "Missing: $cmd" >&2; exit 69; }; done
  [[ -f "$CASE_DIR/system/decomposeParDict" ]] || { echo "Missing decomposeParDict" >&2; exit 64; }
  grep -Eq "numberOfSubdomains[[:space:]]+$NPROC;" "$CASE_DIR/system/decomposeParDict" || { echo "MPI count mismatch" >&2; exit 64; }
fi
[[ ! -e "$CASE_DIR/log.blockMesh" && ! -e "$CASE_DIR/log.foamRun" && ! -e "$CASE_DIR/RUN_STARTED" ]] || { echo "Case has already been started; refusing to overwrite" >&2; exit 73; }
shopt -s nullglob
procs=("$CASE_DIR"/processor[0-9]*)
(( ${#procs[@]} == 0 )) || { echo "Processor directories already exist" >&2; exit 73; }
[[ ! -d "$CASE_DIR/constant/polyMesh" ]] || { echo "Mesh already exists; refusing to overwrite" >&2; exit 73; }
cd "$CASE_DIR"
printf 'Started UTC: %s\nMPI processes: %s\n' "$(date -u +%FT%TZ)" "$NPROC" > RUN_STARTED
trap 'rc=$?; printf "Exit: %s\nFinished UTC: %s\n" "$rc" "$(date -u +%FT%TZ)" >> RUN_STARTED; exit "$rc"' EXIT
blockMesh -dict "$FOAM_TUTORIALS/resources/blockMesh/pitzDaily" > log.blockMesh 2>&1
checkMesh > log.checkMesh 2>&1
if (( NPROC == 1 )); then
  foamRun > log.foamRun 2>&1
else
  decomposePar > log.decomposePar 2>&1
  mpirun -np "$NPROC" foamRun -parallel > log.foamRun 2>&1
  # Reconstruction takes place after the measured solver run.
  reconstructPar -latestTime > log.reconstructPar 2>&1
fi
printf 'SUCCESS\n' > RUN_SUCCESS

• Exit codes: 64 for usage errors (folder, number of processes other than 1, 2, 4 or different from the decomposeParDict), 69 for missing commands, 73 if the case has already been started, has a mesh or processor* folders.
• RUN_STARTED receives the exit code and end time on exit even on error; RUN_SUCCESS is written only if all commands succeed (set -Eeuo pipefail).
• As in the previous guide, checkMesh always exits with code 0 even when checks fail (checkMesh.C): check by hand grep -E "Mesh OK|Failed" log.checkMesh and that log.foamRun ends with “End”.

Test with the OpenFOAM commands replaced by placeholders (same sequence of files and checks; the second launch of reference and parallel_2 with 4 processes show the refusals):

$ bash run_case.sh ~/OpenFOAM/guide19_run01/reference 1
[exit 0]
$ ls ~/OpenFOAM/guide19_run01/reference
0
EXPERIMENT.txt
RUN_STARTED
RUN_SUCCESS
constant
log.blockMesh
log.checkMesh
log.foamRun
system
[exit 0]
$ cat ~/OpenFOAM/guide19_run01/reference/RUN_STARTED ~/OpenFOAM/guide19_run01/reference/RUN_SUCCESS
Started UTC: 2026-10-31T16:39:11Z
MPI processes: 1
Exit: 0
Finished UTC: 2026-10-31T16:39:11Z
SUCCESS
[exit 0]
$ bash run_case.sh ~/OpenFOAM/guide19_run01/reference 1
Case has already been started; refusing to overwrite
[exit 73]
$ bash run_case.sh ~/OpenFOAM/guide19_run01/upwind 1
[exit 0]
$ bash run_case.sh ~/OpenFOAM/guide19_run01/limited 1
[exit 0]
$ bash run_case.sh ~/OpenFOAM/guide19_run01/relax_06 1
[exit 0]
$ bash run_case.sh ~/OpenFOAM/guide19_run01/parallel_2 4
MPI count mismatch
[exit 64]
$ bash run_case.sh ~/OpenFOAM/guide19_run01/parallel_2 2
[exit 0]
$ bash run_case.sh ~/OpenFOAM/guide19_run01/parallel_4 4
[exit 0]
$ ls ~/OpenFOAM/guide19_run01/parallel_4
0
EXPERIMENT.txt
RUN_STARTED
RUN_SUCCESS
constant
log.blockMesh
log.checkMesh
log.decomposePar
log.foamRun
log.reconstructPar
processor0
processor1
processor2
processor3
system
[exit 0]

• Before the test write down which quantities decide the comparison and with what tolerance, for example: pressure drop to 4 significant digits, x_r/h to 2 decimals, flow balance within 10⁻⁴. Install the monitors with add_observables.py on both copies before the run.
• After the run: python3 residual_history.py ~/OpenFOAM/guide19_run01 --cases reference relax_06 --plot (residuals at the same iteration and the “stop@” column) and python3 observables.py ~/OpenFOAM/guide19_run01 --cases reference relax_06 --plot (variation of the observables over the last 100 iterations and x_r at every write time).
• If the observables in relax_06 are still changing, prepare a new campaign with more iterations, for example python3 prepare_campaign.py --out ~/OpenFOAM/guide19_run02 --steps 2000, instead of restarting the cases already used. At convergence both variants must give the same values (section 6): a persistent difference indicates a convergence problem, not a physical effect of the relaxation.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Parse solver logs and optionally plot MPI speedup; no CFD execution performed."""
import argparse
import csv
from pathlib import Path
import re

TIME_RE = re.compile(r"(?:^|\n)\s*Time\s*=\s*([0-9.eE+-]+)")
CLOCK_RE = re.compile(r"ExecutionTime\s*=\s*([0-9.eE+-]+)\s*s\s+ClockTime\s*=\s*([0-9.eE+-]+)\s*s")
RESID_RE = re.compile(r"Solving for\s+([^,]+),\s+Initial residual\s*=\s*([0-9.eE+-]+)")
CONT_RE = re.compile(r"time step continuity errors\s*:\s*sum local\s*=\s*([0-9.eE+-]+)")


def parse_case(base, row):
    name = row['case']
    p = base / name
    result = dict(row)
    result.update(success=str((p/'RUN_SUCCESS').exists()), last_step='', clock_s='', cpu_s='',
                  last_Ux_initial='', last_p_initial='', last_continuity_local='', log_present='false')
    log = p/'log.foamRun'
    if not log.is_file():
        return result
    s = log.read_text(errors='replace')
    result['log_present']='true'
    times = TIME_RE.findall(s)
    clocks = CLOCK_RE.findall(s)
    resids = RESID_RE.findall(s)
    conts = CONT_RE.findall(s)
    if times:
        result['last_step'] = times[-1]
    if clocks:
        result['cpu_s'], result['clock_s'] = clocks[-1]
    for name_field, key in [('Ux','last_Ux_initial'),('p','last_p_initial')]:
        values = [res for field,res in resids if field.strip() == name_field]
        if values:
            result[key] = values[-1]
    if conts:
        result['last_continuity_local'] = conts[-1]
    return result


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument('campaign', type=Path)
    ap.add_argument('--plot', action='store_true', help='Requires matplotlib')
    args = ap.parse_args()
    base = args.campaign.expanduser().resolve()
    with (base/'campaign_manifest.csv').open(newline='', encoding='utf-8') as f:
        manifest = list(csv.DictReader(f))
    rows = [parse_case(base, row) for row in manifest]
    out=base/'results_summary.csv'
    if out.exists():
        ap.error('Refusing to overwrite results_summary.csv; rename or archive it first')
    with out.open('w',newline='',encoding='utf-8') as f:
        wr=csv.DictWriter(f,fieldnames=list(rows[0]));wr.writeheader();wr.writerows(rows)
    print('Wrote:',out)
    if args.plot:
        try:
            import matplotlib.pyplot as plt
        except ImportError as exc:
            ap.error(f'matplotlib is required for --plot: {exc}')
        pts = []
        for r in rows:
            if r['case'] in ('reference','parallel_2','parallel_4') and r['success']=='True' and r['clock_s']:
                pts.append((int(r['mpi_processes']),float(r['clock_s']),r['last_step'],int(r['steps'])))
        if len(pts) == 3 and all(float(x[2]) >= x[3] for x in pts):
            pts.sort()
            t1=next(t for n,t,_,_ in pts if n==1)
            fig,ax=plt.subplots(figsize=(6.8,4.2))
            ax.plot([p[0] for p in pts],[t1/p[1] for p in pts],marker='o',label='Measured speedup')
            ax.plot([1,2,4],[1,2,4],ls='--',color='gray',label='Ideal speedup')
            ax.set(xlabel='MPI processes',ylabel='S(n) = T(1)/T(n)',title='OpenFOAM 14: fixed-mesh MPI study')
            ax.grid(alpha=.25);ax.legend();fig.tight_layout()
            png=base/'speedup.png'
            if png.exists():
                ap.error('Refusing to overwrite speedup.png')
            fig.savefig(png,dpi=180)
            print('Wrote:',png)
        else:
            print('Speedup plot skipped: need 3 successful runs at the requested step count.')

if __name__ == '__main__':
    main()

• It reads the last Time, the last ExecutionTime/ClockTime pair, the last initial residuals of Ux and p and the last local continuity error; it writes results_summary.csv and plots the speedup only if the three cases of the MPI comparison succeeded and reached the requested iterations.
• The speedup uses ClockTime, which OpenFOAM 14 prints in whole seconds (section 10), and a single run per case: read it as an indication; for the median and millisecond resolution use benchmark_repeats.py and benchmark_stats.py.
• The summary is written before the plot and never overwritten: to also get speedup.png use --plot on the first run, otherwise rename results_summary.csv and repeat.

Test on the synthetic logs produced by the placeholders (the integer ClockTime of the synthetic log reproduces the resolution of the real one):

$ python3 analyze_campaign.py ~/OpenFOAM/guide19_run01 --plot
Wrote: /home/user/OpenFOAM/guide19_run01/results_summary.csv
Wrote: /home/user/OpenFOAM/guide19_run01/speedup.png
[exit 0]
$ cat ~/OpenFOAM/guide19_run01/results_summary.csv
case,scheme,mpi_processes,steps,notes,success,last_step,clock_s,cpu_s,last_Ux_initial,last_p_initial,last_continuity_local,log_present
reference,bounded Gauss linearUpwind grad(U),1,600,Reference scheme,True,600,12,12.2238,3.80923e-05,0.000327463,0.000547599,true
upwind,bounded Gauss upwind,1,600,First-order upwind,True,600,12,12.2238,2.49677e-06,1.69327e-05,2.02735e-05,true
limited,bounded Gauss limitedLinearV 1,1,600,Limited vector scheme,True,600,12,12.2238,9.74942e-06,0.000103747,0.000167696,true
relax_06,bounded Gauss linearUpwind grad(U),1,600,U relaxation 0.6,True,600,12,12.2238,0.00384015,0.0116862,0.0198363,true
parallel_2,bounded Gauss linearUpwind grad(U),2,600,Reference MPI 2,True,600,7,7.69173,3.80923e-05,0.000327463,0.000547599,true
parallel_4,bounded Gauss linearUpwind grad(U),4,600,Reference MPI 4,True,600,5,5.76251,3.80923e-05,0.000327463,0.000547599,true
[exit 0]
$ python3 analyze_campaign.py ~/OpenFOAM/guide19_run01 --plot
usage: analyze_campaign.py [-h] [--plot] campaign
analyze_campaign.py: error: Refusing to overwrite results_summary.csv; rename or archive it first
[exit 2]

The script reads the files of the add_observables.py monitors and reports, at the last iteration and as relative variation over the last --window iterations: inlet and outlet flow rates and their imbalance, kinematic pressure drop (in Pa with --rho), Ux at the three probes, and the reattachment point x_r on the lower wall at every write time. x_r is the last sign change of the x component of the wall shear stress before the region with the downstream sign: the computation does not depend on the sign convention of wallShearStress and excludes the small corner eddy near the step. The columns of the .xy file are recognised from the header written by setWriter.C, which shortens names that are too long with “…”.

#!/usr/bin/env python3
"""Physical observables of the Guide 19 campaign from the add_observables.py monitors.

For every case (campaign_manifest.csv or --cases) it reads postProcessing/
and reports at the last iteration, and their variation over the last
--window iterations (default 100):

  * inlet and outlet volumetric flow rates and their imbalance;
  * pressure drop between inlet and outlet (kinematic, and in Pa with --rho);
  * Ux at the three probes;
  * reattachment point x_r of the separated flow on the lower wall, found
    where the x-wall-shear-stress changes sign for the last time
    (independent of the sign convention), and x_r/h with h = --h;
    x_r is reported for every write time, to show whether it is still moving.

Outputs observables_summary.csv and, with --plot, observables.png.
Existing outputs are never overwritten.

Usage: python3 observables.py CAMPAIGN_DIR [--cases ...] [--rho 1.2] [--h 0.0254] [--window 100] [--plot]
"""
import argparse, csv, glob, math, os, re, sys


def read_dat(case, name):
    out = []
    for d in sorted(glob.glob(os.path.join(case, 'postProcessing', name, '*')), key=lambda p: float(os.path.basename(p))):
        for f in glob.glob(os.path.join(d, '*.dat')):
            for line in open(f, encoding='utf-8'):
                if line.startswith('#') or not line.strip():
                    continue
                v = line.split()
                out.append((float(v[0]), float(v[1])))
    return sorted(out)


def read_probe(case, field='U'):
    rows = []
    for d in sorted(glob.glob(os.path.join(case, 'postProcessing', 'probes', '*')), key=lambda p: float(os.path.basename(p))):
        f = os.path.join(d, field)
        if not os.path.isfile(f):
            continue
        for line in open(f, encoding='utf-8'):
            if line.startswith('#') or not line.strip():
                continue
            t = float(line.split()[0])
            vecs = re.findall(r'\(([^)]*)\)', line)
            rows.append((t, [[float(c) for c in v.split()] for v in vecs]))
    return rows


def reattachment(xyfile, h):
    head, data = None, []
    for line in open(xyfile, encoding='utf-8'):
        if line.startswith('#'):
            head = line[1:].split()
            continue
        if line.strip():
            data.append([float(v) for v in line.split()])
    ix, it = 1, 4
    if head:
        if 'x' in head:
            ix = head.index('x')
        cand = [i for i, h_ in enumerate(head) if h_.startswith('wallShe')]
        if cand:
            it = cand[0]
    pts = sorted((r[ix], r[it]) for r in data if r[ix] > 1e-3)
    if len(pts) < 5:
        return None
    xmax = pts[-1][0]
    tail = [t for x, t in pts if x > 0.9 * xmax]
    sd = 1 if sum(tail) > 0 else -1
    xr = None
    for (x0, t0), (x1, t1) in zip(pts[:-1], pts[1:]):
        if t0 * sd < 0 <= t1 * sd:
            xr = x0 + (x1 - x0) * (-t0) / (t1 - t0) if t1 != t0 else x0
    return xr


def variation(series, window):
    if not series:
        return None, None
    tl = series[-1][0]
    w = [v for t, v in series if t >= tl - window]
    mean = sum(w) / len(w)
    return series[-1][1], (max(w) - min(w)) / abs(mean) if mean else float('inf')


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('campaign')
    ap.add_argument('--cases', nargs='*')
    ap.add_argument('--rho', type=float, help='density [kg/m3] to convert the kinematic pressure drop to Pa')
    ap.add_argument('--h', type=float, default=0.0254, help='step height [m]')
    ap.add_argument('--window', type=float, default=100)
    ap.add_argument('--plot', action='store_true')
    a = ap.parse_args()
    base = os.path.abspath(os.path.expanduser(a.campaign))
    if a.cases:
        cases = a.cases
    else:
        with open(os.path.join(base, 'campaign_manifest.csv'), encoding='utf-8') as f:
            cases = [r['case'] for r in csv.DictReader(f)]
    out_csv = os.path.join(base, 'observables_summary.csv')
    out_png = os.path.join(base, 'observables.png')
    for p in (out_csv,) + ((out_png,) if a.plot else ()):
        if os.path.exists(p):
            sys.exit('refusing to overwrite ' + p)
    rows, hist = [], {}
    for c in cases:
        cd = os.path.join(base, c)
        qi, qo = read_dat(cd, 'inletFlowRate'), read_dat(cd, 'outletFlowRate')
        pi, po = read_dat(cd, 'inletPressure'), read_dat(cd, 'outletPressure')
        if not (qi and qo and pi and po):
            print(f'{c}: observables not found (were they installed with add_observables.py before the run?)')
            continue
        po_d = dict(po)
        dp = [(t, v - po_d[t]) for t, v in pi if t in po_d]
        dp_last, dp_var = variation(dp, a.window)
        imb = (qo[-1][1] + qi[-1][1]) / abs(qi[-1][1])
        pr = read_probe(cd, 'U')
        probes = []
        for k in range(len(pr[-1][1]) if pr else 0):
            last, var = variation([(t, v[k][0]) for t, v in pr], a.window)
            probes.append((last, var))
        xr = []
        for d in sorted(glob.glob(os.path.join(cd, 'postProcessing', 'lowerWallShearStress', '*')), key=lambda p: float(os.path.basename(p))):
            files = glob.glob(os.path.join(d, '*.xy'))
            if files:
                xr.append((float(os.path.basename(d)), reattachment(files[0], a.h)))
        hist[c] = dict(dp=dp, xr=xr)
        row = dict(case=c, last_iteration=int(dp[-1][0]), Q_inlet_m3_s=qi[-1][1], Q_outlet_m3_s=qo[-1][1], imbalance=imb,
                   dp_kinematic_m2_s2=dp_last, dp_Pa=(dp_last * a.rho if a.rho else ''), dp_rel_variation=dp_var)
        for k, (last, var) in enumerate(probes):
            row[f'Ux_probe{k}'] = last
            row[f'Ux_probe{k}_rel_variation'] = var
        row['x_r_m'] = xr[-1][1] if xr else ''
        row['x_r_over_h'] = xr[-1][1] / a.h if xr and xr[-1][1] else ''
        row['x_r_history'] = ' '.join(f'{int(t)}:{(v / a.h):.3f}' for t, v in xr if v)
        rows.append(row)
    if not rows:
        sys.exit('no observables found')
    print(f'variations are (max - min)/|mean| over the last {a.window:g} iterations')
    print(f'{"case":12s} {"it":>4s} {"Q_in m3/s":>11s} {"imbal.":>9s} {"dp m2/s2":>9s} {"var":>8s}'
          + (f' {"dp Pa":>8s}' if a.rho else '') + f' {"Ux0":>7s} {"Ux1":>7s} {"Ux2":>7s} {"x_r/h":>6s}  x_r/h at write times')
    for r in rows:
        print(f'{r["case"]:12s} {r["last_iteration"]:4d} {r["Q_inlet_m3_s"]:11.4e} {r["imbalance"]:9.2e} {r["dp_kinematic_m2_s2"]:9.4f} {r["dp_rel_variation"]:8.1e}'
              + (f' {r["dp_Pa"]:8.3f}' if a.rho else '')
              + ''.join(f' {r.get(f"Ux_probe{k}", float("nan")):7.3f}' for k in range(3))
              + f' {r["x_r_over_h"] if r["x_r_over_h"] == "" else format(r["x_r_over_h"], "6.3f"):>6}  {r["x_r_history"]}')
    keys = list(dict.fromkeys(k for r in rows for k in r))
    with open(out_csv, 'x', newline='', encoding='utf-8') as f:
        w = csv.DictWriter(f, fieldnames=keys)
        w.writeheader()
        w.writerows(rows)
    print('Wrote', out_csv)
    if a.plot:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        fig, axs = plt.subplots(1, 2, figsize=(11, 4.2), dpi=150)
        for c, hh in hist.items():
            if c.startswith('parallel'):
                continue
            axs[0].plot([t for t, _ in hh['dp']], [v for _, v in hh['dp']], lw=1.2, label=c)
            axs[1].plot([t for t, v in hh['xr'] if v], [v / a.h for t, v in hh['xr'] if v], marker='o', label=c)
        axs[0].set(xlabel='SIMPLE iteration', ylabel='p_inlet - p_outlet [m2/s2]', title='pressure drop (kinematic)')
        axs[1].set(xlabel='SIMPLE iteration (write times)', ylabel='x_r / h', title='reattachment point on the lower wall')
        for ax in axs:
            ax.grid(alpha=0.3)
            ax.legend(fontsize=8)
        fig.tight_layout()
        fig.savefig(out_png)
        print('Wrote', out_png)


if __name__ == '__main__':
    main()

To test the scripts without a run, make_synthetic_campaign.py writes logs and monitor files with OpenFOAM 14 names and format from analytic models (exponential decays, a moving reattachment point); it is also the generator used by the foamRun placeholders in the tests on this page. Synthetic data, not OpenFOAM results.

#!/usr/bin/env python3
"""SYNTHETIC test data for the Guide 19 parsers - NOT OpenFOAM results.

  log CASE    print to stdout a foamRun-style log for a prepared case,
              using endTime, div(phi,U), the U relaxation factor and
              numberOfSubdomains of that case to vary an imposed residual
              decay and an imposed cost per iteration;
  post CASE   write CASE/postProcessing files with the names and layout
              of the Guide 19 observables (add_observables.py): flow rates,
              mean pressures, probes and the wall-shear-stress graph.

All numbers come from analytic toy models (exponential decays, a moving
reattachment point); they only exercise residual_history.py,
observables.py and benchmark_stats.py. Every file states SYNTHETIC.
"""
import math, os, re, sys


def rd(case, rel):
    p = os.path.join(case, rel)
    return open(p, encoding='utf-8').read() if os.path.isfile(p) else ''


def case_params(case):
    cd, fs, fv, dp = (rd(case, 'system/' + f) for f in ('controlDict', 'fvSchemes', 'fvSolution', 'decomposeParDict'))
    steps = int(float(re.search(r'(?m)^\s*endTime\s+([0-9.eE+-]+)', cd).group(1)))
    scheme = re.search(r'div\(phi,U\)\s+([^;]+);', fs).group(1)
    relax = float(re.search(r'(?m)^\s*U\s+([0-9.]+);', fv).group(1))
    m = re.search(r'numberOfSubdomains\s+(\d+)', dp)
    n = int(m.group(1)) if m else 1
    tau = 60.0
    if 'upwind;' in scheme + ';' and 'linearUpwind' not in scheme:
        tau = 42.0
    elif 'limitedLinearV' in scheme:
        tau = 52.0
    tau *= 0.9 / relax * (1.0 if relax >= 0.9 else 1.25)
    return steps, scheme, relax, n, tau


def log(case):
    steps, scheme, relax, n, tau = case_params(case)
    cost = {1: 0.0200, 2: 0.0118, 4: 0.0079, 8: 0.0068}.get(n, 0.02 / n)
    w = sys.stdout.write
    w('/*---------------------------------------------------------------------------*\\\n'
      '  SYNTHETIC foamRun log written by make_synthetic_campaign.py - NOT OpenFOAM output\n'
      '\\*---------------------------------------------------------------------------*/\n')
    w(f'Build  : 14 (synthetic)\nExec   : foamRun{" -parallel" if n > 1 else ""}\nnProcs : {n}\n\n'
      'Create time\n\nCreate mesh for time = 0\n\nSelecting solver incompressibleFluid\n\nStarting time loop\n\n')
    wall = 0.6 * n ** 0.5
    cum = 0.0
    for it in range(1, steps + 1):
        e = math.exp(-it / tau)
        osc = 1 + 0.15 * math.sin(it / 7.0)
        rUx = min(1.0, 0.9 * e * osc + 2e-6)
        rUy = min(1.0, 1.0 * e * osc + 3e-6)
        rp = min(1.0, 0.8 * math.exp(-it / (1.3 * tau)) * osc + 5e-6)
        re_ = min(1.0, 0.3 * e + 1e-6)
        rk = min(1.0, 0.6 * e + 1e-6)
        nUx = 1 if it > 20 else 2
        nP = max(2, int(17 * math.exp(-it / 400.0)) + (it % 3))
        cont = 1.2 * math.exp(-it / (1.3 * tau)) + 1e-8
        cum += 0.4 * cont
        wall += cost * (1 + 0.05 * math.sin(it))
        w(f'Time = {it}s\n\n')
        w(f'smoothSolver:  Solving for Ux, Initial residual = {rUx:.6g}, Final residual = {0.05 * rUx:.6g}, No Iterations {nUx}\n')
        w(f'smoothSolver:  Solving for Uy, Initial residual = {rUy:.6g}, Final residual = {0.05 * rUy:.6g}, No Iterations {nUx + 1}\n')
        w(f'GAMG:  Solving for p, Initial residual = {rp:.6g}, Final residual = {0.08 * rp:.6g}, No Iterations {nP}\n')
        w(f'time step continuity errors : sum local = {cont:.6g}, global = {0.1 * cont:.6g}, cumulative = {cum:.6g}\n')
        w(f'smoothSolver:  Solving for epsilon, Initial residual = {re_:.6g}, Final residual = {0.05 * re_:.6g}, No Iterations 3\n')
        w(f'smoothSolver:  Solving for k, Initial residual = {rk:.6g}, Final residual = {0.05 * rk:.6g}, No Iterations 3\n')
        w(f'ExecutionTime = {0.97 * wall:.6g} s  ClockTime = {int(wall)} s\n\n')
    w('End\n\n')
    if os.environ.get('SYNTH_SLEEP'):
        import time
        time.sleep(float(os.environ['SYNTH_SLEEP']) * cost * steps)   # imitate the run time (test harness only)


def post(case):
    steps, scheme, relax, n, tau = case_params(case)
    pp = os.path.join(case, 'postProcessing')
    os.makedirs(pp, exist_ok=False)
    tag = '# SYNTHETIC data written by make_synthetic_campaign.py - NOT OpenFOAM output\n'
    xr_final = {42.0: 0.150, 52.0: 0.163}.get(round(tau / (0.9 / relax * (1.0 if relax >= 0.9 else 1.25)), 1), 0.168)

    def series(name, typ, cols, fun):
        d = os.path.join(pp, name, '0')
        os.makedirs(d)
        with open(os.path.join(d, typ + '.dat'), 'w') as f:
            f.write(tag + '# Time        \t' + '\t'.join(cols) + '\n')
            for it in range(1, steps + 1):
                f.write(f'{it}\t' + '\t'.join(f'{v:.8e}' for v in fun(it)) + '\n')

    Q = 2.54e-4
    series('inletFlowRate', 'surfaceFieldValue', ['sum(phi)'], lambda it: [-Q])
    series('outletFlowRate', 'surfaceFieldValue', ['sum(phi)'], lambda it: [Q * (1 + 0.3 * math.exp(-it / (1.3 * tau)))])
    series('inletPressure', 'surfaceFieldValue', ['areaAverage(p)'], lambda it: [9.6 + 4.0 * math.exp(-it / tau) * math.cos(it / 9.0)])
    series('outletPressure', 'surfaceFieldValue', ['areaAverage(p)'], lambda it: [0.0])
    d = os.path.join(pp, 'probes', '0')
    os.makedirs(d)
    pts = [(0.05, -0.0125, 0.0), (0.15, -0.02, 0.0), (0.25, 0.0, 0.0)]
    with open(os.path.join(d, 'U'), 'w') as f:
        f.write(tag)
        for i, p in enumerate(pts):
            f.write(f'# Probe {i} ({p[0]:g} {p[1]:g} {p[2]:g})\n')
        f.write('#       Probe' + ''.join(f'{i:>27d}' for i in range(len(pts))) + '\n#        Time\n')
        for it in range(1, steps + 1):
            e = math.exp(-it / tau)
            vals = [(-1.6 + 1.5 * e, 0.1, 0), (-0.4 + 2.0 * e, 0.05, 0), (6.1 - 3.0 * e, 0.0, 0)]
            f.write(f'{it:>13d}' + ''.join(f'    ({u:.6g} {v:.6g} {z:g})' for u, v, z in vals) + '\n')
    for it in range(100, steps + 1, 100):
        d = os.path.join(pp, 'lowerWallShearStress', str(it))
        os.makedirs(d)
        xr = xr_final - 0.08 * math.exp(-it / (1.5 * tau))
        with open(os.path.join(d, 'patchCutLayerAverage.xy'), 'w') as f:
            f.write('#     distance            x            y            z    wallShe...    wallShe...    wallShe...\n')
            for i in range(120):
                x = -0.0206 + (0.29 + 0.0206) * i / 119
                if x < 0:
                    tau_x = -0.30
                elif x < 0.012:
                    tau_x = -0.02 * math.sin(math.pi * x / 0.012)          # corner eddy
                else:
                    tau_x = 0.25 * math.sin(math.pi * (x - 0.012) / (xr - 0.012)) if x < xr else -0.35 * (1 - math.exp(-(x - xr) / 0.02))
                y = 0.0 if x < 0 else -0.0254
                f.write(f'{x + 0.0206:13.6g}{x:13.6g}{y:13.6g}{0:13.6g}{tau_x:13.6g}{0:13.6g}{0:13.6g}\n')


if __name__ == '__main__':
    if len(sys.argv) != 3 or sys.argv[1] not in ('log', 'post'):
        sys.exit('usage: make_synthetic_campaign.py log|post CASE')
    (log if sys.argv[1] == 'log' else post)(sys.argv[2])

$ for c in reference upwind limited relax_06; do bash run_case.sh ~/OpenFOAM/guide19_run02/$c 1; done
[exit 0]
$ ls ~/OpenFOAM/guide19_run02/reference/postProcessing ~/OpenFOAM/guide19_run02/reference/postProcessing/lowerWallShearStress
/home/user/OpenFOAM/guide19_run02/reference/postProcessing:
inletFlowRate
inletPressure
lowerWallShearStress
outletFlowRate
outletPressure
probes

/home/user/OpenFOAM/guide19_run02/reference/postProcessing/lowerWallShearStress:
100
200
300
400
500
600
[exit 0]
$ head -4 ~/OpenFOAM/guide19_run02/reference/postProcessing/lowerWallShearStress/600/patchCutLayerAverage.xy
#     distance            x            y            z    wallShe...    wallShe...    wallShe...
            0      -0.0206            0            0         -0.3            0            0
   0.00261008   -0.0179899            0            0         -0.3            0            0
   0.00522017   -0.0153798            0            0         -0.3            0            0
[exit 0]
$ python3 observables.py ~/OpenFOAM/guide19_run02 --cases parallel_2
parallel_2: observables not found (were they installed with add_observables.py before the run?)
no observables found
[exit 1]
$ python3 observables.py ~/OpenFOAM/guide19_run02 --cases reference upwind limited relax_06 --rho 1.2 --plot
variations are (max - min)/|mean| over the last 100 iterations
case           it   Q_in m3/s    imbal.  dp m2/s2      var    dp Pa     Ux0     Ux1     Ux2  x_r/h  x_r/h at write times
reference     600 -2.5400e-04  1.37e-04    9.5999  1.4e-04   11.520  -1.600  -0.400   6.100  6.596  100:5.567 200:6.261 300:6.491 400:6.559 500:6.592 600:6.596
upwind        600 -2.5400e-04  5.07e-06    9.6000  3.5e-06   11.520  -1.600  -0.400   6.100  5.885  100:5.256 200:5.769 300:5.872 400:5.882 500:5.884 600:5.885
limited       600 -2.5400e-04  4.19e-05    9.6000  3.8e-05   11.520  -1.600  -0.400   6.100  6.396  100:5.524 200:6.171 300:6.321 400:6.388 500:6.394 600:6.396
relax_06      600 -2.5400e-04  4.96e-03    9.5851  8.1e-03   11.502  -1.593  -0.390   6.086  6.501  100:4.857 200:5.635 300:6.077 400:6.297 500:6.422 600:6.501
Wrote /home/user/OpenFOAM/guide19_run02/observables_summary.csv
Wrote /home/user/OpenFOAM/guide19_run02/observables.png
[exit 0]

• The tutorial geometry is that of the experiment of R. W. Pitz, J. W. Daily, “Combustion in a turbulent mixing layer formed at a rearward-facing step”, AIAA Journal 21, 1565–1570, 1983, doi:10.2514/3.8290. A quantitative comparison with the measurements requires their inlet conditions, different from the 10 m/s of the tutorial: here x_r/h serves as a numerical observable to compare schemes and convergence.
• In backward-facing step flows the standard k–ε model tends to underpredict the recirculation length compared with measurements (a classic reference for the reattachment of a shear layer is D. M. Driver, H. L. Seegmiller, “Features of a reattaching turbulent shear layer in divergent channel flow”, AIAA Journal 23, 163–171, 1985, doi:10.2514/3.8890); the numerical diffusion of upwind acts in the same direction. The two causes add up: a shorter x_r with upwind is not a sign of higher accuracy.
• The pressure of incompressibleFluid is kinematic (p/ρ, m²/s²): to obtain Pa multiply by the density of the fluid being represented (--rho 1.2 for air); the case only defines ν = 10⁻⁵ m²/s.

The thresholds are a starting point to adapt to how the result will be used; they must be written in the protocol before seeing the data.

• “There are not enough slots available in the system”: Open MPI sees fewer cores than the requested processes; reduce -np or check with nproc and lscpu.
• Open MPI refuses to start as the root user: use a normal user, as the package also recommends.
• MPI count mismatch from run_case.sh: the number of processes passed does not match the case numberOfSubdomains (test in section 12).
• analyze_campaign.py stops with “Refusing to overwrite results_summary.csv”: rename the previous file (section 14).
• Scripts edited on Windows that fail with $'\r': command not found: convert the line endings to LF (Notepad++: Edit › EOL Conversion › Unix).
• Very different times between two identical repetitions: other processes on the server, frequency scaling, processes not pinned to cores; timing.csv records the load average before each run.

Campaign: guide19_run01            Date: ____    Operator: ____
OpenFOAM: Foundation 14, commit/package: ____   (foamVersion)
Server: ____ physical cores, CPU ____, RAM ____, filesystem ____
Case: pitzDailySteady, kEpsilon, nu 1e-05, inlet U 10 m/s
Mesh: resources/blockMesh/pitzDaily, 12225 cells, checkMesh: ____
Variants: campaign_manifest.csv (attached), sha256sum of system/* (attached)
Iterations: 600, residualControl disabled
Observables defined beforehand: dp [m2/s2], x_r/h, flow balance, probe Ux
Acceptance thresholds: ____ (section 16)
Timings: benchmark_repeats.py, N = __ repetitions, median, load average before
Results: results_summary.csv, residual_history.csv, observables_summary.csv,
         timing.csv, plots (attached)
Conclusions and limits: ____

The fingerprints of the configuration files are obtained with sha256sum ~/OpenFOAM/guide19_run01/*/system/* > config_sha256.txt: they document that the cases did not change between preparation and run.

Finite volumes and pressure–velocity coupling

• J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, Springer, 2020, doi:10.1007/978-3-319-99693-6.
• F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, Springer, 2016, doi:10.1007/978-3-319-16874-6.
• S. V. Patankar, D. B. Spalding, Int. J. Heat Mass Transfer 15, 1787–1806, 1972, doi:10.1016/0017-9310(72)90054-3.
• J. P. Van Doormaal, G. D. Raithby, Numer. Heat Transfer 7, 147–163, 1984, doi:10.1080/01495728408961817.
• C. M. Rhie, W. L. Chow, AIAA J. 21, 1525–1532, 1983, doi:10.2514/3.8284.

Convection schemes

• R. F. Warming, R. M. Beam, AIAA J. 14, 1241–1249, 1976, doi:10.2514/3.61457.
• P. K. Sweby, SIAM J. Numer. Anal. 21, 995–1011, 1984, doi:10.1137/0721062.

Linear solvers and parallel computing

• K. Stüben, J. Comput. Appl. Math. 128, 281–309, 2001, doi:10.1016/S0377-0427(00)00516-1.
• F. Pellegrini, J. Roman, LNCS 1067, 493–498, 1996, doi:10.1007/3-540-61142-8_588.
• G. M. Amdahl, AFIPS 1967, 483–485, doi:10.1145/1465482.1465560.
• J. L. Gustafson, Commun. ACM 31, 532–533, 1988, doi:10.1145/42411.42415.
• A. H. Karp, H. P. Flatt, Commun. ACM 33, 539–543, 1990, doi:10.1145/78607.78614.

Turbulence, case and verification

• B. E. Launder, D. B. Spalding, Comput. Methods Appl. Mech. Eng. 3, 269–289, 1974, doi:10.1016/0045-7825(74)90029-2.
• R. W. Pitz, J. W. Daily, AIAA J. 21, 1565–1570, 1983, doi:10.2514/3.8290.
• D. M. Driver, H. L. Seegmiller, AIAA J. 23, 163–171, 1985, doi:10.2514/3.8890.
• P. J. Roache, J. Fluids Eng. 116, 405–413, 1994, doi:10.1115/1.2910291.
• I. B. Celik et al., J. Fluids Eng. 130, 078001, 2008, doi:10.1115/1.2960953.

OpenFOAM 14 sources cited in the boxes

• applications/modules/incompressibleFluid/momentumPredictor.C
• applications/modules/incompressibleFluid/correctPressure.C
• applications/solvers/foamRun/foamRun.C
• src/finiteVolume/…/boundedConvectionScheme.C
• src/finiteVolume/…/linearUpwind.C
• src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
• src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
• src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C
• src/finiteVolume/…/singleRegionConvergenceControl.C
• src/OSspecific/POSIX/cpuTime/cpuTime.C
• src/OpenFOAM/global/clock/clock.C
• src/sampling/sampledSet/writers/setWriter.C