OpenFOAM guide/Matrices in OpenFOAM

Linear equations resulting from the discretization method are solved in OpenFOAM using matrix methods. OpenFOAM has a specialized strategy for handling these matrices.

1 Overview

The discretization process results in the equation:

$\mathbf{A} \boldsymbol \Psi = \mathbf B$

Where:

• A is the discretization matrix, which holds the matrix coefficients;
• $\boldsymbol \Psi$ is the vector array of variables being solved for; and
• B is the discretization source term.

The A component is an N x N square matrix, where N is the number of equations produced by the discretization process. OpenFOAM has developed a specialized strategy to handle this matrix. Important components of this strategy include:

• lduMatrix (lower - diagonal - upper - matrix storage class); and
• lduAddressing (an indexing class);
• fvMatrix (contains components specific to the finite volume method).

2 lduMatrix

The lduMatrix class stores the matrix coefficients in three arrays:

scalarField *lowerPtr_;
scalarField *diagPtr_;
scalarField *upperPtr_;

These arrays are Fields that contain all the non-zero entries in the coefficient matrix.

The access of the coefficients goes over the functions:

scalarField& upper()
scalarField& lower()
scalarField& diag()

The size of the arrays lowerPtr_ and upperPtr_ are equal to the number of internal faces. The size of the array diagPtr_ is equal to the number of cells.

2.1 Diagonals

Every cell has one and only one diagonal coefficient, none of which are zero. Therefore the diagonal list contains N values, where N is the number of cells in the mesh. These are denoted:

• As discretization coefficients, $A_p \neq 0\,\!$; or
• As matrix coefficients, $A_{ii} \neq 0\,\!$.

2.2 Upper and lower triangles

The upper and lower triangles have (N - 1)N possible values, however most are zero. The non-zero values correspond to cell-pairs that influence one another. OpenFOAM uses a small computational molecule, therefore only adjacent cells influence one another. In other words, there are two non-zero coefficients for every internal face in the mesh: one that appears in the lower triangle; and one that appears in the upper triangle.

For instance, the next figure shows a simple 3 x 3 mesh:

It consists in 9 cell centers, 12 internal faces and 12 external (boundary) faces. The numbers of the cell centers are shown in black and the number of the faces are shown in blue.

The non-zero coefficients for OpenFOAM would be:

Non-zero coefficients

i	-1-	-2-	-3-	-4-	-5-	-6-	-7-	-8-	-9-
-1-	X	N				N
-2-	O	X	N		N
-3-		O	X	N
-4-			O	X	N				N
-5-		O		O	X	N		N
-6-	O				O	X	N
-7-						O	X	N
-8-					O		O	X	N
-9-				O				O	X

Non-zero coefficients

i	-1-	-2-	-3-	-4-	-5-	-6-	-7-	-8-	-9-
-1-	$A_{1,1}$	$A_{1,2}$				$A_{1,6}$
-2-	$A_{2,1}$	$A_{2,2}$	$A_{2,3}$		$A_{2,5}$
-3-		$A_{3,2}$	$A_{3,3}$	$A_{3,4}$
-4-			$A_{4,3}$	$A_{4,4}$	$A_{4,5}$				$A_{4,9}$
-5-		$A_{5,2}$		$A_{5,4}$	$A_{5,5}$	$A_{5,6}$		$A_{5,8}$
-6-	$A_{6,1}$				$A_{6,5}$	$A_{6,6}$	$A_{6,7}$
-7-						$A_{7,6}$	$A_{7,7}$	$A_{7,8}$
-8-					$A_{8,5}$		$A_{8,7}$	$A_{8,8}$	$A_{8,9}$
-9-				$A_{9,4}$				$A_{9,8}$	$A_{9,9}$

Some mathematical operations require distinguishing between the two coefficients as they are sometimes opposite in sign (for instance when fluxes are involved). OpenFOAM's mesh design anticipated this requirement, which is why every face has an owner cell and a neighbour cell. Owner, neighbour and diagonal coefficients are defined as:

$A_{ij} = \begin{cases} \text{owner}, & i > j \\ \text{diagonal}, & i = j \\ \text{neighbour}, & i < j \end{cases}$

Where:

• i is the matrix row number (corresponds to cell index); and
• j is the matrix column number (corresponds to cell index).

With this definition, the lower triangle has the owner coefficients; and the upper triangle has the neighbour coefficients. They are denoted:

• As discretization coefficients, $A_{r,i}\,\!$, where r indicates "related cells";
• As matrix coefficients:
  • Owner coefficients are: $\left \lbrack A_{ij}\right \rbrack_{i>j} = A_o\,\!$; and
  • Neighbour coefficients are: $\left \lbrack A_{ij}\right \rbrack_{i<j} = A_n\,\!$.

