OpenFOAM: campagne parametriche, automazione e calcolo parallelo

La differenza fra le pressioni medie sulle patch inlet e outlet raccoglie tutto ciò che accade nel dominio: la perdita dovuta al cilindro, l’eventuale attrito sulle pareti laterali, la differenza di flusso di quantità di moto fra due profili di velocità non uniformi e la scia non ancora riassorbita alla sezione di uscita. Per questo cambia anche spostando l’uscita più a valle, a parità di cilindro.

La forza sul solo cilindro si ottiene integrando pressione e sforzi viscosi sulla sua superficie, con il function object forces (o forceCoeffs per i coefficienti adimensionali) applicato alla patch del cilindro. Il nome esatto della patch si legge in constant/polyMesh/boundary dopo blockMesh e mirrorMesh. Nei casi incomprimibili il campo p è cinematico e la forza va riportata in newton con una densità di riferimento dichiarata, esattamente come la Δp di questa guida: le parole chiave della versione 14 sono nella documentazione API di forces.

# nomi e tipi delle patch dopo la generazione della mesh
less constant/polyMesh/boundary

Per un cilindro di diametro D in una corrente di velocità U e viscosità cinematica ν, il parametro che governa il regime è Re = U·D/ν. La viscosità del caso si legge in constant/physicalProperties, il diametro dalla geometria di system/blockMeshDict. Nella campagna la velocità varia di un fattore 4 (da 0,005 a 0,020 m/s), quindi anche Re varia di un fattore 4.

Per un cilindro circolare in corrente indefinita la scia è stazionaria fino a Re ≈ 47; oltre compare il distacco periodico dei vortici. Se l’intervallo della campagna attraversa questa soglia, alcuni casi saranno stazionari e altri periodici: per questi ultimi l’ultimo campione registrato dipende dalla fase dell’oscillazione e il confronto va fatto su medie temporali su un numero intero di periodi. Il confinamento del dominio sposta la soglia, quindi va controllata nei log e nei campioni, non presunta.

# viscosità e geometria del caso
foamDictionary constant/physicalProperties -entry nu -value
less system/blockMeshDict

foamDictionary interpreta il dizionario come lo fa OpenFOAM, quindi è più affidabile di sed o grep su file che contengono macro, commenti e sottodizionari. Le voci annidate si indicano con la barra.

foamDictionary 0/U -keywords                         # voci di primo livello
foamDictionary 0/U -entry Uinlet -value             # valore della velocità
foamDictionary 0/U -entry boundaryField/inlet       # condizione sulla patch inlet
foamDictionary -expand 0/U                          # file con le macro $Uinlet sostituite
foamDictionary system/controlDict -entry endTime -value
foamDictionary system/controlDict -entry functions -keywords

L’ultima riga è il modo più rapido per verificare, dopo la preparazione, che nel controlDict di ogni caso compaiano pIn, pOut, qIn e qOut. Leggere sempre anche l’Allrun del tutorial (cat "$FOAM_TUTORIALS/incompressibleFluid/cylinder/Allrun"): è la descrizione ufficiale della sequenza che gli script replicano.

config/velocities.csv: una riga per caso, intestazione obbligatoria, punto decimale.

case_id,U_m_s
U005,0.005
U010,0.010
U015,0.015
U020,0.020

config/measurement_functions.foam: il blocco functions che lo script accoda al controlDict di ogni caso.

// Tutorial package: quantities measured during the solver run.
// This block is appended to the clean OpenFOAM-14 cylinder controlDict.
functions
{
    pIn
    {
        type surfaceFieldValue;
        libs ("libfieldFunctionObjects.so");
        patch inlet;
        operation areaAverage;
        fields (p);
        writeControl writeTime;
        writeFields false;
        log true;
    }
    pOut
    {
        type surfaceFieldValue;
        libs ("libfieldFunctionObjects.so");
        patch outlet;
        operation areaAverage;
        fields (p);
        writeControl writeTime;
        writeFields false;
        log true;
    }
    qIn
    {
        type surfaceFieldValue;
        libs ("libfieldFunctionObjects.so");
        patch inlet;
        operation sum;
        fields (phi);
        writeControl writeTime;
        writeFields false;
        log true;
    }
    qOut
    {
        type surfaceFieldValue;
        libs ("libfieldFunctionObjects.so");
        patch outlet;
        operation sum;
        fields (phi);
        writeControl writeTime;
        writeFields false;
        log true;
    }
}

