OpenFOAM 14, parte 6: fvSolution, solutori lineari e algoritmi di soluzione

Fino alla versione 10 della distribuzione Foundation ogni combinazione di fisica e algoritmo aveva un eseguibile dedicato: simpleFoam per i flussi incomprimibili stazionari con SIMPLE, pisoFoam e pimpleFoam per quelli transitori. Dalla versione 11 questi eseguibili sono stati sostituiti da moduli solutore caricati da un unico programma, foamRun; il modulo si sceglie in system/controlDict con la voce solver incompressibleFluid; (oppure con foamRun -solver incompressibleFluid).

La conseguenza per fvSolution è importante: in v14 non è più il nome dell’eseguibile a dire se l’algoritmo è SIMPLE o PIMPLE. Lo stesso modulo incompressibleFluid gestisce entrambi i casi: il comportamento stazionario si ottiene dichiarando ddtSchemes { default steadyState; } in fvSchemes e scrivendo un sotto-dizionario SIMPLE in fvSolution, come nel tutorial pitzDailySteady della guida v14. Le voci dei solutori lineari (solvers), invece, sono rimaste sostanzialmente le stesse: per questo un vecchio fvSolution è spesso ancora leggibile, ma il suo contesto (eseguibile, file dei modelli, nomi delle dimensioni) va verificato.

foamToC -solvers          # list the solver modules available to foamRun
foamInfo incompressibleFluid
foamDictionary -entry solver -value system/controlDict

Secondo la User Guide v14 (§4.6), fvSolution contiene un sotto-dizionario solvers con le impostazioni del solutore lineare di ciascuna equazione, un sotto-dizionario relaxationFactors e un sotto-dizionario di controllo dell’algoritmo, SIMPLE per i casi stazionari oppure PIMPLE per quelli transitori. Altri moduli possono leggere ulteriori voci (per esempio i moduli multifase leggono controlli specifici per la frazione di fase), ma per incompressibleFluid questi tre blocchi sono tutto ciò che serve.

Il file non viene letto una sola volta. Le impostazioni di un solutore lineare vengono cercate nel dizionario solvers ogni volta che un’equazione viene risolta, e con runTimeModifiable attivo (il comportamento abituale, vedi Parte 4) OpenFOAM rilegge i file di system/ quando vengono modificati durante l’esecuzione. Si può quindi, per esempio, abbassare un fattore di rilassamento mentre il calcolo gira; la modifica ha effetto all’iterazione successiva e il log lo segnala con una riga che indica la rilettura del dizionario. È comodo, ma rende meno riproducibile la cronologia della simulazione: annota sempre la modifica e il tempo in cui l’hai fatta.

Errore tipico: se manca il blocco di un campo che il modulo deve risolvere, l’esecuzione si ferma con un FOAM FATAL IO ERROR del tipo «keyword p is undefined in dictionary .../system/fvSolution/solvers». Il messaggio indica sia il nome cercato sia il percorso del dizionario: è la prima cosa da leggere.

Consideriamo la diffusione stazionaria monodimensionale d/dx(Γ dφ/dx) = 0 su un segmento lungo 3 diviso in tre celle uguali (Δx = 1, Γ = 1), con φ = 1 al contorno sinistro e φ = 0 al contorno destro. Integrando sulla cella P e approssimando il flusso diffusivo su ciascuna faccia con la differenza tra i centri, Γ(φ_N − φ_P)/d, si ottiene un’equazione per cella. Per una faccia interna d = Δx = 1 e il coefficiente vale 1; per una faccia di contorno il centro della cella dista d = Δx/2 dalla faccia e il coefficiente vale 2.

cell 1:  3*phi1 - 1*phi2           = 2   (2 = 2*1, from the left boundary value)
cell 2: -1*phi1 + 2*phi2 - 1*phi3  = 0
cell 3:          -1*phi2 + 3*phi3  = 0   (right boundary value is 0)

    | 3 -1  0 |        | 2 |
A = |-1  2 -1 |    b = | 0 |      x = (5/6, 1/2, 1/6)
    | 0 -1  3 |        | 0 |

La soluzione, x = (0,833; 0,5; 0,167), coincide con il profilo lineare esatto 1 − x/3 valutato nei centri delle celle (0,5; 1,5; 2,5). Nella notazione usata qui, a_P è il termine diagonale (3, 2, 3), gli a_N sono i coefficienti fuori diagonale (−1) e b_P raccoglie i contributi noti, qui i valori al contorno. Il valore assoluto della diagonale è maggiore o uguale alla somma dei valori assoluti fuori diagonale in ogni riga: la matrice è diagonalmente dominante, ed è simmetrica perché il coefficiente che lega P a N è lo stesso che lega N a P. Con un milione di celle esaedriche la struttura è identica: ogni riga ha un elemento diagonale più uno per ciascuna faccia interna della cella (al massimo sei), quindi circa 7 milioni di elementi non nulli contro i 10¹² di una matrice densa.

