OpenFOAM: metodi numerici avanzati, convergenza e prestazioni

Prima di confrontare schemi e tempi conviene sapere su che cosa si lavora. extra/case_estimates.py legge il dizionario resources/blockMesh/pitzDaily (lo stesso usato da run_case.sh), applica le regole di gradazione di blockMesh, comprese le sezioni multiple come $posY, e legge le condizioni d’ingresso del caso. I numeri principali:

Lo script e la sua esecuzione sul caso preparato sono nel riquadro della sezione 7 (qualità della mesh); il valore di ω all’ingresso del tutorial, 440,15 s⁻¹, coincide con ε/(C_μk): le condizioni d’ingresso dei diversi modelli di turbolenza sono coerenti fra loro.

Contenuto dei file (senza intestazioni), da confrontare con la copia sul server. Il caso completo è nel repository della versione 14. system/blockMeshDict è una copia di tutorials/resources/blockMesh/pitzDaily (lo dice il commento nel file): Allrun e run_case.sh usano la versione delle risorse con -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 (usati solo da altri modelli di turbolenza)

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)
        );
    }
);

Dettagli utili: il modello kEpsilon è quello attivo (il commento elenca i modelli provati con lo stesso caso); residualControl fermerebbe il calcolo prima delle 2000 iterazioni di endTime; consistent yes attiva la variante SIMPLEC (riquadro seguente); wallDist { method meshWave; } serve ai modelli che usano la distanza dalla parete; la funzione streamlinesLine scrive le linee di corrente a ogni salvataggio.

Il modulo incompressibleFluid costruisce l’equazione della quantità di moto in momentumPredictor.C, la rilassa e la risolve con il gradiente di pressione dell’iterazione precedente:

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, con consistent yes, il coefficiente dell’equazione di pressione non è 1/A ma quello «consistente» di 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);

• Con SIMPLEC la correzione di velocità tiene conto dei coefficienti delle celle vicine (H1): per questo il tutorial non sotto-rilassa la pressione (nessuna voce fields { p … }, quindi p.relax() non agisce) e usa 0,9 per tutte le equazioni, mentre il SIMPLE classico di 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) richiede tipicamente 0,3 su p e 0,7 su U.
• L’interpolazione dei flussi phiHbyA dai valori di cella segue l’idea di C. M. Rhie, W. L. Chow (AIAA Journal 21, 1525–1532, 1983, doi:10.2514/3.8284), che evita le oscillazioni a scacchiera di pressione su griglie co-locate.
• Il modello di turbolenza (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) è risolto dopo la pressione: nel log compaiono, per ogni iterazione, Ux, Uy, p, la riga di continuità, epsilon e k.

• Errore iterativo. È la differenza fra la soluzione all’iterazione corrente e quella che il sistema discreto avrebbe a convergenza. Non si misura con il residuo di un singolo solutore lineare, ma osservando come cambiano le grandezze d’interesse: se la caduta di pressione e il punto di riattacco variano ancora dello 0,1 % ogni 100 iterazioni, l’errore iterativo è di quell’ordine (sezione 15, observables.py).
• Errore di discretizzazione. Dipende da mesh e schemi; a convergenza iterativa si stima con almeno tre mesh in rapporto costante, l’ordine osservato e il 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 (Guida 4). Il confronto fra schemi su una sola mesh mostra quanto l’errore di discretizzazione dipenda dallo schema, non quale schema sia «giusto».
• Regola pratica. L’errore iterativo deve essere almeno un ordine di grandezza più piccolo delle differenze che si vogliono attribuire agli schemi: altrimenti il confronto misura soltanto quanto ogni caso sia lontano dalla convergenza dopo 600 iterazioni.
• Testi di riferimento: J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, 4ª 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 (con esempi in OpenFOAM).

bounded. In boundedConvectionScheme.C la variante bounded sottrae dal termine convettivo il prodotto della divergenza del flusso per il campo:

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

Cioè discretizza ∇·(φU) − (∇·φ)U. A convergenza ∇·φ = 0 e il termine sparisce; durante le iterazioni di un caso stazionario, quando la continuità non è ancora soddisfatta, evita che l’errore di massa agisca come una sorgente spuria di quantità di moto. È un accorgimento per la stabilità del ciclo SIMPLE, indipendente dallo schema di interpolazione che segue.