Gli script sono distribuiti con licenza GNU GPL versione 3 o successiva; i file del tutorial OpenFOAM non sono inclusi nel pacchetto e restano sotto la loro licenza. Il testo completo dei quattro script è riportato più avanti, nelle sezioni in cui vengono usati.

file scripts/*.sh config/velocities.csv      # segnala "with CRLF line terminators" e "with BOM"
grep -c $'\r' scripts/run_campaign.sh          # righe con CR: deve essere 0
sed -i 's/\r$//' scripts/*.sh                 # correzione se dos2unix non è installato

Il CSV è più tollerante degli script: prepare_campaign.sh elimina da solo un eventuale \r finale da ogni campo, quindi un CSV con fine riga Windows viene letto correttamente. Non tollera invece il BOM che alcuni programmi (Excel in primis) aggiungono all’inizio dei file UTF-8: l’intestazione diventa «\ufeffcase_id», non viene riconosciuta come tale, non supera il controllo sui nomi e lo script si ferma con ERROR: invalid row e codice 3. In Notepad++ scegliere Codifica → UTF-8 (senza BOM) prima di salvare.

bash -n controlla soltanto la sintassi: non esegue nulla, quindi non scopre variabili d’ambiente mancanti o percorsi errati. Quei controlli li fanno gli script stessi all’avvio.

#!/usr/bin/env bash
# GPL-3.0-or-later. Prepare fresh, isolated cases from the official
# OpenFOAM Foundation 14 incompressibleFluid/cylinder tutorial.
set -Eeuo pipefail
root=$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")/.." && pwd -P)
source_case="${FOAM_TUTORIALS:-}/incompressibleFluid/cylinder"
case_root="$root/work/cases"

if ! command -v foamVersion >/dev/null 2>&1 || ! command -v foamDictionary >/dev/null 2>&1; then
  echo "ERROR: source the OpenFOAM Foundation 14 environment first." >&2; exit 2
fi
if [[ "$(foamVersion)" != *14* ]]; then
  echo "ERROR: this package targets OpenFOAM Foundation 14, got: $(foamVersion)" >&2; exit 2
fi
if [[ ! -f "$source_case/0/U" || ! -f "$source_case/system/blockMeshDict" ]]; then
  echo "ERROR: missing official tutorial at $source_case" >&2; exit 2
fi
if [[ ! -d "$case_root" ]]; then mkdir -p -- "$case_root"; fi

while IFS=, read -r name speed || [[ -n "${name:-}" ]]; do
  name=${name%$'\r'}; speed=${speed%$'\r'}
  [[ "$name" == case_id || -z "$name" ]] && continue
  if [[ ! "$name" =~ ^[A-Za-z0-9_-]+$ || ! "$speed" =~ ^(0|[0-9]+)(\.[0-9]+)?([eE][-+]?[0-9]+)?$ ]] ||
     ! awk -v v="$speed" 'BEGIN{exit !(v>0 && v<=1)}'; then
    echo "ERROR: invalid row: $name,$speed" >&2; exit 3
  fi
  case_dir="$case_root/$name"
  if [[ -e "$case_dir" ]]; then
    echo "SKIP existing case $case_dir (no overwrite)"; continue
  fi
  mkdir -- "$case_dir"
  # Copy only case INPUTS, never any pre-existing results or solver logs.
  cp -a -- "$source_case/0" "$source_case/constant" "$source_case/system" "$case_dir/"
  if [[ -d "$case_dir/constant/polyMesh" ]]; then
    echo "ERROR: source unexpectedly contains a mesh; review case $case_dir" >&2; exit 4
  fi
  foamDictionary "$case_dir/0/U" -entry Uinlet -set "$speed" >/dev/null
  cat "$root/config/measurement_functions.foam" >> "$case_dir/system/controlDict"
  printf 'case_id,U_m_s\n%s,%s\n' "$name" "$speed" > "$case_dir/parameters.csv"
  (cd "$case_dir" && sha256sum 0/U system/controlDict system/fvSolution system/fvSchemes constant/physicalProperties constant/momentumTransport > input.sha256)
  printf 'PREPARED\n' > "$case_dir/status.txt"
  echo "CREATED $name U=$speed m/s"
done < "$root/config/velocities.csv"

• set -Eeuo pipefail: il primo comando che fallisce, una variabile non definita o un errore in una pipe interrompono lo script.
• Tre controlli preliminari, tutti con codice di uscita 2: ambiente OpenFOAM caricato (foamVersion e foamDictionary disponibili), versione che contiene «14», tutorial presente con 0/U e system/blockMeshDict.
• Ogni riga del CSV è validata: nome fatto solo di lettere, cifre, trattino e trattino basso; velocità numerica e compresa fra 0 (escluso) e 1 m/s. Una riga non valida ferma tutto con codice 3. Per velocità superiori a 1 m/s va modificata la condizione v>0 && v<=1 nella riga con awk.
• Se la directory del caso esiste, viene stampato SKIP e si passa alla riga successiva: è il meccanismo che impedisce le sovrascritture.
• Si copiano solo 0, constant e system; se nella copia comparisse già una constant/polyMesh lo script si ferma con codice 4, perché la mesh deve essere generata nel caso, non ereditata.
• foamDictionary … -set scrive la velocità, il blocco delle misure viene accodato al controlDict, poi si scrivono parameters.csv, l’impronta input.sha256 e lo stato PREPARED.

L’impronta copre sei file: 0/U, system/controlDict, system/fvSolution, system/fvSchemes, constant/physicalProperties e constant/momentumTransport. Le modifiche a questi file dopo la preparazione bloccano il caso; quelle ad altri file (per esempio system/blockMeshDict o system/decomposeParDict) non vengono rilevate e vanno quindi documentate a mano.

# verifica manuale di un caso preparato
cd work/cases/U015
cat parameters.csv status.txt
sha256sum -c input.sha256
foamDictionary 0/U -entry Uinlet -value

#!/usr/bin/env bash
# GPL-3.0-or-later. Sequential campaign, optional MPI per case.
# Never deletes, force-recreates or overwrites existing case results.
set -Eeuo pipefail
root=$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")/.." && pwd -P)
case_root="$root/work/cases"
np=1
resume=0
while (($#)); do
  case "$1" in
    --np) (($#>=2)) || { echo '--np requires an argument' >&2; exit 2; }; np=$2; shift 2 ;;
    --resume) resume=1; shift ;;
    *) echo "Usage: bash scripts/run_campaign.sh [--np N] [--resume]" >&2; exit 2 ;;
  esac
done
[[ "$np" =~ ^[1-9][0-9]*$ ]] || { echo 'ERROR: --np must be positive integer' >&2; exit 2; }
command -v foamRun >/dev/null 2>&1 || { echo 'ERROR: source OpenFOAM first' >&2; exit 2; }
[[ -d "$case_root" ]] || { echo 'ERROR: run prepare_campaign.sh first' >&2; exit 2; }

execute_case() (
  set -Eeuo pipefail
  case_dir="$1"
  cd "$case_dir"
  name=$(basename "$case_dir")
  if [[ -f .completed ]]; then echo "$name: SKIP completed"; return 0; fi
  if [[ -f .run_started && "$resume" -eq 0 ]]; then
    echo "$name: SKIP unfinished previous run; inspect logs and use --resume explicitly" >&2
    return 8
  fi
  if [[ -f .nproc && "$(cat .nproc)" != "$np" ]]; then
    echo "$name: ERROR, processor count changed; use the original --np" >&2
    return 5
  fi
  if [[ ! -f .nproc ]]; then printf '%s\n' "$np" > .nproc; fi
  stamp=$(date -u +%Y%m%dT%H%M%SZ)_$$
  mkdir -p logs
  printf 'RUNNING\n' > status.txt
  : > .run_started
  fail() { printf 'FAILED\n' > status.txt; echo "$name: FAILED; see logs/$stamp.*" >&2; }
  # Catch accidental edits after preparation: use a fresh identifier for revisions.
  if ! sha256sum -c input.sha256 > "logs/$stamp.input-check.log" 2>&1; then
    echo "$name: ERROR: inputs differ from recorded hashes" >&2
    printf 'FAILED\n' > status.txt
    return 7
  fi
  step() {
    tag=$1; shift
    echo "$name: $tag"
    # Explicit check: bash errexit is disabled inside functions called by `if`.
    if "$@" > "logs/$stamp.$tag.log" 2>&1; then
      return 0
    else
      rc=$?; fail; return "$rc"
    fi
  }
  # Stage stamps prevent repeated mirrorMesh or skipped partial preparations.
  # Before retrying a failed *blockMesh* inspect mesh state: never force-clean.
  if [[ ! -f .block_mesh_done ]]; then
    if [[ -f constant/polyMesh/points ]]; then
      echo 'ERROR: mesh exists without stage marker; inspect manually' >&2; return 4
    fi
    step blockMesh blockMesh || return $?
    : > .block_mesh_done
  fi
  if [[ ! -f .mirror_mesh_done ]]; then
    step mirrorMesh mirrorMesh || return $?
    : > .mirror_mesh_done
  fi
  if [[ ! -f .mesh_checked ]]; then
    step checkMesh checkMesh || return $?
    : > .mesh_checked
  fi
  if ((np>1)); then
    command -v mpirun >/dev/null 2>&1 || { echo 'ERROR: mpirun missing' >&2; return 3; }
    if [[ ! -d processor0 ]]; then step decomposePar decomposePar || return $?; fi
    step foamRun mpirun -np "$np" foamRun -parallel || return $?
    step reconstructPar reconstructPar -latestTime || return $?
  else
    [[ ! -d processor0 ]] || { echo 'ERROR: parallel case cannot be resumed in serial' >&2; return 4; }
    step foamRun foamRun || return $?
  fi
  # Completion here means exit status 0 + data presence, NOT physical convergence.
  for measurement in pIn pOut qIn qOut; do
    find "postProcessing/$measurement" -type f -name '*.dat' -print -quit 2>/dev/null | grep -q . || {
      echo "ERROR: missing measurement $measurement" >&2; fail; return 6;
    }
  done
  printf 'SOLVER_FINISHED\n' > status.txt
  : > .completed
  echo "$name: SOLVER_FINISHED (check physical quality independently)"
)
failures=0
for case_dir in "$case_root"/*; do
  [[ -d "$case_dir" ]] || continue
  if ! execute_case "$case_dir"; then failures=$((failures + 1)); fi
done
if ((failures)); then echo "Campaign has $failures failed case(s)" >&2; exit 1; fi

La funzione execute_case ha il corpo fra parentesi tonde: viene eseguita in una subshell, quindi il cd nel caso e le impostazioni di errore non si propagano al ciclo principale. Un caso che fallisce restituisce un codice diverso da zero, il ciclo conta l’errore e passa al caso successivo; alla fine lo script esce con 1 se almeno un caso è fallito.

La funzione step controlla esplicitamente il codice di ogni comando perché, come dice il commento nello script, in Bash errexit non agisce dentro funzioni chiamate in un if. Ogni fase scrive in logs/<data>_<pid>.<fase>.log, dove la data è in UTC.

Se la connessione cade senza tmux, status.txt resta su RUNNING e il marcatore .run_started è presente: al lancio successivo il caso viene saltato con codice 8 finché non si sceglie esplicitamente --resume.

decomposePar legge system/decomposeParDict e crea tanti sottodomini quanti indicati da numberOfSubdomains. Lo script passa --np solo a mpirun: il dizionario non viene né creato né modificato. Prima di una campagna con --np 4 verificare quindi in ogni caso che il file esista e che il numero coincida, altrimenti decomposePar fallisce oppure crea un numero di processor* diverso dai processi MPI e foamRun -parallel si ferma.

for d in work/cases/*/; do
  printf '%s: ' "$d"
  foamDictionary "$d/system/decomposeParDict" -entry numberOfSubdomains -value 2>/dev/null || echo 'MANCANTE'