OpenFOAM non usa un formato sparso generico (come CSR), ma il formato LDU della classe lduMatrix, costruito direttamente sulla topologia della mesh: un vettore diag con un valore per cella, un vettore upper e un vettore lower con un valore per ogni faccia interna. Ogni faccia interna ha una cella owner (proprietaria, indice minore) e una cella neighbour (vicina): il coefficiente upper della faccia sta nella riga dell’owner e nella colonna del neighbour, il coefficiente lower nella posizione trasposta. Se la matrice è simmetrica, lower non viene nemmeno allocato.

Le facce di contorno non generano elementi fuori diagonale: la condizione al contorno (Parte 3) contribuisce con una parte implicita alla diagonale e con una parte esplicita al termine noto b. Per questo un fixedValue rafforza la diagonale della cella adiacente (come il 3 nell’esempio a tre celle), mentre un zeroGradient non aggiunge nulla. Le patch di processore nei calcoli paralleli (vedi il caso pratico sul calcolo parallelo MPI) sono trattate come «interfacce» che scambiano i valori tra sottodomini a ogni prodotto matrice-vettore: è il motivo per cui il numero di iterazioni di un solutore lineare può cambiare leggermente con il numero di processori.

Il «può» della verifica di previsione è preciso. In OpenFOAM molti schemi di ordine superiore sono implementati come correzione differita: la parte implicita, che finisce nella matrice, è quella dello schema upwind, mentre la differenza tra lo schema scelto e l’upwind viene calcolata con i valori correnti ed entra nel termine noto b. È il caso di linearUpwind: passando da bounded Gauss upwind a bounded Gauss linearUpwind grad(U) la matrice A dell’equazione della quantità di moto resta strutturalmente quella dell’upwind e cambia soprattutto b. Passando invece a Gauss linear (differenze centrate) cambiano i pesi di interpolazione e quindi i coefficienti di A: la dominanza diagonale si indebolisce e i solutori iterativi possono richiedere più iterazioni o diventare instabili su mesh grossolane.

Il secondo enunciato vale in entrambe le direzioni: un solutore lineare diverso, se converge alla stessa tolleranza, produce la stessa soluzione discreta a meno dell’errore algebrico residuo. Se cambiando solo il solutore lineare cambia in modo apprezzabile il risultato finale (per esempio la caduta di pressione), le tolleranze erano troppo larghe: è un test utile per verificarlo.

La tabella qui sopra elenca le voci essenziali. La User Guide v14 ne documenta altre, con i relativi valori predefiniti:

Il solutore diagonal serve per le equazioni esplicite, la cui matrice è solo diagonale: non itera e non usa tolleranze. Per vedere quali nomi accetta la tua installazione si può provocare di proposito un errore su una copia (per esempio solver xyz;): OpenFOAM risponde con l’elenco dei solutori validi per quel tipo di matrice. Per controllare che cosa il modulo leggerà davvero dopo l’espansione delle macro, usa foamDictionary -expand -entry solvers system/fvSolution.

Quando il modulo deve risolvere, per esempio, l’equazione di k, cerca la voce k nel dizionario solvers con questa precedenza: prima una chiave con il nome esatto; se non c’è, le chiavi tra virgolette che sono espressioni regolari, esaminate dall’ultima alla prima nell’ordine del file. L’espressione deve corrispondere al nome per intero: "(U|k|epsilon)" corrisponde a k ma non a kFinal, ed è per questo che serve il secondo gruppo "(U|k|epsilon)Final".

La macro $p; copia nel blocco corrente tutte le voci del blocco p; le voci scritte dopo la macro le sovrascrivono (qui relTol 0;). Anche $U; funziona sebbene non esista una chiave letterale U: la ricerca della macro accetta le chiavi regex, quindi trova "(U|k|epsilon)". Il risultato espanso si controlla così:

foamDictionary -expand -entry solvers/pFinal system/fvSolution
foamDictionary -expand -entry solvers system/fvSolution     # whole block, macros expanded

Nel blocco PIMPLE di questo estratto nOuterCorrectors non compare: vale il predefinito 1, cioè la modalità «PISO» (un solo passaggio esterno per passo temporale con due correzioni di pressione), adatta quando il numero di Courant resta sotto 1. In un caso laminare come il canale di riferimento k ed epsilon non esistono: il gruppo regex è innocuo ma inutile, e conviene ridurlo a U e UFinal.

Il residuo «grezzo» di un sistema è il vettore b − A x, che ha le unità dell’equazione e dipende dalla dimensione delle celle: un suo valore assoluto non direbbe nulla. OpenFOAM lo rende adimensionale dividendolo per un fattore di normalizzazione costruito con il valore medio x̄ del campo corrente:

            sum_cells | b - A x |
r = ---------------------------------------------------
     sum_cells ( | A x - A xbar | + | b - A xbar | ) + small

Il denominatore misura quanto il problema «si discosta» da una soluzione uniforme. Due conseguenze pratiche. Primo: se il campo iniziale è uniforme (x = x̄, come internalField uniform 0), il primo termine del denominatore si annulla e il secondo coincide con il numeratore, quindi il residuo iniziale del primo passo vale 1: non è un errore, è la normalizzazione. Secondo: il residuo è una somma (norma L1) su tutte le celle divisa per una somma analoga, quindi è una media; residui piccoli possono nascondere poche celle con errori grandi, motivo in più per i controlli fisici delle sezioni successive.