With the above example the entries of the lowerPtr_, upperPtr_ and diagPtr_ array are:

$lower() = A_{2,1}, A_{3,2}, A_{6,1}, A_{5,2}, A_{4,3}, A_{6,5}, A_{5,4}, A_{7,6}, A_{8,5}, A_{9,4}, A_{8,7}, A_{9,8}$

$upper() = A_{1,2}, A_{2,3}, A_{1,6}, A_{2,5}, A_{3,4}, A_{5,6}, A_{4,5}, A_{6,7}, A_{5,8}, A_{4,9}, A_{7,8}, A_{8,9}$

$diag() = A_{1,1}, A_{2,2}, A_{3,3}, A_{4,4}, A_{5,5}, A_{6,6}, A_{7,7}, A_{8,8}, A_{9,9}$

Note that the ordering of the entries of the arrays lowerPtr_ and upperPtr_ are stored according to the face number and of the array diagPtr_ according to the cell number.

With this consistently applied across the mesh, a direction for positive surface-normal vectors is defined: a positive surface-normal vector (such as fluxes) points away from the owner cell, and towards the neighbour cell.

2.3 lduAddressing

The lower, diagonal and upper coefficients are all scalarFields (basically a list), but they are indexed differently. The diagonal coefficients are indexed by cell index; whereas the upper and lower triangles are indexed by face index. In order to get the indices to match, an addressing array is provided for the lower and upper triangles to translate their face index into a corresponding cell index. This addressing array is called an lduAddressing array.

The lduAdressing class (see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1lduAddressing.html) contains two function which return the label list containing the owner (lowerAddr()) and neighbor (upperAddr()) of each face. For our example above the two arrays look like:

lowerAddr() = 1, 2, 1, 2, 3, 5, 4, 6, 5, 4, 7, 8
upperAddr() = 2, 3, 6, 5, 4, 6, 5, 7, 8, 9, 8, 9

The above arrays are ordered according to the face number and contain the owner (lowerAddr()) and neighbor (upperAddr()) of the face.

2.3.1 When do you need it?

lduAddressing is needed anytime you need to perform an operation on all rows of an lduMatrix that involves both diagonal and off-diagonal coefficients.

2.4 How do you use it?

To use lduAddressing:

• separate operations involving the lower triangle from operations involving the upper triangle;
• loop over all the faces in the mesh;
• upper and lower triangle coefficients use the loop index - e.g. upper[facei] = ...
• diagonal coefficients use the lduAddressing coefficient of the face index:
  • diag[l[facei]] -= lower[facei], where l is the lduAddressing for the lower triangle; and
  • diag[u[facei]] -= upper[facei], where u is the lduAddressing for the upper triangle.

2.4.1 Example

For example, the negSumDiag operation is frequently required. This is where the diagonal is the negative sum of its neighbours:

$A_p = -\sum\limits_r A_r$

The guts of its code is:

for (register label face=0; face<l.size(); face++)
{
    Diag[l[face]] -= Lower[face];
    Diag[u[face]] -= Upper[face];
}

Where:

• Diag is the diagonal coefficient;
• l is the lduAddressing for the lower triangle;
• u is the lduAddressing for the upper triangle;
• Lower is the lower triangle; and
• Upper is the upper triangle.

2.5 Diagonal, symmetric, and asymmetric matrices

It is mandatory to define the diagonal; the others are not required. If only one of the lower or upper triangles are defined, the matrix is assumed to be symmetric. If both are defined, it will assume the matrix is asymmetric, even if they are equal.

lowerPtr	diagPtr	upperPtr	Result
Yes/No	No	Yes/No	Invalid matrix - error thrown
No	Yes	No	Diagonal matrix
Yes	Yes	No	Symmetric matrix
No	Yes	Yes	Symmetric matrix
Yes	Yes	Yes	Asymmetric matrix

Caution: Accessing the opposite pointer without a const modifier will convert the matrix to an asymmetric matrix. To modify the off-diagonal of a symmetric matrix, first test which pointer is active using hasUpper() and hasLower().

2.6 OpenFOAM code example

For example, in calculating the diffusion of a quantity $\phi\,\!$, the governing equation is:

$\boldsymbol \nabla \bullet \Gamma \boldsymbol \nabla \mathbf U = 0$

This is the Laplacian, which is calculated in the finite volume method using fvm::laplacian. The discretization coefficients are:

$A_r = - \frac{\Gamma S_f}{\delta}$

