OpenFOAM 14 Part 6: fvSolution, Linear Solvers and Solution Algorithms

Up to version 10 of the Foundation distribution each combination of physics and algorithm had a dedicated executable: simpleFoam for steady incompressible flows with SIMPLE, pisoFoam and pimpleFoam for transient ones. From version 11 these executables were replaced by solver modules loaded by a single program, foamRun; the module is selected in system/controlDict with the entry solver incompressibleFluid; (or with foamRun -solver incompressibleFluid).

The consequence for fvSolution matters: in v14 it is no longer the executable name that says whether the algorithm is SIMPLE or PIMPLE. The same incompressibleFluid module handles both cases: steady behaviour is obtained by declaring ddtSchemes { default steadyState; } in fvSchemes and writing a SIMPLE sub-dictionary in fvSolution, as in the v14 guide’s pitzDailySteady tutorial. The linear-solver entries (solvers), on the other hand, have stayed substantially the same: that is why an old fvSolution is often still readable, but its context (executable, model files, dimension names) must be checked.

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

According to the v14 User Guide (§4.6), fvSolution contains a solvers sub-dictionary with the linear-solver settings of each equation, a relaxationFactors sub-dictionary and an algorithm-control sub-dictionary, SIMPLE for steady cases or PIMPLE for transient ones. Other modules may read further entries (for instance multiphase modules read phase-fraction-specific controls), but for incompressibleFluid these three blocks are all that is needed.

The file is not read only once. A linear solver’s settings are looked up in the solvers dictionary every time an equation is solved, and with runTimeModifiable active (the usual behaviour, see Part 4) OpenFOAM re-reads the system/ files when they are modified during the run. You can therefore, for example, lower a relaxation factor while the calculation is running; the change takes effect at the next iteration and the log reports the re-reading of the dictionary. It is convenient, but it makes the history of the simulation less reproducible: always record the change and the time at which you made it.

Typical error: if the block for a field that the module must solve is missing, the run stops with a FOAM FATAL IO ERROR of the form “keyword p is undefined in dictionary .../system/fvSolution/solvers”. The message gives both the name being looked up and the dictionary path: it is the first thing to read.

Consider steady one-dimensional diffusion d/dx(Γ dφ/dx) = 0 on a segment of length 3 divided into three equal cells (Δx = 1, Γ = 1), with φ = 1 at the left boundary and φ = 0 at the right boundary. Integrating over cell P and approximating the diffusive flux through each face by the difference between centres, Γ(φ_N − φ_P)/d, gives one equation per cell. For an internal face d = Δx = 1 and the coefficient is 1; for a boundary face the cell centre is d = Δx/2 from the face and the coefficient is 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 |

The solution, x = (0.833, 0.5, 0.167), coincides with the exact linear profile 1 − x/3 evaluated at the cell centres (0.5, 1.5, 2.5). In the notation used here, a_P is the diagonal term (3, 2, 3), the a_N are the off-diagonal coefficients (−1) and b_P collects the known contributions, here the boundary values. The magnitude of the diagonal is greater than or equal to the sum of the off-diagonal magnitudes in every row: the matrix is diagonally dominant, and it is symmetric because the coefficient linking P to N is the same as the one linking N to P. With one million hexahedral cells the structure is identical: each row has one diagonal entry plus one per internal face of the cell (at most six), so about 7 million non-zeros versus the 10¹² entries of a dense matrix.

OpenFOAM does not use a generic sparse format (such as CSR) but the LDU format of the lduMatrix class, built directly on the mesh topology: a diag array with one value per cell, and an upper and a lower array with one value per internal face. Each internal face has an owner cell (lower index) and a neighbour cell: the face’s upper coefficient sits in the owner’s row and the neighbour’s column, the lower coefficient in the transposed position. If the matrix is symmetric, lower is not even allocated.