Su maxIter: se nel log compare No Iterations 1000 (o il valore che hai impostato), il solutore non ha raggiunto né tolerance né relTol ed è stato fermato dal tetto. Non è un messaggio d’errore, ma è quasi sempre il sintomo di un problema: matrice mal condizionata, solutore inadatto, mesh di cattiva qualità o una soluzione che sta divergendo.

Molti metodi iterativi, dopo le prime iterazioni, riducono il residuo di un fattore circa costante ρ a ogni iterazione (il fattore di convergenza). Nella tabella il primo passo porta il residuo da 2,0 × 10⁻³ a 8,0 × 10⁻⁴, cioè ρ ≈ 0,4. Se ρ restasse 0,4, dopo n iterazioni il rapporto sarebbe 0,4ⁿ, e si possono stimare le iterazioni necessarie per ciascun criterio:

La stima mostra la regola generale: ogni decade in più di riduzione costa circa ln 10 / |ln ρ| iterazioni (qui circa 2,5). Il valore reale nella tabella scende più in fretta (1,5 × 10⁻⁵ corrisponde a un rapporto 0,0075), perché ρ non è costante: con GAMG in particolare ogni ciclo riduce il residuo molto più di un fattore 0,4. Il modo corretto di misurare ρ sul proprio caso è leggere nel log residuo iniziale, finale e numero di iterazioni: ρ ≈ (r_final/r_initial)^1/N. Con le righe di log della sezione 7: per p, (0,00008/0,014)^1/5 ≈ 0,36.

Ecco come appare, in forma tipica, un passo di un caso transitorio incompressibleFluid 3D con PIMPLE in modalità PISO (nOuterCorrectors 1, nCorrectors 2). I numeri sono illustrativi; la struttura è quella che troverai nel tuo log.

Courant Number mean: 0.12 max: 0.45
Time = 0.0025s

smoothSolver:  Solving for Ux, Initial residual = 0.0032, Final residual = 0.00021, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0051, Final residual = 0.00030, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.0047, Final residual = 0.00028, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.014, Final residual = 0.00008, No Iterations 5
time step continuity errors : sum local = 2.1e-07, global = -3.4e-09, cumulative = -3.4e-09
GAMG:  Solving for p, Initial residual = 0.0021, Final residual = 9.1e-08, No Iterations 9
time step continuity errors : sum local = 3.0e-10, global = 1.2e-12, cumulative = -3.4e-09
ExecutionTime = 0.41 s  ClockTime = 0 s

• Le tre righe Ux, Uy, Uz sono il predittore della quantità di moto: l’equazione vettoriale viene risolta componente per componente con la stessa matrice. In un caso 2D con patch empty la componente normale al piano non viene risolta e la sua riga manca.
• Le due righe p sono le due correzioni di pressione (nCorrectors 2); la seconda usa le impostazioni di pFinal (sezione 8), quindi va fino a tolerance e richiede più iterazioni.
• time step continuity errors è l’indicatore di bilancio di massa richiesto dal riquadro «Controlli fisici»: sum local è la media pesata del valore assoluto della divergenza dei flussi, global la sua somma con segno, cumulative la somma nel tempo. Valori che crescono passo dopo passo indicano che la pressione non è risolta abbastanza.
• Con un precondizionatore il nome compare come prefisso: DICPCG: Solving for p, DILUPBiCGStab: Solving for Ux.

Per confrontare due esecuzioni (come richiede il Laboratorio A) servono numeri, non impressioni. Alcuni modi, dal più semplice al più strutturato:

# 1. every pressure solve: initial residual and iteration count
grep "Solving for p" log.foamRun | awk -F'[=,]' '{print $3, $NF}' | head

# 2. total linear iterations for p over the whole run
grep "Solving for p" log.foamRun | awk '{s += $NF} END {print s}'

# 3. continuity errors only
grep "continuity errors" log.foamRun | tail -5

# 4. split the log into one file per quantity (logs/p_0, logs/Ux_0, ...)
foamLog log.foamRun

In alternativa la funzione residuals scrive i residui iniziali in un file di testo durante l’esecuzione: aggiungi #includeFunc residuals(p, U) al sotto-dizionario functions di controlDict (Parte 4) e troverai i dati in postProcessing/residuals/; controlla le opzioni con foamInfo residuals. Ricorda che questi file registrano il residuo iniziale di ogni soluzione: è la cronologia «attraverso i passi», non quella «dentro una soluzione».

Nell’implementazione Foundation il suffisso Final viene aggiunto al nome del campo quando l’algoritmo si trova nell’ultima iterazione: per la velocità e gli altri campi durante l’ultimo ciclo esterno PIMPLE (con nOuterCorrectors 1, quindi, a ogni passo); per la pressione nell’ultima correzione di pressione dell’ultimo ciclo esterno (e, se ci sono correttori non ortogonali, nell’ultimo di essi). La User Guide v14 descrive lo scopo: risolvere l’equazione con impostazioni separate per il correttore finale, tipicamente relTol 0, in modo che ogni passo temporale converga fino a tolerance.