done

Un dizionario minimo, da copiare in system/ di ogni caso prima del lancio (il file non è coperto da input.sha256, quindi l’aggiunta va annotata):

FoamFile
{
    format      ascii;
    class       dictionary;
    object      decomposeParDict;
}

numberOfSubdomains 4;

method          scotch;

Con T1 il tempo di calcolo in seriale e Tn quello con n processi, lo speed-up è S = T1/Tn e l’efficienza E = S/n. Un’efficienza vicina a 1 indica un buon uso dei core; valori bassi indicano che comunicazione e sincronizzazione pesano più del calcolo. Come ordine di grandezza pratico, il parallelo comincia a rendere quando ogni processo gestisce almeno decine di migliaia di celle: il numero di celle si legge nel log di checkMesh.

# celle della mesh e tempi registrati da foamRun
grep -m1 "cells:" work/cases/U015/logs/*.checkMesh.log
grep ExecutionTime work/cases/U015/logs/*.foamRun.log | tail -n 1
# carico del server prima di scegliere --np
nproc; free -h; uptime

Non lanciare più processi dei core fisici liberi: con il server sovraccarico i processi MPI si attendono a vicenda e il calcolo diventa più lento del seriale.

Ogni function object scrive in postProcessing/<nome>/<tempo di inizio>/ un file di testo (tipicamente surfaceFieldValue.dat): righe di intestazione che iniziano con #, poi una riga per ogni scrittura con il tempo nella prima colonna e il valore nella seconda. Se il calcolo riparte da un tempo successivo, nasce una nuova sottodirectory con il nuovo tempo di inizio.

collect_results.py legge tutti i *.dat di ciascuna misura, scarta intestazioni e righe non numeriche e prende il campione con il tempo più grande: per questo gestisce anche le ripartenze. Controlla poi che i quattro tempi coincidano e, se non è così, marca il caso INCOMPLETE_DATA.

Con writeControl writeTime si ottiene un campione solo agli istanti di scrittura dei campi. Per una storia più fitta, utile a vedere se Δp si è stabilizzata, si può usare writeControl timeStep; con writeInterval 1; (o un altro intervallo). La modifica va fatta in config/measurement_functions.foam prima della preparazione: dopo, cambierebbe il controlDict registrato nell’impronta e il caso verrebbe bloccato.

ls work/cases/U015/postProcessing/*/*/
tail -n 3 work/cases/U015/postProcessing/pIn/*/*.dat

Con normali uscenti, la portata entrante ha segno negativo. Con valori di esempio Qin = −1,000·10⁻³ m³/s e Qout = +0,998·10⁻³ m³/s la somma vale −2·10⁻⁶ m³/s e lo squilibrio è εQ = 100 · 2·10⁻⁶ / 1,000·10⁻³ = 0,2 %.

In un incomprimibile a densità costante la continuità vale a ogni passo, anche nei transitori: lo squilibrio riflette quindi soprattutto la tolleranza del solutore della pressione. Un valore grande di solito non è fisico: segnala una patch aperta non inclusa nella somma, un segno sbagliato o due campioni presi a tempi diversi. Se il dominio ha altre patch con flusso, vanno aggiunte al bilancio.

#!/usr/bin/env python3
# GPL-3.0-or-later. No third-party dependencies. Run on Windows or Linux.
"""Read the last numerical records of four OpenFOAM surfaceFieldValue outputs.

Values are NOT convergence/validation certificates. CFD runs are required
before real measurements exist. Assumes OpenFOAM Foundation 14 cylinder case:
p is kinematic pressure [m2/s2]; phi is volumetric flux [m3/s].
"""
import argparse
import csv
import math
from pathlib import Path


def read_measurement(case: Path, name: str):
    files = sorted((case / 'postProcessing' / name).glob('**/*.dat'))
    samples = []
    for file in files:
        for raw in file.read_text(encoding='utf-8', errors='replace').splitlines():
            line = raw.strip()
            if not line or line.startswith('#'):
                continue
            fields = line.split()
            if len(fields) < 2:
                continue
            try:
                t, value = float(fields[0]), float(fields[1])
            except ValueError:
                continue
            if math.isfinite(t) and math.isfinite(value):
                samples.append((t, value))
    if not samples:
        raise ValueError(f'No numerical rows found for {name}')
    return max(samples, key=lambda r: r[0])


