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.
Contents
1 Overview
The discretization process results in the equation:
Where:
- A is the discretization matrix, which holds the matrix coefficients;
- 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:
- As discretization coefficients, ; or
- As matrix coefficients, .
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:
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:
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, , where r indicates "related cells"^{[1]};
- As matrix coefficients:
- Owner coefficients are: ; and
- Neighbour coefficients are: .
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 , the governing equation is:
This is the Laplacian, which is calculated in the finite volume method using fvm::laplacian. The discretization coefficients are:
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 above - this is because fvm.lower are matrix coefficients, whereas 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:
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.
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, from . It handles what happens on the boundaries, and also implements specific operations that occur frequently in finite volume algorithms, such as fvMatrix;;H.