Perché conta nei casi transitori: l’errore algebrico lasciato a fine passo non viene più corretto e si accumula nel tempo come errore di continuità. Le correzioni intermedie, invece, verranno comunque superate dalla successiva: è inutile pagarle a caro prezzo. Nel log della sezione 7 la differenza è visibile: 5 iterazioni con relTol 0.01, 9 con relTol 0.

Se in un caso transitorio la voce manca, l’esecuzione si interrompe alla prima soluzione finale con un errore del tipo «keyword pFinal is undefined in dictionary .../system/fvSolution/solvers». Nei casi stazionari con SIMPLE (come il File didattico A della sezione 18 e il tutorial pitzDailySteady) le voci Final normalmente non compaiono.

• Le virgolette doppie sono obbligatorie: è la virgolettatura che dice a OpenFOAM di trattare la chiave come espressione regolare. Senza virgolette, (U|k|epsilon) non è nemmeno una parola valida.
• La corrispondenza è sull’intero nome: "U.*" cattura U, UFinal ma anche, per esempio, un eventuale Ua; ".*" cattura tutto ed è legittimo solo in contesti come relaxationFactors (sezione 19).
• Il punto è un metacarattere: tra virgolette, "alpha.water" corrisponde anche a alphaXwater. Quando conta il punto letterale, la soluzione più sicura è scrivere il nome esatto senza virgolette (alpha.water), che viene confrontato carattere per carattere.
• Precedenza: il nome esatto vince sempre; tra più regex vince l’ultima nel file. Un blocco esplicito U scritto dopo il gruppo prevale quindi sul gruppo per U.

Quali campi vadano coperti dipende dal modello in constant/momentumTransport (Parte 7):

nut, la viscosità turbolenta, è calcolata dal modello e non ha equazione propria: non serve una voce in solvers. Se un campo risolto non è coperto, l’errore arriva alla prima soluzione di quel campo, con il messaggio «keyword omega is undefined in dictionary …/solvers».

Riprendiamo le tre celle della sezione 2 e aggiungiamo un flusso convettivo F = 1 diretto verso +x, discretizzato con l’upwind: sulla faccia di uscita di ogni cella il valore trasportato è quello della cella stessa (contributo +F sulla diagonale), sulla faccia di ingresso è quello della cella a monte (contributo −F sul coefficiente del vicino a monte). La riga della cella 2 diventa (2 + F) φ₂ − (1 + F) φ₁ − φ₃ = 0:

diffusion only            diffusion + upwind convection (F = 1)
| 3 -1  0 |               | 4 -1  0 |
|-1  2 -1 |               |-2  3 -1 |
| 0 -1  3 |               | 0 -2  4 |
symmetric: A12 = A21      asymmetric: A21 = -2, A12 = -1

L’equazione della quantità di moto contiene sempre il termine convettivo div(phi,U), quindi la sua matrice è asimmetrica; l’equazione della pressione degli algoritmi SIMPLE/PIMPLE è invece un’equazione di Laplace (o di Poisson) ∇·((1/A)∇p) = ∇·(HbyA), con una matrice simmetrica. Ecco perché nei tutorial p usa GAMG o PCG, mentre U, k, epsilon usano smoothSolver o PBiCGStab.

Una precisazione sulla tabella: la User Guide v14 abbina PCG alle matrici simmetriche e PBiCG a quelle asimmetriche, ma presenta PBiCGStab come variante stabilizzata utilizzabile per entrambe le classi; l’output di foamSearch riportato nella guida mostra infatti tutorial che lo usano anche per p. La riga «asimmetrica» della tabella va letta come «adatta anche alle asimmetriche». GAMG e smoothSolver accettano entrambe le classi, a patto di scegliere uno smoother compatibile (sezione 11).

Secondo la guida v14, usare un solutore simmetrico su una matrice asimmetrica (o viceversa) produce un FOAM FATAL IO ERROR. La forma tipica è la seguente; l’elenco esatto dipende dalla versione e dalle librerie caricate:

--> FOAM FATAL IO ERROR:
Unknown asymmetric matrix solver PCG

Valid asymmetric matrix solvers are :
(
GAMG
PBiCG
PBiCGStab
smoothSolver
)

file: .../system/fvSolution/solvers/U from line ... to line ...

Leggilo come insegna il riquadro: la parola chiave sbagliata è solver, il file e il blocco sono indicati nell’ultima riga, la classe di matrice («asymmetric») è scritta nel messaggio e le alternative valide sono elencate. La correzione è una sola: cambiare solver (e di conseguenza preconditioner, perché DIC vale solo per matrici simmetriche: per PBiCGStab si usa DILU).

Gauss-Seidel. Si percorrono le celle in ordine e si aggiorna ciascuna risolvendo la propria equazione con i valori più recenti dei vicini: φ_P ← (b_P − Σ a_N φ_N) / a_P. Ogni aggiornamento propaga l’informazione di una sola cella, quindi gli errori che oscillano da cella a cella (alta frequenza) spariscono in poche passate, mentre quelli estesi su tutto il dominio decadono lentissimamente. symGaussSeidel esegue una passata in avanti e una all’indietro, eliminando la direzionalità legata alla numerazione delle celle: è lo smoother preferito nei tutorial per U e i campi di turbolenza. La voce nSweeps (predefinito 1) fissa quante passate fare prima di ricalcolare il residuo.