def main():
    p = argparse.ArgumentParser()
    p.add_argument('--root', type=Path, default=Path(__file__).resolve().parents[1])
    p.add_argument('--rho', type=float, default=1000.0,
                   help='Fluid density kg/m3, default 1000; NOT inferred from tutorial')
    p.add_argument('--output', type=Path, help='CSV destination; default ROOT/results.csv')
    args = p.parse_args()
    if args.rho <= 0 or not math.isfinite(args.rho):
        p.error('rho must be a positive finite number')
    case_root = args.root / 'work' / 'cases'
    if not case_root.is_dir():
        p.error(f'No cases found: {case_root}')
    fields = ['case_id', 'U_m_s', 'status', 'time', 'pIn_m2_s2', 'pOut_m2_s2',
              'DeltaP_Pa', 'qIn_m3_s', 'qOut_m3_s', 'imbalance_percent', 'notes']
    records = []
    for case in sorted(p for p in case_root.iterdir() if p.is_dir()):
        parameters = case / 'parameters.csv'
        if not parameters.is_file():
            continue
        with parameters.open(encoding='utf-8', newline='') as f:
            cfg = next(csv.DictReader(f))
        row = {key: '' for key in fields}
        row.update(case_id=case.name, U_m_s=cfg['U_m_s'])
        status_file = case / 'status.txt'
        row['status'] = status_file.read_text(encoding='utf-8').strip() if status_file.exists() else 'UNKNOWN'
        if not (case / '.completed').exists():
            row['notes'] = 'No .completed marker; do not treat as finished'
        else:
            try:
                out = {n: read_measurement(case, n) for n in ('pIn', 'pOut', 'qIn', 'qOut')}
                times = [out[n][0] for n in out]
                if max(times) - min(times) > 1e-7 * max(1.0, max(abs(t) for t in times)):
                    raise ValueError('Measurement times differ; inspect files')
                pin, pout, qin, qout = (out[n][1] for n in ('pIn', 'pOut', 'qIn', 'qOut'))
                row['time'] = format(max(times), '.10g')
                row['pIn_m2_s2'] = format(pin, '.10g')
                row['pOut_m2_s2'] = format(pout, '.10g')
                row['DeltaP_Pa'] = format((pin-pout)*args.rho, '.10g')
                row['qIn_m3_s'] = format(qin, '.10g')
                row['qOut_m3_s'] = format(qout, '.10g')
                denom = max(abs(qin), abs(qout))
                row['imbalance_percent'] = format(100*abs(qin+qout)/denom, '.6g') if denom else ''
                row['notes'] = 'Last recorded sample, NOT a convergence test'
            except (ValueError, OSError) as exc:
                row['status'] = 'INCOMPLETE_DATA'
                row['notes'] = str(exc)
        records.append(row)
    output = args.output or args.root / 'results.csv'
    output.parent.mkdir(parents=True, exist_ok=True)
    with output.open('w', encoding='utf-8', newline='') as f:
        writer = csv.DictWriter(f, fieldnames=fields)
        writer.writeheader()
        writer.writerows(records)
    print(f'Wrote {len(records)} row(s): {output}')
    for row in records:
        print(row['case_id'], row['status'], row['DeltaP_Pa'], row['notes'])