Boundary faces do not generate off-diagonal entries: the boundary condition (Part 3) contributes an implicit part to the diagonal and an explicit part to the source b. That is why a fixedValue strengthens the diagonal of the adjacent cell (like the 3 in the three-cell example), while a zeroGradient adds nothing. Processor patches in parallel runs (see the worked case on MPI parallel execution) are treated as “interfaces” that exchange values between subdomains at every matrix–vector product: this is why a linear solver’s iteration count can change slightly with the number of processors.

The “may” in the prediction check is precise. In OpenFOAM many higher-order schemes are implemented as a deferred correction: the implicit part, which goes into the matrix, is that of the upwind scheme, while the difference between the chosen scheme and upwind is computed from the current values and enters the source b. This is the case for linearUpwind: switching from bounded Gauss upwind to bounded Gauss linearUpwind grad(U) leaves the momentum matrix A structurally that of upwind and changes mainly b. Switching instead to Gauss linear (central differencing) changes the interpolation weights and therefore the coefficients of A: diagonal dominance weakens and iterative solvers may need more iterations or become unstable on coarse meshes.

The second statement holds both ways: a different linear solver, if it converges to the same tolerance, produces the same discrete solution apart from the remaining algebraic error. If changing only the linear solver noticeably changes the final result (for example the pressure drop), the tolerances were too loose: it is a useful test to check this.

The table above lists the essential entries. The v14 User Guide documents further ones, with their defaults:

The diagonal solver is for explicit equations, whose matrix is diagonal only: it does not iterate and uses no tolerances. To see which names your installation accepts you can deliberately provoke an error on a copy (for example solver xyz;): OpenFOAM replies with the list of valid solvers for that matrix type. To check what the module will actually read after macro expansion, use foamDictionary -expand -entry solvers system/fvSolution.

When the module has to solve, for example, the k equation, it looks up the entry k in the solvers dictionary with this precedence: first a key with the exact name; if there is none, the quoted keys that are regular expressions, examined from last to first in file order. The expression must match the whole name: "(U|k|epsilon)" matches k but not kFinal, which is why the second group "(U|k|epsilon)Final" is needed.

The macro $p; copies into the current block all the entries of the p block; entries written after the macro override them (here relTol 0;). $U; also works although there is no literal key U: the macro lookup accepts regex keys, so it finds "(U|k|epsilon)". The expanded result can be checked like this:

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

In the PIMPLE block of this extract nOuterCorrectors does not appear: the default 1 applies, i.e. “PISO” mode (a single outer pass per time step with two pressure corrections), suitable when the Courant number stays below 1. In a laminar case such as the reference channel k and epsilon do not exist: the regex group is harmless but useless, and it is better to reduce it to U and UFinal.

The “raw” residual of a system is the vector b − A x, which has the units of the equation and depends on cell size: its absolute value would mean nothing. OpenFOAM makes it dimensionless by dividing it by a normalisation factor built with the mean value x̄ of the current field:

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

The denominator measures how far the problem “departs” from a uniform solution. Two practical consequences. First: if the initial field is uniform (x = x̄, like internalField uniform 0), the first term of the denominator vanishes and the second equals the numerator, so the initial residual of the first step is 1: it is not an error, it is the normalisation. Second: the residual is a sum (L1 norm) over all cells divided by a similar sum, so it is an average; small residuals can hide a few cells with large errors, one more reason for the physical checks of the following sections.

On maxIter: if the log shows No Iterations 1000 (or the value you set), the solver reached neither tolerance nor relTol and was stopped by the cap. It is not an error message, but it is almost always the symptom of a problem: ill-conditioned matrix, unsuitable solver, poor mesh quality or a solution that is diverging.

Many iterative methods, after the first iterations, reduce the residual by a roughly constant factor ρ per iteration (the convergence factor). In the table the first step takes the residual from 2.0 × 10⁻³ to 8.0 × 10⁻⁴, i.e. ρ ≈ 0.4. If ρ stayed at 0.4, after n iterations the ratio would be 0.4ⁿ, and the iterations needed for each criterion can be estimated:

The estimate shows the general rule: each extra decade of reduction costs about ln 10 / |ln ρ| iterations (here about 2.5). The actual value in the table falls faster (1.5 × 10⁻⁵ corresponds to a ratio of 0.0075), because ρ is not constant: with GAMG in particular each cycle reduces the residual far more than a factor 0.4. The correct way to measure ρ on your own case is to read the initial residual, final residual and iteration count in the log: ρ ≈ (r_final/r_initial)^1/N. With the log lines of section 7: for p, (0.00008/0.014)^1/5 ≈ 0.36.

Here is what one step of a 3D transient incompressibleFluid case with PIMPLE in PISO mode (nOuterCorrectors 1, nCorrectors 2) typically looks like. The numbers are illustrative; the structure is what you will find in your 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

• The three lines Ux, Uy, Uz are the momentum predictor: the vector equation is solved component by component with the same matrix. In a 2D case with empty patches the component normal to the plane is not solved and its line is missing.
• The two p lines are the two pressure corrections (nCorrectors 2); the second uses the pFinal settings (section 8), so it goes down to tolerance and needs more iterations.
• time step continuity errors is the mass-balance indicator called for by the “Physical checks” box: sum local is the weighted mean of the absolute flux divergence, global its signed sum, cumulative the sum over time. Values that grow step after step indicate that pressure is not solved well enough.
• With a preconditioner its name appears as a prefix: DICPCG: Solving for p, DILUPBiCGStab: Solving for Ux.

To compare two runs (as Lab A requires) you need numbers, not impressions. A few ways, from the simplest to the most structured:

# 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

Alternatively the residuals function writes the initial residuals to a text file during the run: add #includeFunc residuals(p, U) to the functions sub-dictionary of controlDict (Part 4) and you will find the data in postProcessing/residuals/; check its options with foamInfo residuals. Remember that these files record the initial residual of each solve: it is the “across steps” history, not the “inside a solve” one.

In the Foundation implementation the suffix Final is appended to the field name when the algorithm is in its last iteration: for velocity and the other fields during the last PIMPLE outer loop (with nOuterCorrectors 1, therefore, at every step); for pressure in the last pressure correction of the last outer loop (and, if there are non-orthogonal correctors, in the last of them). The v14 User Guide describes the purpose: solving the equation with separate settings for the final corrector, typically relTol 0, so that each time step converges down to tolerance.

Why it matters in transient cases: the algebraic error left at the end of a step is never corrected again and accumulates over time as continuity error. Intermediate corrections, by contrast, will be superseded by the next one anyway: there is no point paying a high price for them. In the log of section 7 the difference is visible: 5 iterations with relTol 0.01, 9 with relTol 0.

If the entry is missing in a transient case, the run stops at the first final solve with an error of the form “keyword pFinal is undefined in dictionary .../system/fvSolution/solvers”. In steady cases with SIMPLE (like Teaching file A in section 18 and the pitzDailySteady tutorial) Final entries normally do not appear.

• Double quotes are required: the quoting is what tells OpenFOAM to treat the key as a regular expression. Without quotes, (U|k|epsilon) is not even a valid word.
• Matching is on the whole name: "U.*" catches U, UFinal but also, for example, any Ua; ".*" catches everything and is legitimate only in contexts such as relaxationFactors (section 19).
• The dot is a metacharacter: in quotes, "alpha.water" also matches alphaXwater. When the literal dot matters, the safest solution is to write the exact unquoted name (alpha.water), which is compared character by character.
• Precedence: the exact name always wins; among several regexes the last one in the file wins. An explicit U block written after the group therefore takes precedence over the group for U.

Which fields must be covered depends on the model in constant/momentumTransport (Part 7):

nut, the turbulent viscosity, is computed by the model and has no equation of its own: it needs no solvers entry. If a solved field is not covered, the error appears at that field’s first solve, with the message “keyword omega is undefined in dictionary …/solvers”.

