OpenFOAM guide/Matrices in OpenFOAM

From OpenFOAMWiki

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


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

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:

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, given a simple 3 x 3 mesh:

-1- -2- -3-
-6- -5- -4-
-7- -8- -9-

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

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} =
\text{owner}, & i > j \\
\text{diagonal}, & i = j \\
\text{neighbour}, & i < j


  • 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:

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

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.

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.

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

3.3 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];


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

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

5 Notes

  1. What I'm calling related cells is conventionally called neighbours. But OpenFOAM has a different meaning for neighbours, so the term related cells is used.