if __name__ == '__main__':
    main()

• --root: cartella del pacchetto (di default quella che contiene scripts/); --rho: densità in kg/m³, deve essere positiva e finita; --output: destinazione del CSV, di default results.csv nella cartella del pacchetto.
• Un caso senza .completed compare nel CSV con lo stato letto da status.txt e la nota «No .completed marker», senza numeri.
• Se i tempi delle quattro misure differiscono oltre una tolleranza relativa di 10⁻⁷, il caso diventa INCOMPLETE_DATA.
• I numeri sono scritti con 10 cifre significative e il punto decimale.

Il file usa la virgola come separatore e il punto come separatore decimale. Un Excel configurato in italiano, aperto con doppio clic, si aspetta il punto e virgola e la virgola decimale e mette tutto in una colonna o trasforma i numeri in date. Usare invece Dati → Da testo/CSV, scegliere il delimitatore «Virgola» e, in Trasforma dati, impostare per le colonne numeriche la lingua «Inglese (Stati Uniti)». In alternativa, LibreOffice Calc chiede separatore e lingua all’apertura.

#!/usr/bin/env python3
# GPL-3.0-or-later. Optional matplotlib package on Windows/Linux.
import argparse
import csv
import math
from pathlib import Path

p = argparse.ArgumentParser()
p.add_argument('csvfile', type=Path, nargs='?', default=Path('results.csv'))
p.add_argument('--output', type=Path, default=Path('DeltaP_vs_U.png'))
a = p.parse_args()
with a.csvfile.open(encoding='utf-8-sig', newline='') as f:
    rows = list(csv.DictReader(f))