$A_p = \sum\limits_r A_r$

To implement this in OpenFOAM:

GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf ( gamma*mesh.magSf() );
fvMatrix<Type>& fvm;
fvm.upper() = deltaCoeffs.internalField() * gamma*mesh.magSf().internalField();
fvm.negSumDiag();

Note:

• fvm.lower is never mentioned - this is because it is equal to fvm.upper and a symmetric matrix is being created.
• fvm.upper is opposite in sign to $A_r$ above - this is because fvm.lower are matrix coefficients, whereas $A_r$ are discretization coefficients.

2.7 Examples how the matrix is build

2.7.1 Laplacian

In the following example we show how the coefficients for the matrix resulting from the discretization of the Laplacian equation are derived. Similar information can be found in the very good book of which describes the details of the discretization in OpenFOAM.

The Laplacian equation reads:

$\boldsymbol \nabla \bullet \Gamma \boldsymbol \nabla \mathbf \phi = 0$

Note the for the explanation here we assume a constant $\Gamma$.

Being $\phi$ the variable solved for and $\Gamma$ the diffusivity.

The above equation is integrated over the control volume and together with the Gauss theorem we get

$\int \boldsymbol \nabla \bullet \Gamma \boldsymbol \nabla \mathbf \phi dV = \int \Gamma \boldsymbol \nabla \phi \mathbf \bullet \boldsymbol n dS = 0$

Substituting the integral over the surface of the volume with a summation over the cell faces we get:

$\sum_f (\Gamma \boldsymbol \nabla \phi )_f\mathbf \bullet \boldsymbol dS = 0$

The the scalar product of the gradient of $\phi$ with the face Area vector $\boldsymbol dS = 0$ is discretized as follows:

$(\Gamma \boldsymbol \nabla \phi )_f\mathbf \bullet \boldsymbol dS = \Gamma \frac{\phi_N - \phi_O}{\delta}S_f$

The face normal vector in OpenFoam is pointing from the face owner (the cell with the lower index) to the face neighbor (the cell with the higher index). This means that if the value $\phi_N$ is higher than $\phi_O$ the scalar product $( \Gamma\boldsymbol \nabla \phi )_f\mathbf \bullet \boldsymbol dS = \Gamma \frac{\phi_N - \phi_O}{\delta}S_f$ is positive and vice versa.

Finally we get:

$\sum_f \Gamma \frac{\phi_{N,f }- \phi_O}{\delta_f}S_f = 0$

Were $\phi_{N,f } , \phi_O, \delta_f, S_f$ are the variable stored at the face neighbor cell center, the variable stored at the owner cell center, the distance from the owner cell center to the neighbor cell center and the face are, respectively. Finally we get for the matrix coefficients for the case where the value of $\phi$ at the boundary is specified:

Non-zero coefficients

i	-1-	-2-	-3-	-4-	-5-	-6-	-7-	-8-	-9-
-1-	$- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_3}{\delta_3}+ \frac{S_{14}}{\delta_{14}} + \frac{S_{25}}{\delta_{25}} )$	$\Gamma\frac{S_1}{\delta_1}$				$\Gamma\frac{S_3}{\delta_3}$
-2-	$\Gamma\frac{S_1}{\delta_1}$	$- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_2}{\delta_2} + \frac{S_4}{\delta_4}+ \frac{S_{15}}{\delta_{15}} )$	$\Gamma\frac{S_2}{\delta_2}$		$\Gamma\frac{S_4}{\delta_4}$
-3-		$\Gamma\frac{S_2}{\delta_2}$	$- \Gamma(\frac{S_2}{\delta_2}+ \frac{S_5}{\delta_5} +\frac{S_{16}}{\delta_{16}} +\frac{S_{17}}{\delta_{17}} )$	$\Gamma\frac{S_5}{\delta_5}$
-4-			$\Gamma\frac{S_5}{\delta_5}$	$- \Gamma(\frac{S_5}{\delta_5}+ \frac{S_7}{\delta_7} + \frac{S_{10}}{\delta_{10}} + \frac{S_{18}}{\delta_{18}} )$	$\Gamma\frac{S_7}{\delta_7}$				$\Gamma\frac{S_{10}}{\delta_{10}}$
-5-		$\Gamma\frac{S_4}{\delta_4}$		$\Gamma\frac{S_7}{\delta_7}$	$- \Gamma(\frac{S_4}{\delta_4}+ \frac{S_6}{\delta_6} + \frac{S_7}{\delta_7} + \frac{S_{9}}{\delta_{9}} )$	$\Gamma\frac{S_6}{\delta_6}$		$\Gamma\frac{S_9}{\delta_9}$
-6-	$\Gamma\frac{S_3}{\delta_3}$				$\Gamma\frac{S_6}{\delta_6}$	$- \Gamma(\frac{S_3}{\delta_3}+ \frac{S_6}{\delta_6} + \frac{S_{8}}{\delta_{8}} +\frac{S_{24}}{\delta_{24}} )$	$\Gamma\frac{S_8}{\delta_8}$
-7-						$\Gamma\frac{S_8}{\delta_8}$	$- \Gamma( \frac{S_8}{\delta_8} + \frac{S_{11}}{\delta_{11}} +\frac{S_{22}}{\delta_{22}} +\frac{S_{23}}{\delta_{23}} )$	$\Gamma\frac{S_{11}}{\delta_{11}}$
-8-					$\Gamma\frac{S_9}{\delta_9}$		$\Gamma\frac{S_{11}}{\delta_{11}}$	$- \Gamma(\frac{S_9}{\delta_9}+ \frac{S_{11}}{\delta_{11}} + \frac{S_{12}}{\delta_{12}} + \frac{S_{21}}{\delta_{21}} )$	$\Gamma\frac{S_{12}}{\delta_{12}}$
-9-				$\Gamma\frac{S_{10}}{\delta_{10}}$				$\Gamma\frac{S_{12}}{\delta_{12}}$	$- \Gamma(\frac{S_{10}}{\delta_{10}} + \frac{S_{12}}{\delta_{12}} + \frac{S_{19}}{\delta_{19}} + \frac{S_{20}}{\delta_{20}} )$