Take the three cells of section 2 again and add a convective flux F = 1 towards +x, discretised with upwind: on each cell’s outflow face the transported value is the cell’s own (contribution +F on the diagonal), on the inflow face it is that of the upstream cell (contribution −F on the upstream neighbour’s coefficient). The row of cell 2 becomes (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

The momentum equation always contains the convective term div(phi,U), so its matrix is asymmetric; the pressure equation of the SIMPLE/PIMPLE algorithms is instead a Laplace (or Poisson) equation ∇·((1/A)∇p) = ∇·(HbyA), with a symmetric matrix. That is why in the tutorials p uses GAMG or PCG, while U, k, epsilon use smoothSolver or PBiCGStab.

A clarification on the table: the v14 User Guide pairs PCG with symmetric matrices and PBiCG with asymmetric ones, but presents PBiCGStab as a stabilised variant usable for both classes; indeed the foamSearch output shown in the guide includes tutorials that use it for p as well. The table’s “asymmetric” row should be read as “suitable for asymmetric too”. GAMG and smoothSolver accept both classes, provided a compatible smoother is chosen (section 11).

According to the v14 guide, using a symmetric solver on an asymmetric matrix (or vice versa) produces a FOAM FATAL IO ERROR. The typical form is the following; the exact list depends on the version and the libraries loaded:

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

Read it as the box teaches: the wrong keyword is solver, the file and block are given in the last line, the matrix class (“asymmetric”) is written in the message and the valid alternatives are listed. There is only one fix: change solver (and therefore preconditioner, because DIC is for symmetric matrices only: with PBiCGStab use DILU).

Gauss-Seidel. The cells are visited in order and each is updated by solving its own equation with the latest neighbour values: φ_P ← (b_P − Σ a_N φ_N) / a_P. Each update propagates information by one cell only, so errors that oscillate from cell to cell (high frequency) vanish within a few sweeps, while those spread over the whole domain decay very slowly. symGaussSeidel performs a forward and a backward sweep, removing the directional bias tied to cell numbering: it is the preferred smoother in the tutorials for U and turbulence fields. The nSweeps entry (default 1) sets how many sweeps are made before the residual is recalculated.

DIC and DILU. They are incomplete factorisations: an approximation M ≈ A = L D U is built with the same sparsity as A (no fill-in) and in which, in the “diagonal” variant used by OpenFOAM, only the diagonal is modified. Solving with M costs about two substitution sweeps, and the Krylov method works on the preconditioned system M⁻¹A x = M⁻¹b, which is much better conditioned. DIC (incomplete Cholesky) requires a symmetric matrix; DILU (incomplete LU) is the version for asymmetric matrices; FDIC is a faster version of DIC that precomputes reciprocal quantities. The same objects can act as smoothers.

DICGaussSeidel applies DIC smoothing followed by Gauss-Seidel and is the classic choice as the GAMG smoother for pressure: being DIC-based, it is for symmetric matrices (for an asymmetric equation solved with GAMG use GaussSeidel or symGaussSeidel). The v14 guide also lists GAMG itself among the preconditioners (a multigrid cycle preconditioning PCG), diagonal (generally not used) and none.

Decompose the error e = x − x_exact into Fourier modes. A Gauss-Seidel iteration strongly damps modes whose wavelength is a few cells, but a mode spread over n cells is reduced only by a factor of about 1 − c/n² per sweep: on a mesh of 200 cells per side tens of thousands of sweeps would be needed. The idea of multigrid is that an error that is smooth on the fine grid becomes oscillatory on a coarser grid: after a few smoothing sweeps (pre-sweeps) the residual is transferred to a level with fewer cells, where the same smoother is effective again; going down to the coarsest level (about nCellsInCoarsestLevel cells) the system is so small it is solved almost exactly; the correction is then interpolated upwards and refined with further sweeps (post-sweeps). The result is a cycle count almost independent of mesh size and a cost proportional to the number of cells.

Comparison with PCG: the conjugate-gradient iteration count grows as the square root of the condition number, which for a Laplacian grows as (L/h)²; halving the cell size therefore roughly doubles the PCG iterations, while those of GAMG stay almost the same. That is why GAMG dominates for pressure on medium-to-large meshes, whereas on small meshes (such as the teaching channel with a few thousand cells) the difference is minimal. “Algebraic” means that the coarse levels are not geometrically built meshes but cell agglomerations decided from the matrix coefficients (in the tutorials with the faceAreaPair agglomerator, which pairs cells across their largest faces).

A note on reading the guide: the v14 User Guide writes this entry as directSolveCoarset; the name used in the code and in dictionaries is directSolveCoarsest, but before using it check its spelling in the tutorials with foamSearch. For the first experiment the advice of the box still holds: these defaults are fine in the vast majority of cases.

Implicit equation relaxation (Patankar style) divides the diagonal by α and adds to the source a term containing the previous value:

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

At convergence φ_P = φ_P,old and the two added terms cancel: the final solution is the same, only the path changes. With α = 0.7 the diagonal grows by 1/0.7 ≈ 1.43 times: the matrix becomes more diagonally dominant (iterative solvers converge better) and each iteration takes a shorter step. In OpenFOAM, before α is applied, the diagonal is also raised to at least the sum of the magnitudes of the off-diagonal coefficients: that is why the v14 guide indicates for transient cases a factor 1 for all equations (equations { ".*" 1; }, as in Teaching file B), which slows nothing down but ensures what the guide calls “diagonal equality”, generally required for the convergence of a matrix equation.

Field example. If in a cell p_old = 10 m²/s² and the pressure equation proposes p_candidate = 20 m²/s², with α = 0.3 the accepted value is 10 + 0.3 × 10 = 13 m²/s². Indications from the v14 guide: steady tutorials use 0.3 for p and 0.7 for U and for turbulence fields (grouped with a regex); values as high as 0.9 can be stable in some cases, while values well below 0.2 make convergence prohibitively slow. A field or equation absent from relaxationFactors is not relaxed.

Write the discrete momentum equation of each cell separating the diagonal A from the rest: A U_P = H(U) − ∇p, where H collects the neighbour contributions and sources. Hence U = HbyA − (1/A)∇p, with HbyA = H/A. Imposing continuity ∇·U = 0 gives the pressure equation. One SIMPLE iteration of incompressibleFluid proceeds as follows:

• Momentum predictor: the U equation is assembled with the current pressure, relaxed (factor equations/U) and solved with the solvers/U settings (lines Ux, Uy, Uz in the log).
• Pressure equation: ∇·((1/A)∇p) = ∇·(HbyA), solved with solvers/p, repeated nNonOrthogonalCorrectors more times if requested.
• Flux correction: the face flux phi is corrected with the pressure gradient just computed, and becomes conservative up to the pressure tolerance.
• Pressure relaxation (factor fields/p, e.g. 0.3) and velocity correction U = HbyA − (1/A)∇p.
• Model equations (turbulence etc.), then the next iteration.

In the SIMPLEC formulation (consistent yes;) 1/(A − ΣH₁) is used instead of 1/A, where ΣH₁ is the sum of the neighbour coefficients: the response of the neighbouring cells is approximately accounted for, the pressure correction is more “consistent” and pressure relaxation becomes largely unnecessary. The v14 guide warns, however, that in the pitzDailySteady tutorial SIMPLEC convergence is very sensitive to the relaxation factors, which there are 0.9 and carefully tuned: they are not suitable for standard SIMPLE, and SIMPLEC is generally reliable only on simple geometries.

At the start of each steady iteration the module compares the initial residual of the first solve of each field listed in residualControl with the given threshold (for a vector such as U the solved components are considered). When all listed fields are below threshold, the simulation writes the results and ends before endTime, with a message such as SIMPLE solution converged in 285 iterations (the number reported by the v14 guide for pitzDailySteady). Fields not listed are not checked. If the thresholds are never reached, the run continues up to the controlDict endTime (Part 4), where in steady cases “time” is only an iteration counter.

Mind the consistency between the two levels: the initial residual of an iteration cannot drop much below the reduction that the linear solver achieves at each iteration. With tolerance 1e-7 for p and residualControl p 1e-5 the combination is sensible; with a linear tolerance looser than the residualControl threshold, the threshold might never be reached.

With nOuterCorrectors 1 PIMPLE reduces to the PISO algorithm: one momentum predictor and nCorrectors pressure corrections per step, without reassembling the momentum equation. It is accurate and cheap if the maximum Courant number stays roughly below 1 (see the deltaT control in Part 4). According to the v14 guide nCorrectors is typically 2 or 3 and nOuterCorrectors typically 1.

With nOuterCorrectors > 1 the whole system is reassembled and solved several times within the same step: the non-linearities (the convective flux phi appearing in the U matrix) are updated, and time steps with a Courant number well above 1 become possible. In this regime the outer iterations behave like mini steady iterations and relaxation may be needed even in transient runs. To avoid wasting outer loops when the step has already converged, the Foundation implementation provides an outerCorrectorResidualControl sub-dictionary (per-field thresholds that end the outer loops early): before using it, look for an example in the tutorials with grep -r outerCorrectorResidualControl $FOAM_TUTORIALS.

momentumPredictor no; skips the solution of the momentum equation and goes straight to pressure, using the previous step’s velocity to build HbyA. The v14 guide indicates it can be switched off for some flows, for example at low Reynolds number and multiphase. The consequence is immediately visible in the log: the Solving for Ux/Uy/Uz lines disappear.

Take a PIMPLE step with nOuterCorrectors 2, nCorrectors 2, nNonOrthogonalCorrectors 1 and momentumPredictor yes, in 3D. The number of linear systems solved is:

If each p solve costs on average 5 GAMG iterations and the final one 12, the step needs about 7 × 5 + 12 = 47 GAMG cycles; setting nOuterCorrectors to 1 brings the pressure solves down to 4 and roughly halves the work, whereas tightening the p relTol from 0.01 to 0.001 would keep 8 solves but make each longer. They are two different levers on two different levels: to know which to use you must first read in the log which level weighs most (section 7).

The diffusive flux through a face requires the face-normal gradient, (∇p)_f · S_f. On an orthogonal mesh the vector d joining the two cell centres is parallel to the normal S_f and the gradient is obtained exactly from the difference (p_N − p_P)/|d|, an implicit part that enters the matrix. If d and S_f form an angle θ (the non-orthogonality), S_f is split as Δ + k with Δ parallel to d: the Δ part stays implicit, the k · (∇p)_f part must be computed from the gradient of the current solution and enters b as an explicit correction. Each non-orthogonal corrector solves the pressure equation again with the updated gradient.

On a single-block blockMesh mesh such as the channel of Part 2 θ = 0 everywhere and no correctors are needed: 0 is the right choice. checkMesh reports Mesh non-orthogonality Max: … average: … and flags faces beyond 70° as severely non-orthogonal. The v14 guide indicates typical values of 0 for steady cases and 1 for transient ones; on highly distorted meshes 1–2 correctors are combined with limited Laplacian schemes (Gauss linear limited corrected 0.33, Part 5). In steady cases the non-orthogonality is in any case “absorbed” by the subsequent SIMPLE iterations, which is why 0 is enough there.

• The entries are read only if the pressure field needs a reference, i.e. if no patch imposes a value (no fixedValue or equivalent type). In the reference channel, with fixedValue 0 at the outlet (Part 3), they would be ignored even if present.
• pRefCell is the index of a cell (0 is the first cell of the mesh numbering), pRefValue the pressure value assigned there, in the units of the p field (m²/s² for incompressibleFluid). The v14 guide usually places them in the SIMPLE or PIMPLE sub-dictionary.
• The classic example is the lid-driven cavity: walls only, zeroGradient for p everywhere, hence pRefCell 0; pRefValue 0;.
• If they are needed and missing, the module stops with an error reporting the absence of the pRefCell entry (or pRefPoint, the alternative giving a location instead of an index) in the algorithm dictionary. Without any reference the pressure matrix would be singular: the level of p would drift and the linear solver would struggle to converge.

• Header. version 2.0; is missing, which is optional in v14; class dictionary and object fvSolution must match the file.
• solvers. Only p and U: exactly what a laminar case needs, with no turbulence fields. No Final entries, because with SIMPLE they are not needed.
• SIMPLE. Orthogonal mesh → 0 correctors; no pRefCell because the outlet fixes the pressure; residualControl with fairly strict thresholds (1e-5 for p, 1e-6 for U, the latter equal to the linear tolerance of U). If the thresholds are not reached, the run goes on to endTime: a situation to recognise in the log, not necessarily an error.
• relaxationFactors. The compact one-line form (fields { p 0.3; }) is syntactically identical to the multi-line form of section 13.

Before using it on the channel check the nature of the case: foamDictionary -entry ddtSchemes/default -value system/fvSchemes must return steadyState. If it returns Euler or backward, the case is transient and File B of the next section is needed. Also check constant/momentumTransport: with simulationType laminar; File A is complete as it is.

Mixing the two files produces typical errors: a SIMPLE dictionary in a case with Euler is not read as the transient algorithm; a transient case without pFinal stops at the last correction of the first step; 0.3/0.7 relaxation in a transient case with nOuterCorrectors 1 makes each time step inconsistent, because the end-of-step solution no longer satisfies the equations. In File B momentumPredictor does not appear: the default applies, i.e. the predictor is solved.

On Teaching file A the two -value queries return the value only, without semicolon: GAMG and 0.01. Without -value you get the full entry (solver GAMG;). Other useful forms, all documented in §4.7 of the v14 guide:

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 counts how many times each value of the searched entry appears in all files with that name below the given directory. The v14 guide shows as an example 64 occurrences of GAMG, 26 of PCG, 5 of smoothSolver and 2 of PBiCGStab for solvers/p/solver: the numbers in your installation may differ, but the indication remains that GAMG is the most common choice for pressure. The same technique works for any entry, for example foamSearch -c "$FOAM_TUTORIALS" fvSolution SIMPLE/consistent to find out which tutorials use SIMPLEC.

A possible complete sequence for steps 3 and 4, with the tutorial path of the v14 guide (adapt it to the one found with find). The name step06_relTol for the second clone is an arbitrary choice.

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

Collect the results in a table with three comparison rows (SIMPLE iterations to convergence, total p linear iterations, final ExecutionTime) and, separately, a physical control quantity, for example the pressure difference between inlet and outlet. A justified conclusion has the form: “with a tighter relTol each p solve needs on average N more iterations, the number of SIMPLE iterations changes by M and the control quantity varies by less than X%: the cost increases without an appreciable change in the solution”. Note that in steady cases a tighter relTol may even reduce the number of outer iterations: this is exactly the kind of effect the comparison must measure, not assume.

OpenFOAM fatal errors always have the same structure: a line --> FOAM FATAL ERROR or FOAM FATAL IO ERROR, the message, then the file with the lines involved and the function that raised the error. To isolate it in a long log:

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

Example 2 produces an explicit fatal error (section 10) at the first U solve, so it is easy to diagnose. Example 1 is more insidious: the v14 User Guide lists tolerance and relTol as required, but depending on the implementation an omitted entry may not be reported at all and may silently take a default value (for relTol, zero, i.e. no relative criterion). In that case every p solve would run down to tolerance and the only symptom would be a higher-than-expected iteration count in the log. That is why the correct check is not “does the case start?” but reading the expanded dictionary with foamDictionary -expand -entry solvers/p system/fvSolution and comparing it with the documentation.

The following commands carry out the data-gathering part without modifying anything. pitzDailySteady (steady) and pitzDaily (transient) are the natural pair, both cited by the v14 guide; the path must be confirmed as in Lab 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

For each tutorial fill in a sheet with: the list of blocks in solvers and, for each, solver, preconditioner/smoother, tolerance, relTol; which actual names the regexes cover (comparing them with the model in constant/momentumTransport); algorithm and controls (nCorrectors, nOuterCorrectors, nNonOrthogonalCorrectors, consistent, residualControl, pRefCell/pRefValue); relaxation factors. For the last task use the form “change → level affected (section 16) → predicted effect on iterations → how I measure it in the log (section 7)”.