DIC e DILU. Sono fattorizzazioni incomplete: si costruisce un’approssimazione M ≈ A = L D U che ha la stessa struttura sparsa di A (nessun riempimento) e in cui, nella variante «diagonale» usata da OpenFOAM, viene modificata solo la diagonale. Risolvere con M costa quanto due passate di sostituzione, e il metodo di Krylov lavora sul sistema precondizionato M⁻¹A x = M⁻¹b, molto meglio condizionato. DIC (Cholesky incompleta) richiede una matrice simmetrica; DILU (LU incompleta) è la versione per matrici asimmetriche; FDIC è una versione più veloce di DIC che precalcola le quantità reciproche. Gli stessi oggetti possono fungere da smoother.

DICGaussSeidel applica uno smoothing DIC seguito da Gauss-Seidel ed è la scelta classica come smoother di GAMG per la pressione: essendo basato su DIC, è per matrici simmetriche (per un’equazione asimmetrica risolta con GAMG si usa GaussSeidel o symGaussSeidel). La guida v14 elenca inoltre GAMG stesso fra i precondizionatori (un ciclo multigrid per precondizionare PCG), diagonal (in genere non usato) e none.

Scomponiamo l’errore e = x − x_esatto in modi di Fourier. Un’iterazione di Gauss-Seidel smorza di molto i modi la cui lunghezza d’onda è di poche celle, ma un modo esteso su n celle si riduce solo di un fattore circa 1 − c/n² per passata: su una mesh di 200 celle per lato servirebbero decine di migliaia di passate. L’idea del multigrid è che un errore liscio sulla griglia fine diventa oscillante su una griglia più grossolana: dopo poche passate di smoothing (pre-sweeps) il residuo viene trasferito a un livello con meno celle, dove lo stesso smoother è di nuovo efficace; scendendo fino al livello più grossolano (circa nCellsInCoarsestLevel celle) il sistema è così piccolo da essere risolto quasi esattamente; la correzione viene poi interpolata verso l’alto e rifinita con altre passate (post-sweeps). Il risultato è un numero di cicli quasi indipendente dalla dimensione della mesh e un costo proporzionale al numero di celle.

Confronto con PCG: il numero di iterazioni del gradiente coniugato cresce come la radice del numero di condizionamento, e per un laplaciano questo cresce come (L/h)²; dimezzare la dimensione delle celle raddoppia quindi circa le iterazioni di PCG, mentre quelle di GAMG restano quasi uguali. È per questo che GAMG domina per la pressione su mesh medio-grandi, mentre su mesh piccole (come il canale didattico da poche migliaia di celle) la differenza è minima. «Algebrico» significa che i livelli grossolani non sono mesh costruite geometricamente ma agglomerazioni di celle decise dai coefficienti della matrice (nei tutorial con l’agglomeratore faceAreaPair, che unisce le celle a coppie attraverso le facce più grandi).

Nota per la lettura della guida: la User Guide v14 scrive questa voce come directSolveCoarset; il nome usato nel codice e nei dizionari è directSolveCoarsest, ma prima di usarla verificane la grafia nei tutorial con foamSearch. Per il primo esperimento il consiglio del riquadro resta valido: questi valori predefiniti vanno bene nella stragrande maggioranza dei casi.

Il rilassamento implicito delle equazioni (alla Patankar) divide la diagonale per α e aggiunge al termine noto un contributo che contiene il valore precedente:

(aP/alpha) * phiP + sum(aN * phiN) = bP + ((1 - alpha)/alpha) * aP * phiP_old

A convergenza φ_P = φ_P,old e i due termini aggiunti si cancellano: la soluzione finale è la stessa, cambia solo il percorso. Con α = 0,7 la diagonale cresce di 1/0,7 ≈ 1,43 volte: la matrice diventa più diagonalmente dominante (i solutori iterativi convergono meglio) e ogni iterazione fa un passo più corto. In OpenFOAM, prima di applicare α, la diagonale viene anche portata almeno al valore della somma dei moduli dei coefficienti fuori diagonale: per questo la guida v14 indica per i casi transitori un fattore 1 per tutte le equazioni (equations { ".*" 1; }, come nel File didattico B), che non rallenta nulla ma assicura quella che la guida chiama «uguaglianza diagonale», in generale necessaria alla convergenza di un’equazione matriciale.

Esempio sul campo. Se in una cella p_vecchio = 10 m²/s² e la soluzione dell’equazione della pressione propone p_candidato = 20 m²/s², con α = 0,3 il valore accettato è 10 + 0,3 × 10 = 13 m²/s². Indicazioni della guida v14: nei tutorial stazionari si usano 0,3 per p e 0,7 per U e per i campi di turbolenza (raggruppati con una regex); valori fino a 0,9 possono essere stabili in alcuni casi, mentre valori molto sotto 0,2 rendono la convergenza proibitivamente lenta. Un campo o un’equazione assenti da relaxationFactors non vengono rilassati.