The source vector resulting from the fixed value boundary condition will look like:

source vector for the fixed value boundary conditons

i	-1-	-2-	-3-	-4-	-5-	-6-	-7-	-8-	-9-
-1-	$\Gamma\frac{S_{14}}{\delta_{14}}\phi_{14} + \Gamma\frac{S_{25}}{\delta_{25}}\phi_{25}$	$\Gamma\frac{S_{15}}{\delta_{15}}\phi_{15}$	$\Gamma\frac{S_{16}}{\delta_{16}}\phi_{16} + \Gamma\frac{S_{17}}{\delta_{17}}\phi_{17}$	$\Gamma\frac{S_{18}}{\delta_{18}}\phi_{18}$	$0$	$\Gamma\frac{S_{24}}{\delta_{24}}\phi_{24}$	$\Gamma\frac{S_{22}}{\delta_{22}}\phi_{22} + \Gamma\frac{S_{23}}{\delta_{23}}\phi_{23}$	$\Gamma\frac{S_{21}}{\delta_{21}}\phi_{21}$	$\Gamma\frac{S_{19}}{\delta_{19}}\phi_{19} + \Gamma\frac{S_{20}}{\delta_{20}}\phi_{20}$

For a constant gradient $q$ boundary conditions the matrix looks like:

Non-zero coefficients

i	-1-	-2-	-3-	-4-	-5-	-6-	-7-	-8-	-9-
-1-	$- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_3}{\delta_3} )$	$\Gamma\frac{S_1}{\delta_1}$				$\Gamma\frac{S_3}{\delta_3}$
-2-	$\Gamma\frac{S_1}{\delta_1}$	$- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_2}{\delta_2} + \frac{S_4}{\delta_4} )$	$\Gamma\frac{S_2}{\delta_2}$		$\Gamma\frac{S_4}{\delta_4}$
-3-		$\Gamma\frac{S_2}{\delta_2}$	$- \Gamma(\frac{S_2}{\delta_2}+ \frac{S_5}{\delta_5} )$	$\Gamma\frac{S_5}{\delta_5}$
-4-			$\Gamma\frac{S_5}{\delta_5}$	$- \Gamma(\frac{S_5}{\delta_5}+ \frac{S_7}{\delta_7} + \frac{S_{10}}{\delta_{10}} )$	$\Gamma\frac{S_7}{\delta_7}$				$\Gamma\frac{S_{10}}{\delta_{10}}$
-5-		$\Gamma\frac{S_4}{\delta_4}$		$\Gamma\frac{S_7}{\delta_7}$	$- \Gamma(\frac{S_4}{\delta_4}+ \frac{S_6}{\delta_6} + \frac{S_7}{\delta_7} + \frac{S_{9}}{\delta_{9}} )$	$\Gamma\frac{S_6}{\delta_6}$		$\Gamma\frac{S_9}{\delta_9}$
-6-	$\Gamma\frac{S_3}{\delta_3}$				$\Gamma\frac{S_6}{\delta_6}$	$- \Gamma(\frac{S_3}{\delta_3}+ \frac{S_6}{\delta_6} + \frac{S_{8}}{\delta_{8}} )$	$\Gamma\frac{S_8}{\delta_8}$
-7-						$\Gamma\frac{S_8}{\delta_8}$	$- \Gamma( \frac{S_8}{\delta_8} + \frac{S_{11}}{\delta_{11}} )$	$\Gamma\frac{S_{11}}{\delta_{11}}$
-8-					$\Gamma\frac{S_9}{\delta_9}$		$\Gamma\frac{S_{11}}{\delta_{11}}$	$- \Gamma(\frac{S_9}{\delta_9}+ \frac{S_{11}}{\delta_{11}} + \frac{S_{12}}{\delta_{12}} )$	$\Gamma\frac{S_{12}}{\delta_{12}}$
-9-				$\Gamma\frac{S_{10}}{\delta_{10}}$				$\Gamma\frac{S_{12}}{\delta_{12}}$	$- \Gamma(\frac{S_{10}}{\delta_{10}} + \frac{S_{12}}{\delta_{12}} )$