points = []
for r in rows:
    if r['status'] != 'SOLVER_FINISHED' or not r['DeltaP_Pa']:
        continue
    try:
        u, dp = float(r['U_m_s']), float(r['DeltaP_Pa'])
    except ValueError:
        continue
    if math.isfinite(u) and math.isfinite(dp):
        points.append((u, dp))
if not points:
    raise SystemExit('No finished CFD results; refusing to create a misleading plot')
points.sort()
try:
    import matplotlib.pyplot as plt
except ImportError as e:
    raise SystemExit('Install matplotlib: python -m pip install matplotlib') from e
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot([p[0] for p in points], [p[1] for p in points], '-o')
ax.set_xlabel('Inlet velocity [m/s]')
ax.set_ylabel('Pressure difference [Pa]')
ax.set_title('OpenFOAM campaign - final samples (not validated)')
ax.grid(True, alpha=.3)
fig.tight_layout()
fig.savefig(a.output, dpi=175)
print('Saved:', a.output)

Lo script si rifiuta di produrre il grafico se nessun caso è SOLVER_FINISHED con un valore di Δp: un grafico vuoto o fatto di casi incompleti sarebbe fuorviante. Il nome del file di uscita si cambia con --output.

Quando prevalgono gli effetti viscosi la caduta di pressione cresce circa linearmente con la velocità; quando prevalgono quelli inerziali cresce circa con il quadrato. L’esponente locale fra due casi consecutivi è n = ln(ΔP₂/ΔP₁) / ln(U₂/U₁): valori vicini a 1 indicano regime viscoso, vicini a 2 regime inerziale, valori intermedi una transizione. È un controllo di plausibilità fisica che non richiede dati esterni.