linearUpwind. In linearUpwind.C il valore sulla faccia è quello della cella a monte più una correzione esplicita con il gradiente (qui grad(U), calcolato con Gauss linear):

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

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

• È l’estrapolazione del secondo ordine di 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. La parte upwind è implicita, la correzione esplicita: non è limitata, quindi può creare piccole sovraelongazioni e sottoelongazioni dove i gradienti cambiano bruscamente, come nello strato di taglio dopo lo scalino.
• upwind. Solo il valore a monte: limitato e robusto, ma del primo ordine, con errore di troncamento equivalente a una diffusività |U|Δx/2 nella direzione del flusso; nel caso vale da 1 a 30 volte la viscosità turbolenta d’ingresso (sezione 1). Accorcia tipicamente le zone di ricircolo.
• limitedLinearV 1. Interpolazione lineare limitata con il limitatore TVD di tipo limitedLinear (coefficiente 1 = limitazione più forte, verso upwind); la «V» indica la versione per vettori, che calcola il limitatore una sola volta nella direzione di massima variazione del vettore invece che componente per componente, preservando la direzione della velocità. Teoria dei limitatori: 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.

• Il punto di riattacco sulla parete inferiore è l’osservabile più sensibile alla diffusione numerica: observables.py lo ricava dallo sforzo a parete a ogni salvataggio (sezione 15). Con upwind ci si attende un ricircolo più corto.
• La caduta di pressione fra ingresso e uscita è una grandezza integrale, meno sensibile; i profili di Ux in sezioni a x = 2h, 4h, 6h (Plot Over Line in ParaView) mostrano dove gli schemi differiscono.
• Solo lo schema di U cambia: k ed ε restano bounded Gauss limitedLinear 1 in tutti i casi, come nel tutorial.

Il residuo stampato nel log è normalizzato in lduMatrixSolver.C: la somma dei moduli di b − Ax è divisa per un fattore che dipende dalla soluzione stessa, così il valore non cambia con la scala delle grandezze:

// --- 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_;

Il criterio di arresto del solutore lineare, in SolverPerformance.C, è «residuo finale < tolerance» oppure «residuo finale < relTol × residuo iniziale»:

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

• Con relTol 0.1 ogni soluzione di p riduce il residuo di un fattore 10 e si ferma: nelle prime iterazioni bastano pochi cicli GAMG, verso la convergenza il residuo iniziale è già piccolo e il lavoro per iterazione cambia. Per questo il numero di iterazioni lineari (No Iterations) va registrato: residual_history.py le somma per caso.
• pcorr (relTol 0) serve solo alla correzione iniziale dei flussi; GAMG è un multigrid algebrico (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) con levigatore GaussSeidel; per U, k ed ε il tutorial usa smoothSolver con symGaussSeidel, adatto a matrici non simmetriche e diagonalmente dominanti.

residualControl. In singleRegionConvergenceControl.C il ciclo si ferma quando, per tutti i campi elencati, il residuo iniziale della prima soluzione dell’iterazione è sotto la soglia; per i vettori conta la componente più grande (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());

Lo script legge log.foamRun di ogni caso iterazione per iterazione e calcola: il residuo iniziale della prima soluzione di ogni campo (quello di residualControl), l’iterazione in cui i criteri del tutorial (p 10⁻², U 10⁻³, k ed ε 10⁻³, modificabili con --criteria) sarebbero stati soddisfatti, le iterazioni lineari totali e medie di p, l’errore di continuità finale e il costo medio per iterazione dagli ExecutionTime successivi (escluse le prime 10). Scrive residual_history.csv e, con --plot, il grafico dei residui di Ux e p.

#!/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()