Scriviamo l’equazione discreta della quantità di moto di ogni cella separando la diagonale A dal resto: A U_P = H(U) − ∇p, dove H raccoglie i contributi dei vicini e le sorgenti. Da qui U = HbyA − (1/A)∇p, con HbyA = H/A. Imponendo la continuità ∇·U = 0 si ottiene l’equazione della pressione. Un’iterazione SIMPLE di incompressibleFluid procede così:

• Predittore della quantità di moto: si assembla l’equazione di U con la pressione corrente, la si rilassa (fattore equations/U) e la si risolve con le impostazioni solvers/U (righe Ux, Uy, Uz nel log).
• Equazione della pressione: ∇·((1/A)∇p) = ∇·(HbyA), risolta con solvers/p, ripetuta nNonOrthogonalCorrectors volte in più se richiesto.
• Correzione dei flussi: il flusso alle facce phi viene corretto con il gradiente di pressione appena calcolato, e diventa conservativo a meno della tolleranza della pressione.
• Rilassamento della pressione (fattore fields/p, per esempio 0,3) e correzione della velocità U = HbyA − (1/A)∇p.
• Equazioni dei modelli (turbolenza ecc.), poi l’iterazione successiva.

Nella formulazione SIMPLEC (consistent yes;) al posto di 1/A si usa 1/(A − ΣH₁), dove ΣH₁ è la somma dei coefficienti dei vicini: si tiene conto in modo approssimato della risposta delle celle vicine, la correzione di pressione è più «coerente» e il rilassamento della pressione diventa in gran parte superfluo. La guida v14 avverte però che nel tutorial pitzDailySteady la convergenza di SIMPLEC è molto sensibile ai fattori di rilassamento, che lì valgono 0,9 e sono stati tarati con cura: non sono adatti al SIMPLE standard, e SIMPLEC è in genere affidabile solo su geometrie semplici.

All’inizio di ogni iterazione stazionaria il modulo confronta il residuo iniziale della prima soluzione di ciascun campo elencato in residualControl con la soglia indicata (per un vettore come U si considerano le componenti risolte). Quando tutti i campi elencati sono sotto soglia, la simulazione scrive i risultati e termina prima di endTime, con un messaggio del tipo SIMPLE solution converged in 285 iterations (è il numero riportato dalla guida v14 per pitzDailySteady). I campi non elencati non vengono controllati. Se le soglie non vengono mai raggiunte, l’esecuzione prosegue fino a endTime di controlDict (Parte 4), dove nei casi stazionari il «tempo» è solo un contatore di iterazioni.

Attenzione alla coerenza tra i due livelli: il residuo iniziale di un’iterazione non può scendere molto sotto la riduzione che il solutore lineare produce a ogni iterazione. Con tolerance 1e-7 per p e residualControl p 1e-5 la combinazione è sensata; con una tolerance lineare più larga della soglia di residualControl la soglia potrebbe non essere mai raggiunta.

Con nOuterCorrectors 1 PIMPLE si riduce all’algoritmo PISO: un predittore della quantità di moto e nCorrectors correzioni di pressione per passo, senza riassemblare l’equazione della quantità di moto. È accurato ed economico se il numero di Courant massimo resta indicativamente sotto 1 (vedi il controllo di deltaT nella Parte 4). Secondo la guida v14 nCorrectors vale tipicamente 2 o 3 e nOuterCorrectors tipicamente 1.

Con nOuterCorrectors > 1 l’intero sistema viene riassemblato e risolto più volte nello stesso passo: le non linearità (il flusso convettivo phi che compare nella matrice di U) vengono aggiornate, e diventano possibili passi temporali con Courant ben superiore a 1. In questo regime le iterazioni esterne si comportano come mini-iterazioni stazionarie e può servire il rilassamento anche nel transitorio. Per non sprecare cicli esterni quando il passo è già convergente, l’implementazione Foundation prevede un sotto-dizionario outerCorrectorResidualControl (soglie per campo che interrompono in anticipo i cicli esterni): prima di usarlo, cercane un esempio nei tutorial con grep -r outerCorrectorResidualControl $FOAM_TUTORIALS.

momentumPredictor no; salta la soluzione dell’equazione della quantità di moto e passa subito alla pressione, usando la velocità del passo precedente per costruire HbyA. La guida v14 indica che si può disattivare per alcuni flussi, per esempio a basso numero di Reynolds e multifase. Nel log la conseguenza è immediata: spariscono le righe Solving for Ux/Uy/Uz.

Prendiamo un passo PIMPLE con nOuterCorrectors 2, nCorrectors 2, nNonOrthogonalCorrectors 1 e momentumPredictor yes, in 3D. Il numero di sistemi lineari risolti è:

Se ogni soluzione di p costa in media 5 iterazioni GAMG e quella finale 12, il passo richiede circa 7 × 5 + 12 = 47 cicli GAMG; portando nOuterCorrectors a 1 le soluzioni di pressione scendono a 4 e il lavoro circa si dimezza, mentre portando relTol di p da 0,01 a 0,001 resterebbero 8 soluzioni ma ognuna più lunga. Sono due leve diverse su due livelli diversi: per sapere quale usare bisogna prima leggere nel log quale livello pesa di più (sezione 7).