# esponenti locali da results.csv (solo libreria standard)
import csv, math
rows = [r for r in csv.DictReader(open("results.csv", encoding="utf-8"))
        if r["status"] == "SOLVER_FINISHED" and r["DeltaP_Pa"]]
pts = sorted((float(r["U_m_s"]), float(r["DeltaP_Pa"])) for r in rows)
for (u1, p1), (u2, p2) in zip(pts, pts[1:]):
    if p1 > 0 and p2 > 0:
        print(f"{u1:g} -> {u2:g} m/s: n = {math.log(p2/p1)/math.log(u2/u1):.3f}")

Con --resume lo script accetta i casi che hanno .run_started ma non .completed. Le fasi già concluse (.block_mesh_done, .mirror_mesh_done, .mesh_checked) non vengono ripetute; in parallelo, se processor0 esiste, decomposePar viene saltato e foamRun -parallel riparte dai sottodomini esistenti. L’impronta degli input viene ricontrollata a ogni lancio.

Conseguenza pratica: se durante un calcolo si modifica il controlDict (per esempio cambiando endTime o stopAt per fermarlo), il file non corrisponde più all’impronta e al --resume successivo il caso si ferma con codice 7. È il comportamento voluto: la modifica va annotata e, di norma, si crea un caso nuovo. Rigenerare input.sha256 è possibile ma deve essere una decisione consapevole e scritta nel registro della campagna.

# processi ancora attivi e loro stato
pgrep -af foamRun
find work/cases -name status.txt -exec sh -c 'echo "$1: $(cat "$1")"' _ {} \;
ls -a work/cases/U015

Eseguiti nella cartella del pacchetto a campagna conclusa, registrano versioni e impronte e producono un archivio senza le directory processor*, che occupano molto spazio e si ricostruiscono.

{ date -u; uname -a; foamVersion; mpirun --version 2>&1 | head -n 1; } > campaign_environment.txt
find work/cases -path '*/processor*' -prune -o -type f -print0 | sort -z | xargs -0 sha256sum > campaign_files.sha256
tar czf "cylinder_campaign_$(date -u +%Y%m%d).tar.gz" --exclude='processor*' config scripts work results.csv campaign_environment.txt campaign_files.sha256
sha256sum cylinder_campaign_*.tar.gz

L’impronta dell’archivio, annotata nel registro, permette di dimostrare in seguito che i dati analizzati sono quelli prodotti. Se results.csv non esiste ancora, toglierlo dalla riga di tar.

• A. Al secondo lancio di prepare_campaign.sh le quattro righe originali producono SKIP existing case … e solo la nuova produce CREATED U025 U=0.025 m/s. Verifica: foamDictionary work/cases/U025/0/U -entry Uinlet -value deve restituire 0.025.
• B. I casi vivono in work/cases dentro la cartella del pacchetto, quindi due campagne indipendenti si ottengono con due copie del pacchetto (per esempio cylinder_campaign_np1 e cylinder_campaign_np4). Il marcatore .nproc impedisce comunque di mescolare le due modalità nello stesso caso. Per la versione con 4 processi preparare prima il decomposeParDict.
• C. Oltre a imbalance_percent, calcolare l’esponente locale fra casi consecutivi e controllare che il tempo della misura sia lo stesso in tutti i casi: confrontare un caso a regime con uno ancora in transitorio produce una curva senza significato.