Prova sui log sintetici scritti dal generatore di prova del riquadro della sezione 15 (modelli esponenziali, non risultati OpenFOAM): il formato delle righe è quello del log di foamRun in OpenFOAM 14 (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

Per le equazioni (relaxationFactors/equations) OpenFOAM applica un rilassamento implicito in fvMatrix.C: rende la matrice diagonalmente dominante, divide la diagonale per α e aggiunge al termine noto la differenza, così la soluzione rilassata resta coerente con l’equazione originale a convergenza:

// 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();

• A convergenza U non cambia più fra due iterazioni e il termine aggiunto si annulla: il rilassamento implicito cambia il percorso, non la soluzione finale. Una differenza fra reference e relax_06 dopo 600 iterazioni misura quindi soltanto la diversa distanza dalla convergenza.
• Con SIMPLEC il coefficiente consistente rAtU è costruito dalla matrice rilassata: abbassare α da 0,9 a 0,6 rallenta anche la correzione di pressione. È atteso un numero di iterazioni nettamente più alto per raggiungere gli stessi residui; residual_history.py lo misura con la colonna «stop@».
• Il commento del tutorial («0.9 is more stable but 0.95 more convergent») indica il verso opposto: con SIMPLEC si può salire verso 0,95 se il caso è stabile.

Lo script ricostruisce le dimensioni delle celle con le regole di blockMesh: per ogni spigolo del blocco legge la gradazione (simpleGrading o le 12 voci di edgeGrading, anche con sezioni multiple «frazione di lunghezza, frazione di celle, rapporto»), calcola la progressione geometrica delle celle e riporta minimi, massimi e rapporto d’aspetto massimo per blocco; poi le grandezze fisiche d’ingresso e la diffusione numerica dell’upwind. Solo libreria standard di Python.

#!/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

• Il totale, 12 225 celle, coincide con la somma dei blocchi del dizionario (540 + 4860 + 5400 + 675 + 750): è il numero da ritrovare nel log.checkMesh («cells: 12225»).
• Il rapporto d’aspetto riportato è un limite superiore (cella più lunga su cella più bassa del blocco); quello effettivo, cella per cella, è nella riga «Max aspect ratio» di checkMesh. Il blocco 4, rastremato e con la gradazione $posYR, è quello con le celle più allungate.
• Dal log.checkMesh registrare nel protocollo: numero di celle, «Max non-orthogonality», «Max skewness», «Max aspect ratio» e l’eventuale riga «Failed … mesh checks». Le facce dei blocchi rastremati non sono ortogonali: la correzione Gauss linear corrected dei laplaciani compensa la non ortogonalità, con nNonOrthogonalCorrectors 0.

run_case.sh genera sempre la mesh dal dizionario delle risorse (blockMesh -dict $FOAM_TUTORIALS/resources/blockMesh/pitzDaily), quindi una mesh raffinata va preparata in una copia nuova e lanciata a mano, usando il system/blockMeshDict del caso:

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

Il rapporto di gradazione è definito fra l’ultima e la prima cella del tratto: raddoppiando le celle con lo stesso rapporto le celle estreme si dimezzano e la progressione resta geometricamente simile. Con 4 volte le celle servono più iterazioni SIMPLE per convergere: 600 iterazioni non bastano in generale, e il confronto con la mesh base va fatto a convergenza iterativa (sezione 3).

• In OpenFOAM 14 foamRun stampa Time = seguito dal nome del tempo con l’unità (foamRun.C, runTime.userTimeName()): nel log si legge «Time = 600s». Le espressioni regolari degli script accettano sia «600» sia «600s».
• Con writeInterval 100 e 600 iterazioni si scrivono le cartelle 100…600: sono iterazioni, non secondi. foamListTimes le elenca; nei casi paralleli sono dentro processorN/ finché reconstructPar -latestTime non ricostruisce l’ultima.
• Il tutorial originale, con residualControl attivo, si ferma da solo quando i residui scendono sotto le soglie e scrive quell’ultima iterazione; nelle copie della campagna il controllo è tolto e tutte le varianti fanno esattamente lo stesso numero di iterazioni.

• scotch partiziona il grafo delle celle minimizzando le facce di confine fra processori, senza coefficienti geometrici (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). Il manuale di OpenFOAM 14 (sezione 3.4) descrive anche simple e hierarchical, geometrici.
• In log.decomposePar leggere, per ogni processore, il numero di celle e di facce di processore: la loro somma rispetto alle facce interne indica quanta comunicazione serve a ogni iterazione. Con circa 3000 celle per processo (4 processi) una parte non trascurabile delle celle confina con un altro sottodominio.
• decomposePar -cellProc scrive il campo cellProc (numero del processore di ogni cella), da visualizzare in ParaView per vedere i sottodomini.
• Usare solo core fisici: lscpu mostra «Core(s) per socket» e «Thread(s) per core». Con Open MPI, --bind-to core fissa ogni processo a un core e riduce la variabilità fra ripetizioni; «There are not enough slots available in the system» significa che si chiedono più processi dei core disponibili (--oversubscribe lo forza, ma il tempo misurato non è più uno speedup significativo).
• Il risultato parallelo deve coincidere con quello seriale entro la tolleranza dei solutori: i residui iniziali del log possono differire un poco perché GAMG e i levigatori lavorano per sottodominio; caduta di pressione e punto di riattacco devono restare uguali alle cifre significative scelte.

ParaView legge direttamente i casi decomposti: nel lettore OpenFOAM scegliere «Decomposed Case» in Case Type. È utile per controllare i sottodomini (con cellProc) e i tempi intermedi, che reconstructPar -latestTime non ricostruisce. Per ricostruire tutto: reconstructPar senza opzioni, in una copia dopo il benchmark, non nella cartella usata per misurare i tempi.

La riga stampata alla fine di ogni iterazione viene da foamRun.C:

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

ExecutionTime è il tempo CPU del processo: cpuTime.C usa clock() della libreria C. Info scrive solo dal processo principale, quindi nel caso parallelo è il tempo CPU del rank 0, non la somma dei processi (la guida è stata corretta su questo punto). Con MPI può avvicinarsi al tempo trascorso, perché l’attesa attiva dei messaggi consuma 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 è un time_t: secondi interi dall’avvio di foamRun (clock.C). Su una corsa di 10 s l’errore di arrotondamento arriva al 10 %, e uno speedup calcolato da due valori interi può variare di molto da una prova all’altra. analyze_campaign.py usa proprio questo valore; per misure accurate extra/benchmark_repeats.py misura il tempo trascorso di ogni fase con time.perf_counter() e ripete le corse.

Le leggi da usare per leggere i numeri. Legge di Amdahl, S(n) = 1/(f + (1 − f)/n) con f frazione non parallelizzabile (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); scalabilità a carico crescente (J. L. Gustafson, «Reevaluating Amdahl’s law», Communications of the ACM 31, 532–533, 1988, doi:10.1145/42411.42415): per un problema di dimensione fissa vale Amdahl, per mesh che crescono con i processi conta l’efficienza «weak scaling»; metrica di Karp–Flatt, e(n) = (1/S − 1/n)/(1 − 1/n), frazione seriale ricavata dalle misure (A. H. Karp, H. P. Flatt, «Measuring parallel processor performance», Communications of the ACM 33, 539–543, 1990, doi:10.1145/78607.78614): se e cresce con n il limite non è il codice seriale ma il costo della comunicazione.

benchmark_repeats.py prende un caso preparato e mai avviato (per esempio parallel_4), lo copia N volte in una cartella nuova e ripete la stessa sequenza di run_case.sh misurando il tempo trascorso di ogni fase (blockMesh, checkMesh, decomposePar, solutore, reconstructPar). In timing.csv registra anche l’ultimo ExecutionTime/ClockTime del log, il numero di iterazioni, il carico medio del server prima di ogni corsa, il numero di core e il modello di CPU. Il caso sorgente non viene toccato; un errore ferma le ripetizioni.

#!/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 unisce i timing.csv, controlla che tutte le corse abbiano lo stesso numero di iterazioni, e per ogni numero di processi riporta mediana, minimo–massimo, dispersione, millisecondi per iterazione, ClockTime del log, tempo complessivo con tutte le fasi, speedup, efficienza, metrica di Karp–Flatt, e la frazione seriale di Amdahl adattata ai dati con le previsioni per 8 e 16 processi.

#!/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()

Prova con i comandi OpenFOAM sostituiti da segnaposto che scrivono log sintetici e attendono un tempo proporzionale a un costo per iterazione imposto (sezione 15): i numeri mostrano il funzionamento degli script e la differenza fra il ClockTime intero del log e il tempo misurato, non le prestazioni di OpenFOAM.

$ 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()

• Copia solo 0/, constant/ e system/ del tutorial; ogni sostituzione deve trovare esattamente una riga (div(phi,U), endTime, writeInterval, il fattore di U a 0,9) e un solo blocco residualControl, altrimenti si ferma e cancella soltanto la cartella che ha creato.
• Il diff con il tutorial mostra le sole modifiche attese, più una conseguenza estetica: l’espressione ^\s* con re.M comprende anche le righe vuote precedenti, che vengono eliminate davanti a endTime e writeInterval. Non cambia il contenuto dei dizionari.
• writeInterval resta 100 (come nel tutorial): con --steps non multiplo di 100 l’ultima iterazione può non essere scritta su disco; usare multipli di 100.

Prova con il tutorial della versione 14 (FOAM_TUTORIALS impostato come fa etc/bashrc):

$ 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]