Il flusso diffusivo attraverso una faccia richiede il gradiente normale alla faccia, (∇p)_f · S_f. Su una mesh ortogonale il vettore d che unisce i centri delle due celle è parallelo alla normale S_f e il gradiente si ottiene esattamente dalla differenza (p_N − p_P)/|d|, una parte implicita che entra nella matrice. Se d e S_f formano un angolo θ (la non ortogonalità), si scompone S_f = Δ + k con Δ parallelo a d: la parte Δ resta implicita, la parte k · (∇p)_f va calcolata con il gradiente della soluzione corrente ed entra in b come correzione esplicita. Ogni correttore non ortogonale risolve di nuovo l’equazione della pressione con il gradiente aggiornato.

Su una mesh blockMesh a blocco unico come il canale della Parte 2 θ = 0 ovunque e i correttori non servono: 0 è la scelta giusta. checkMesh riporta Mesh non-orthogonality Max: … average: … e segnala le facce oltre 70° come severely non-orthogonal. La guida v14 indica valori tipici di 0 per i casi stazionari e 1 per quelli transitori; su mesh molto distorte si combinano 1–2 correttori con schemi laplaciani limitati (Gauss linear limited corrected 0.33, Parte 5). Nei casi stazionari la non ortogonalità viene comunque «assorbita» dalle iterazioni SIMPLE successive, ed è per questo che vi basta 0.

• Le voci vengono lette solo se il campo della pressione ha bisogno di un riferimento, cioè se nessuna patch impone un valore (nessun fixedValue o tipo equivalente). Nel canale di riferimento, con fixedValue 0 all’uscita (Parte 3), sarebbero ignorate anche se presenti.
• pRefCell è l’indice di una cella (0 è la prima cella della numerazione della mesh), pRefValue il valore di pressione attribuito lì, nelle unità del campo p (m²/s² per incompressibleFluid). La guida v14 li colloca abitualmente nel sotto-dizionario SIMPLE o PIMPLE.
• L’esempio classico è la cavità con coperchio mobile: tutte pareti, zeroGradient per p ovunque, e quindi pRefCell 0; pRefValue 0;.
• Se servono e mancano, il modulo si ferma con un errore che segnala l’assenza della voce pRefCell (o pRefPoint, l’alternativa che indica una posizione invece di un indice) nel dizionario dell’algoritmo. Senza alcun riferimento la matrice della pressione sarebbe singolare: il livello di p deriverebbe e il solutore lineare faticherebbe a convergere.

• Intestazione. Manca version 2.0;, che in v14 è facoltativa; class dictionary e object fvSolution devono corrispondere al file.
• solvers. Solo p e U: è proprio ciò che serve a un caso laminar, senza campi di turbolenza. Nessuna voce Final, perché con SIMPLE non servono.
• SIMPLE. Mesh ortogonale → 0 correttori; nessun pRefCell perché l’uscita fissa la pressione; residualControl con soglie piuttosto severe (1e-5 per p, 1e-6 per U, quest’ultima uguale alla tolerance lineare di U). Se le soglie non vengono raggiunte, l’esecuzione arriva fino a endTime: è una situazione da riconoscere nel log, non necessariamente un errore.
• relaxationFactors. La forma compatta su una riga (fields { p 0.3; }) è sintatticamente identica a quella su più righe della sezione 13.

Prima di usarlo sul canale verifica la natura del caso: foamDictionary -entry ddtSchemes/default -value system/fvSchemes deve restituire steadyState. Se restituisce Euler o backward, il caso è transitorio e serve il File B della sezione successiva. Controlla anche constant/momentumTransport: con simulationType laminar; il File A è completo così com’è.

Mescolare i due file produce errori tipici: un SIMPLE in un caso con Euler non viene letto come algoritmo transitorio; un caso transitorio senza pFinal si ferma all’ultima correzione del primo passo; un rilassamento 0,3/0,7 in un caso transitorio con nOuterCorrectors 1 rende ogni passo temporale inconsistente, perché la soluzione di fine passo non soddisfa più le equazioni. Nel File B momentumPredictor non compare: vale il predefinito, cioè il predittore viene risolto.

Sul File didattico A le due interrogazioni con -value restituiscono solo il valore, senza punto e virgola: GAMG e 0.01. Senza -value si ottiene la voce completa (solver GAMG;). Altre forme utili, tutte documentate nella §4.7 della guida v14:

foamDictionary -entry solvers -keywords system/fvSolution   # p, U (or p, pFinal, "(U|k|epsilon)", ...)
foamDictionary -expand system/fvSolution                   # macros such as $p expanded
foamDictionary -entry solvers/p/relTol -set 0.001 system/fvSolution   # edit ONLY on a copy
diff -u ../step06_steady/system/fvSolution system/fvSolution          # record exactly what changed

