## 4.5 Solution and algorithm control

The equation solvers, tolerances and algorithms are controlled from the fvSolution dictionary in the system directory. Below is an example set of entries from the fvSolution dictionary required for the icoFoam solver.

18 solvers

19 {

20 p

21 {

22 solver PCG;

23 preconditioner DIC;

24 tolerance 1e-06;

25 relTol 0;

26 }

27

28 U

29 {

30 solver smoothSolver;

31 smoother symGaussSeidel;

32 tolerance 1e-05;

33 relTol 0;

34 }

35 }

36

37 PISO

38 {

39 nCorrectors 2;

40 nNonOrthogonalCorrectors 0;

41 pRefCell 0;

42 pRefValue 0;

43 }

44

45

46 // ************************************************************************* //

fvSolution contains a set of subdictionaries that are specific to the solver being run. However, there is a small set of standard subdictionaries that cover most of those used by the standard solvers. These subdictionaries include solvers, relaxationFactors, PISO and SIMPLE which are described in the remainder of this section.

### 4.5.1 Linear solver control

The first sub-dictionary in our example, and one that appears in all solver applications, is solvers. It specifies each linear-solver that is used for each discretised equation; it is emphasised that the term linear-solver refers to the method of number-crunching to solve the set of linear equations, as opposed to application solver which describes the set of equations and algorithms to solve a particular problem. The term ‘linear-solver’ is abbreviated to ‘solver’ in much of the following discussion; we hope the context of the term avoids any ambiguity.

The syntax for each entry within solvers uses a keyword that is the word relating to the variable being solved in the particular equation. For example, icoFoam solves equations for velocity and pressure , hence the entries for U and p. The keyword is followed by a dictionary containing the type of solver and the parameters that the solver uses. The solver is selected through the solver keyword from the choice in OpenFOAM, listed in Table 4.12. The parameters, including tolerance, relTol, preconditioner, etc. are described in following sections.

The solvers distinguish between symmetric matrices and asymmetric matrices. The symmetry of the matrix depends on the structure of the equation being solved and, while the user may be able to determine this, it is not essential since OpenFOAM will produce an error message to advise the user if an inappropriate solver has been selected, e.g.

Valid asymmetric matrix solvers are :

3

(

PBiCG

smoothSolver

GAMG

)

#### 4.5.1.1 Solution tolerances

The sparse matrix solvers are iterative, i.e. they are based on reducing the equation residual over a succession of solutions. The residual is ostensibly a measure of the error in the solution so that the smaller it is, the more accurate the solution. More precisely, the residual is evaluated by substituting the current solution into the equation and taking the magnitude of the difference between the left and right hand sides; it is also normalised to make it independent of the scale of the problem being analysed.

Before solving an equation for a particular field, the initial residual is evaluated based on the current values of the field. After each solver iteration the residual is re-evaluated. The solver stops if either of the following conditions are reached:

- the residual falls below the solver tolerance, tolerance;
- the ratio of current to initial residuals falls below the solver relative tolerance, relTol;
- the number of iterations exceeds a maximum number of iterations, maxIter;

The solver tolerance should represent the level at which the residual is small enough that the solution can be deemed sufficiently accurate. The solver relative tolerance limits the relative improvement from initial to final solution. In transient simulations, it is usual to set the solver relative tolerance to 0 to force the solution to converge to the solver tolerance in each time step. The tolerances, tolerance and relTol must be specified in the dictionaries for all solvers; maxIter is optional.

#### 4.5.1.2 Preconditioned conjugate gradient solvers

There are a range of options for preconditioning of matrices in the conjugate gradient solvers, represented by the preconditioner keyword in the solver dictionary. The preconditioners are listed in Table 4.13.

#### 4.5.1.3 Smooth solvers

The solvers that use a smoother require the smoother to be specified. The smoother options are listed in Table 4.14. Generally GaussSeidel is the most reliable option, but for bad matrices DIC can offer better convergence. In some cases, additional post-smoothing using GaussSeidel is further beneficial, i.e. the method denoted as DICGaussSeidel

The user must also pecify the number of sweeps, by the nSweeps keyword, before the residual is recalculated, following the tolerance parameters.

#### 4.5.1.4 Geometric-algebraic multi-grid solvers

The generalised method of geometric-algebraic multi-grid (GAMG) uses the principle of: generating a quick solution on a mesh with a small number of cells; mapping this solution onto a finer mesh; using it as an initial guess to obtain an accurate solution on the fine mesh. GAMG is faster than standard methods when the increase in speed by solving first on coarser meshes outweighs the additional costs of mesh refinement and mapping of field data. In practice, GAMG starts with the mesh specified by the user and coarsens/refines the mesh in stages. The user is only required to specify an approximate mesh size at the most coarse level in terms of the number of cells nCoarsestCells.

The agglomeration of cells is performed by the algorithm specified by the agglomerator keyword. Presently we recommend the faceAreaPair method. It is worth noting there is an MGridGen option that requires an additional entry specifying the shared object library for MGridGen:

In the experience of OpenCFD, the MGridGen method offers no obvious benefit over the faceAreaPair method. For all methods, agglomeration can be optionally cached by the cacheAgglomeration switch.