Le grandezze fisiche vanno definite prima della prova (sezione 13). Lo script aggiunge al system/functions del tutorial, che conserva com’è (copia in functions.tutorial), le funzioni standard di etc/caseDicts/functions con la sintassi degli argomenti di #includeFunc: portate su ingresso e uscita, pressione media sulle due sezioni, tre sonde di U e p, lo sforzo a parete e il suo profilo lungo x sulla parete inferiore (graphPatchCutLayerAverage). Rifiuta casi già avviati, con mesh o con un functions diverso dal tutorial; controlla tutti i casi prima di modificarne uno.

#!/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()

Il file risultante (extra/functions_guide19 nell’archivio):

#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)

Prova su una campagna nuova (il primo comando rifiuta l’elenco perché contiene un caso già calcolato e non modifica nulla):

$ 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

• Codici d’uscita: 64 per errori d’uso (cartella, numero di processi diverso da 1, 2, 4 o da quello del decomposeParDict), 69 per comandi mancanti, 73 se il caso è già stato avviato, ha una mesh o cartelle processor*.
• RUN_STARTED riceve all’uscita il codice e l’ora di fine anche in caso di errore; RUN_SUCCESS è scritto solo se tutti i comandi riescono (set -Eeuo pipefail).
• Come nella guida precedente, checkMesh termina sempre con codice 0 anche quando i controlli falliscono (checkMesh.C): controllare a mano grep -E "Mesh OK|Failed" log.checkMesh e che log.foamRun termini con «End».