The source vector resulting from the fixed gradient boundary condition will look like:

source vector for the fixed value boundary conditons

i	-1-	-2-	-3-	-4-	-5-	-6-	-7-	-8-	-9-
-1-	$\Gamma S_{14}q + \Gamma S_{25}q$	$\Gamma S_{15}q$	$\Gamma S_{16}q + \Gamma S_{17}q$	$\Gamma S_{18}q$	$0$	$\Gamma S_{24} q$	$\Gamma S_{22}q + \Gamma S_{23}q$	$\Gamma S_{21}q$	$\Gamma S_{19}q + \Gamma S_{20}q$

 
template<class Type, class GType>
tmp<fvMatrix<Type>>
gaussLaplacianScheme<Type, GType>::fvmLaplacian
(
    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = this->mesh();
 
    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
 
    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
    (
        SfGamma & Sn
    );
    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
 
    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected
    (
        SfGammaSn,
        this->tsnGradScheme_().deltaCoeffs(vf),
        vf
    );
    fvMatrix<Type>& fvm = tfvm.ref();
 
    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tfaceFluxCorrection
        = gammaSnGradCorr(SfGammaCorr, vf);
 
    if (this->tsnGradScheme_().corrected())
    {
        tfaceFluxCorrection.ref() +=
            SfGammaSn*this->tsnGradScheme_().correction(vf);
    }
 
    fvm.source() -= mesh.V()*fvc::div(tfaceFluxCorrection())().primitiveField();
 
    if (mesh.fluxRequired(vf.name()))
    {
        fvm.faceFluxCorrectionPtr() = tfaceFluxCorrection.ptr();
    }
 
    return tfvm;
}
 
 
 
template<class Type, class GType>
tmp<fvMatrix<Type>>
gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
(
    const surfaceScalarField& gammaMagSf,
    const surfaceScalarField& deltaCoeffs,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();
 
    fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
    fvm.negSumDiag();
 
    forAll(vf.boundaryField(), patchi)
    {
        const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
        const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
        const fvsPatchScalarField& pDeltaCoeffs =
            deltaCoeffs.boundaryField()[patchi];
 
        if (pvf.coupled())
        {
            fvm.internalCoeffs()[patchi] =
                pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
            fvm.boundaryCoeffs()[patchi] =
               -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
        }
        else
        {
            fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
            fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
        }
    }
 
    return tfvm;
}
 

 
 
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedValueFvPatchField<Type>::gradientInternalCoeffs() const
{
    return -pTraits<Type>::one*this->patch().deltaCoeffs();
}
 
 
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedValueFvPatchField<Type>::gradientBoundaryCoeffs() const
{
    return this->patch().deltaCoeffs()*(*this);
}
 

 
 
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedGradientFvPatchField<Type>::gradientInternalCoeffs() const
{
    return tmp<Field<Type>>::New(this->size(), Zero);
}
 
 
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedGradientFvPatchField<Type>::gradientBoundaryCoeffs() const
{
    return gradient();
}
 

3 fvMatrix

Since OpenFOAM implements multiple numerical strategies such as finite difference, finite area, finite element, and so on, everything that is specific to any one of these strategies is stripped out of the base class lduMatrix. This facilitates code reuse.

fvMatrix is the finite volume extension of lduMatrix. It includes the source, B, as well as a reference to the field, $\boldsymbol \Psi$ from $\mathbf A \boldsymbol \Psi = \mathbf B$. It handles what happens on the boundaries, and also implements specific operations that occur frequently in finite volume algorithms, such as fvMatrix;;H.