foamSearch -c conta quante volte compare ciascun valore della voce cercata in tutti i file con quel nome sotto la cartella indicata. La guida v14 riporta come esempio 64 occorrenze di GAMG, 26 di PCG, 5 di smoothSolver e 2 di PBiCGStab per solvers/p/solver: i numeri della tua installazione possono differire, ma l’indicazione resta che GAMG è la scelta più comune per la pressione. La stessa tecnica funziona per qualsiasi voce, per esempio foamSearch -c "$FOAM_TUTORIALS" fvSolution SIMPLE/consistent per scoprire quali tutorial usano SIMPLEC.

Una possibile sequenza completa per i passi 3 e 4, con il percorso dei tutorial della guida v14 (adattalo a quello trovato con find). Il nome step06_relTol per il secondo clone è una scelta arbitraria.

# --- step 3: baseline
cd "$HOME/OpenFOAM/training_step06/step06_steady"
foamRun > log.foamRun 2>&1
grep "solution converged" log.foamRun            # iterations to residualControl, if reached
grep -c "Solving for p" log.foamRun              # number of pressure solves
grep "Solving for p" log.foamRun | awk '{s += $NF} END {print s}'   # total p iterations
foamDictionary -entry solvers/p/relTol -value system/fvSolution     # note the current value

# --- step 4: second pristine clone, one change only
test -e "$HOME/OpenFOAM/training_step06/step06_relTol" && echo "exists: choose another name"
cp -a "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady" \
   "$HOME/OpenFOAM/training_step06/step06_relTol"
cd "$HOME/OpenFOAM/training_step06/step06_relTol"
foamDictionary -entry solvers/p/relTol -set 0.01 system/fvSolution  # choose a value different from the baseline
diff -u ../step06_steady/system/fvSolution system/fvSolution        # must show ONE changed line
blockMesh > log.blockMesh 2>&1
foamRun > log.foamRun 2>&1

Raccogli i risultati in una tabella con tre righe di confronto (iterazioni SIMPLE fino alla convergenza, iterazioni lineari totali di p, tempo di esecuzione finale ExecutionTime) e, separatamente, una grandezza fisica di controllo, per esempio la differenza di pressione tra ingresso e uscita. Una conclusione motivata ha la forma: «con relTol più stretto ogni soluzione di p richiede in media N iterazioni in più, il numero di iterazioni SIMPLE cambia di M e la grandezza di controllo varia meno di X%: il costo aumenta senza un cambiamento apprezzabile della soluzione». Notare che nei casi stazionari un relTol più stretto può anche ridurre il numero di iterazioni esterne: è proprio il tipo di effetto che il confronto deve misurare, non presumere.

Gli errori fatali di OpenFOAM hanno sempre la stessa struttura: una riga --> FOAM FATAL ERROR o FOAM FATAL IO ERROR, il messaggio, poi il file con le righe coinvolte e la funzione che ha generato l’errore. Per isolarlo da un log lungo:

grep -n "FOAM FATAL" log.foamRun
grep -n -A 15 "FOAM FATAL" log.foamRun | head -40

L’Esempio 2 genera un errore fatale esplicito (sezione 10) alla prima soluzione di U, quindi è facile da diagnosticare. L’Esempio 1 è più insidioso: la User Guide v14 indica tolerance e relTol come richieste, ma a seconda dell’implementazione una voce omessa può non essere segnalata affatto e venire presa silenziosamente con un valore predefinito (per relTol lo zero, cioè nessun criterio relativo). In quel caso ogni soluzione di p andrebbe fino a tolerance e l’unico sintomo sarebbe un numero di iterazioni più alto del previsto nel log. Per questo la verifica corretta non è «il caso parte?», ma leggere il dizionario espanso con foamDictionary -expand -entry solvers/p system/fvSolution e confrontarlo con la documentazione.

I comandi seguenti svolgono la parte di raccolta dati senza modificare nulla. pitzDailySteady (stazionario) e pitzDaily (transitorio) sono le due coppie naturali, entrambe citate dalla guida v14; il percorso va confermato come nel Laboratorio A.

ls "$FOAM_TUTORIALS/incompressibleFluid"
for c in pitzDailySteady pitzDaily; do
    f="$FOAM_TUTORIALS/incompressibleFluid/$c/system/fvSolution"
    echo "== $c"
    foamDictionary -entry solvers -keywords "$f"
    foamDictionary -entry ddtSchemes/default -value "${f%fvSolution}fvSchemes"
    foamDictionary -entry SIMPLE "$f" 2>/dev/null || foamDictionary -entry PIMPLE "$f"
    foamDictionary -entry relaxationFactors "$f"
done

Per ogni tutorial compila una scheda con: elenco dei blocchi in solvers e, per ciascuno, solver, preconditioner/smoother, tolerance, relTol; quali nomi reali coprono le regex (confrontandoli con il modello in constant/momentumTransport); algoritmo e controlli (nCorrectors, nOuterCorrectors, nNonOrthogonalCorrectors, consistent, residualControl, pRefCell/pRefValue); fattori di rilassamento. Per l’ultimo compito usa la forma «modifica → livello interessato (sezione 16) → effetto previsto sulle iterazioni → come lo misuro nel log (sezione 7)».