Prova con i comandi OpenFOAM sostituiti da segnaposto (stessa sequenza di file e di controlli; il secondo lancio di reference e parallel_2 con 4 processi mostrano i rifiuti):

$ 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]

• Prima della prova scrivere quali grandezze decidono il confronto e con quale tolleranza, per esempio: caduta di pressione con 4 cifre significative, x_r/h con 2 decimali, bilancio di portata entro 10⁻⁴. Installare i monitor con add_observables.py su entrambe le copie prima del calcolo.
• Dopo il calcolo: python3 residual_history.py ~/OpenFOAM/guide19_run01 --cases reference relax_06 --plot (residui alla stessa iterazione e colonna «stop@») e python3 observables.py ~/OpenFOAM/guide19_run01 --cases reference relax_06 --plot (variazione delle osservabili nelle ultime 100 iterazioni e x_r a ogni salvataggio).
• Se in relax_06 le osservabili variano ancora, preparare una campagna nuova con più iterazioni, per esempio python3 prepare_campaign.py --out ~/OpenFOAM/guide19_run02 --steps 2000, invece di riprendere i casi già usati. A convergenza le due varianti devono dare gli stessi valori (sezione 6): una differenza persistente indica un problema di convergenza, non un effetto fisico del rilassamento.

#!/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()

• Legge l’ultimo Time, l’ultima coppia ExecutionTime/ClockTime, gli ultimi residui iniziali di Ux e p e l’ultimo errore locale di continuità; scrive results_summary.csv e disegna lo speedup solo se i tre casi del confronto MPI sono riusciti e hanno raggiunto le iterazioni richieste.
• Lo speedup usa ClockTime, che OpenFOAM 14 stampa in secondi interi (sezione 10), e una sola corsa per caso: va letto come un’indicazione; per la mediana e una risoluzione al millisecondo usare benchmark_repeats.py e benchmark_stats.py.
• Il riepilogo viene scritto prima del grafico e non viene mai sovrascritto: per avere anche speedup.png usare --plot già alla prima esecuzione, altrimenti rinominare results_summary.csv e ripetere.

Prova sui log sintetici prodotti dai segnaposto (il ClockTime intero del log sintetico riproduce la risoluzione di quello vero):

$ 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]