Smoothing is specified by the smoother as described in section 4.5.1.3. The number of sweeps used by the smoother at different levels of mesh density are specified by the nPreSweeps, nPostSweeps and nFinestSweeps keywords. The nPreSweeps entry is used as the algorithm is coarsening the mesh, nPostSweeps is used as the algorithm is refining, and nFinestSweeps is used when the solution is at its finest level.

The mergeLevels keyword controls the speed at which coarsening or refinement levels is performed. It is often best to do so only at one level at a time, i.e. set mergeLevels 1. In some cases, particularly for simple meshes, the solution can be safely speeded up by coarsening/refining two levels at a time, i.e. setting mergeLevels 2.

### 4.5.2 Solution under-relaxation

A second sub-dictionary of fvSolution that is often used in OpenFOAM is relaxationFactors which controls under-relaxation, a technique used for improving stability of a computation, particularly in solving steady-state problems. Under-relaxation works by limiting the amount which a variable changes from one iteration to the next, either by modifying the solution matrix and source prior to solving for a field or by modifying the field directly. An under-relaxation factor specifies the amount of under-relaxation, as described below.

- No specified : no under-relaxation.
- : guaranteed matrix diagonal equality/dominance.
- decreases, under-relaxation increases.
- : solution does not change with successive iterations.

An optimum choice of is one that is small enough to ensure stable computation but large enough to move the iterative process forward quickly; values of as high as 0.9 can ensure stability in some cases and anything much below, say, 0.2 are prohibitively restrictive in slowing the iterative process.

The user can specify the relaxation factor for a particular field by specifying first the word associated with the field, then the factor. The user can view the relaxation factors used in a tutorial example of simpleFoam for incompressible, laminar, steady-state flows.

18 solvers

19 {

20 p

21 {

22 solver GAMG;

23 tolerance 1e-06;

24 relTol 0.1;

25 smoother GaussSeidel;

26 nPreSweeps 0;

27 nPostSweeps 2;

28 cacheAgglomeration on;

29 agglomerator faceAreaPair;

30 nCellsInCoarsestLevel 10;

31 mergeLevels 1;

32 }

33

34 ”(U|k|epsilon|omega|f|v2)”

35 {

36 solver smoothSolver;

37 smoother symGaussSeidel;

38 tolerance 1e-05;

39 relTol 0.1;

40 }

41 }

42

43 SIMPLE

44 {

45 nNonOrthogonalCorrectors 0;

46 consistent yes;

47

48 residualControl

49 {

50 p 1e-2;

51 U 1e-3;

52 ”(k|epsilon|omega|f|v2)” 1e-3;

53 }

54 }

55

56 relaxationFactors

57 {

58 equations

59 {

60 U 0.9; // 0.9 is more stable but 0.95 more convergent

61 ”.*” 0.9; // 0.9 is more stable but 0.95 more convergent

62 }

63 }

64

65

66 // ************************************************************************* //

### 4.5.3 PISO and SIMPLE algorithms

Most fluid dynamics solver applications in OpenFOAM use the pressure-implicit split-operator (PISO) or semi-implicit method for pressure-linked equations (SIMPLE) algorithms. These algorithms are iterative procedures for solving equations for velocity and pressure, PISO being used for transient problems and SIMPLE for steady-state.

Both algorithms are based on evaluating some initial solutions and then correcting them. SIMPLE only makes 1 correction whereas PISO requires more than 1, but typically not more than 4. The user must therefore specify the number of correctors in the PISO dictionary by the nCorrectors keyword as shown in the example on page 302.

An additional correction to account for mesh non-orthogonality is available in both SIMPLE and PISO in the standard OpenFOAM solver applications. A mesh is orthogonal if, for each face within it, the face normal is parallel to the vector between the centres of the cells that the face connects, e.g. a mesh of hexahedral cells whose faces are aligned with a Cartesian coordinate system. The number of non-orthogonal correctors is specified by the nNonOrthogonalCorrectors keyword as shown in the examples above and on page 302. The number of non-orthogonal correctors should correspond to the mesh for the case being solved, i.e. 0 for an orthogonal mesh and increasing with the degree of non-orthogonality up to, say, 20 for the most non-orthogonal meshes.

#### 4.5.3.1 Pressure referencing

In a closed incompressible system, pressure is relative: it is the pressure range that matters not the absolute values. In these cases, the solver sets a reference level of pRefValue in cell pRefCell where p is the name of the pressure solution variable. Where the pressure is p_rgh, the names are p_rhgRefValue and p_rhgRefCell respectively. These entries are generally stored in the PISO/SIMPLE sub-dictionary and are used by those solvers that require them when the case demands it. If ommitted, the solver will not run, but give a message to alert the user to the problem.

### 4.5.4 Other parameters

The fvSolutions dictionaries in the majority of standard OpenFOAM solver applications contain no other entries than those described so far in this section. However, in general the fvSolution dictionary may contain any parameters to control the solvers, algorithms, or in fact anything. For a given solver, the user can look at the source code to find the parameters required. Ultimately, if any parameter or sub-dictionary is missing when an solver is run, it will terminate, printing a detailed error message. The user can then add missing parameters accordingly.