Lo script legge i file dei monitor di add_observables.py e riporta, all’ultima iterazione e come variazione relativa nelle ultime --window iterazioni: portate in ingresso e in uscita e loro sbilancio, caduta di pressione cinematica (in Pa con --rho), Ux alle tre sonde, e il punto di riattacco x_r sulla parete inferiore a ogni salvataggio. x_r è l’ultimo cambio di segno della componente x dello sforzo a parete prima della zona con il segno di valle: il calcolo non dipende dalla convenzione di segno di wallShearStress ed esclude il piccolo vortice d’angolo vicino allo scalino. Le colonne del file .xy sono riconosciute dall’intestazione scritta da setWriter.C, che accorcia con «…» i nomi troppo lunghi.

#!/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()

Per provare gli script senza un calcolo, make_synthetic_campaign.py scrive log e file dei monitor con i nomi e il formato di OpenFOAM 14 a partire da modelli analitici (decadimenti esponenziali, un punto di riattacco che si sposta); è anche il generatore usato dai segnaposto di foamRun nelle prove di questa pagina. Dati sintetici, non risultati OpenFOAM.

#!/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]

• La geometria del tutorial è quella dell’esperimento di 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. Il confronto quantitativo con le misure richiede le loro condizioni d’ingresso, diverse dai 10 m/s del tutorial: qui x_r/h serve come osservabile numerica per confrontare schemi e convergenza.
• Nei flussi dietro uno scalino il modello k–ε standard tende a sottostimare la lunghezza di ricircolo rispetto alle misure (un riferimento classico per il riattacco di uno strato di taglio è 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); la diffusione numerica dell’upwind agisce nello stesso verso. Le due cause si sommano: un x_r più corto con upwind non è un’indicazione di maggiore accuratezza.
• La pressione di incompressibleFluid è cinematica (p/ρ, m²/s²): per avere Pa si moltiplica per la densità del fluido che si vuole rappresentare (--rho 1.2 per l’aria); il caso definisce solo ν = 10⁻⁵ m²/s.

Le soglie sono un punto di partenza da adattare all’uso del risultato; vanno scritte nel protocollo prima di vedere i dati.

• «There are not enough slots available in the system»: Open MPI vede meno core dei processi richiesti; ridurre -np o verificare con nproc e lscpu.
• Open MPI rifiuta di partire come utente root: usare un utente normale, come raccomanda anche il pacchetto.
• MPI count mismatch da run_case.sh: il numero di processi passato non coincide con numberOfSubdomains del caso (prova nella sezione 12).
• analyze_campaign.py si ferma con «Refusing to overwrite results_summary.csv»: rinominare il file precedente (sezione 14).
• Script modificati su Windows che falliscono con $'\r': command not found: convertire i fine riga in LF (Notepad++: Modifica › Conversione fine riga › Unix).
• Tempi molto diversi fra due ripetizioni identiche: altri processi sul server, regolazione della frequenza, processi non fissati ai core; timing.csv registra il carico medio prima di ogni corsa.

Campagna: guide19_run01            Data: ____    Operatore: ____
OpenFOAM: Foundation 14, commit/pacchetto: ____   (foamVersion)
Server: ____ core fisici, CPU ____, RAM ____, filesystem ____
Caso: pitzDailySteady, kEpsilon, nu 1e-05, U ingresso 10 m/s
Mesh: resources/blockMesh/pitzDaily, 12225 celle, checkMesh: ____
Varianti: campaign_manifest.csv (allegato), sha256sum di system/* (allegato)
Iterazioni: 600, residualControl disattivato
Osservabili definite prima: dp [m2/s2], x_r/h, bilancio di portata, Ux sonde
Soglie di accettazione: ____ (sezione 16)
Tempi: benchmark_repeats.py, N = __ ripetizioni, mediana, carico medio prima
Risultati: results_summary.csv, residual_history.csv, observables_summary.csv,
           timing.csv, grafici (allegati)
Conclusioni e limiti: ____

Le impronte dei file di configurazione si ottengono con sha256sum ~/OpenFOAM/guide19_run01/*/system/* > config_sha256.txt: documentano che i casi non sono cambiati fra la preparazione e il calcolo.

Volumi finiti e accoppiamento pressione–velocità

• 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.

Schemi convettivi

• 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.

Solutori lineari e parallelo

• 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.

Turbolenza, caso e verifica

• 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.

Sorgenti di OpenFOAM 14 citati nei riquadri

• 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