A little discuss about fvMatrix fvMesh lduMatrix lduMesh in OpenFOAM

A little discuss about fvMatrix fvMesh lduMatrix lduMesh in OpenFOAMfvMatrix inherit lduMatrix, therefore it has the function of lduMatrix. To have a deeper look at the linear system (E.g. AX = b) we are going to solve, we need to know the matrix A of the linear system. E.g. if we have:

1
2
3
4
5
6
fvScalarMatrix TEqn
(
    fvm::ddt(T) - fvm::laplacian(DT, T)
 ==
    fvOptions(T)
);

corresponding to equation:

The fvm in OpenFOAM means implicit disretization. Therefore, a linear system of has to be solved. To dig into the information of this matrix A, OpenFOAM gives us lots of choices.

In fvMatrix.H, there are several access functions called:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
// Member Functions
 
    // Access
 
        const GeometricField<Type, fvPatchField, volMesh>& psi() const
        {
            return psi_;
        }
 
        const dimensionSet& dimensions() const
        {
            return dimensions_;
        }
 
        Field<Type>& source()
        {
            return source_;
        }
 
        const Field<Type>& source() const
        {
            return source_;
        }
 
        //- fvBoundary scalar field containing pseudo-matrix coeffs
        //  for internal cells
        FieldField<Field, Type>& internalCoeffs()
        {
            return internalCoeffs_;
        }
 
        //- fvBoundary scalar field containing pseudo-matrix coeffs
        //  for boundary cells
        FieldField<Field, Type>& boundaryCoeffs()
        {
            return boundaryCoeffs_;
        }
 
 
        //- Declare return type of the faceFluxCorrectionPtr() function
        typedef GeometricField<Type, fvsPatchField, surfaceMesh>
            *surfaceTypeFieldPtr;
 
        //- Return pointer to face-flux non-orthogonal correction field
        surfaceTypeFieldPtr& faceFluxCorrectionPtr()
        {
            return faceFluxCorrectionPtr_;
        }

In lduMatrix, there are several access functions called:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
// Member Functions
 
    // Access to addressing
 
        //- Return the LDU mesh from which the addressing is obtained
        const lduMesh& mesh() const
        {
            return lduMesh_;
        }
 
        //- Return the LDU addressing
        const lduAddressing& lduAddr() const
        {
            return lduMesh_.lduAddr();
        }
 
        //- Return the patch evaluation schedule
        const lduSchedule& patchSchedule() const
        {
            return lduAddr().patchSchedule();
        }
 
        //- Return interfaces
        const LduInterfaceFieldPtrsList<Type>& interfaces() const
        {
            return interfaces_;
        }
 
        //- Return interfaces
        LduInterfaceFieldPtrsList<Type>& interfaces()
        {
            return interfaces_;
        }
 
 
    // Access to coefficients
 
        Field<DType>& diag();
        Field<LUType>& upper();
        Field<LUType>& lower();
        Field<Type>& source();
 
        FieldField<Field, LUType>& interfacesUpper()
        {
            return interfacesUpper_;
        }
 
        FieldField<Field, LUType>& interfacesLower()
        {
            return interfacesLower_;
        }
 
 
        const Field<DType>& diag() const;
        const Field<LUType>& upper() const;
        const Field<LUType>& lower() const;
        const Field<Type>& source() const;
 
        const FieldField<Field, LUType>& interfacesUpper() const
        {
            return interfacesUpper_;
        }
 
        const FieldField<Field, LUType>& interfacesLower() const
        {
            return interfacesLower_;
        }

The lduMesh_ in lduMatrix is defined as:

1
2
//- LDU mesh reference
const lduMesh& lduMesh_;

So, when you call lduAddr() in lduMatrix, it is actually calls the lduAddr() function in lduMesh class (because of the return of the function).

has a type of ```fvScalarMatrix```, which is a specified template class of ```fvMatrix```. This has been discussed before
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

In the constructor of ```fvMatrix```:

```cpp
 template<class Type>
 Foam::fvMatrix<Type>::fvMatrix
 (
     const GeometricField<Type, fvPatchField, volMesh>& psi,
     const dimensionSet& ds
 )
 :
     lduMatrix(psi.mesh()),
     psi_(psi),
     dimensions_(ds),
     source_(psi.size(), Zero),
     internalCoeffs_(psi.mesh().boundary().size()),
     boundaryCoeffs_(psi.mesh().boundary().size()),
     faceFluxCorrectionPtr_(nullptr)
is constructed by ```psi.mesh()```, which actually has a type of ```fvMesh``` in this case. The constructor of ```lduMatrix``` is shown here:
1
2
3
4
5
6
7
8
9

```cpp
 Foam::lduMatrix::lduMatrix(const lduMesh& mesh)
 :
     lduMesh_(mesh),
     lowerPtr_(nullptr),
     diagPtr_(nullptr),
     upperPtr_(nullptr)
 {}

The lduMesh_ which is defined as:

1
2
//- LDU mesh reference
const lduMesh& lduMesh_;

will be initialized with psi.mesh() which has a type of fvMesh. There is no doubt this is logical cause lduMesh is the base class of fvMesh. Therefore we we call lduAddr() in lduMatrix.H :

1
2
3
4
5
//- Return the LDU addressing
const lduAddressing& lduAddr() const
{
    return lduMesh_.lduAddr();
}

It is actually calling the lduAddr() function in class fvMesh:

1
2
3
4
5
6
7
8
9
const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
{
    if (!lduPtr_)
    {
        lduPtr_ = new fvMeshLduAddressing(*this);
    }

    return *lduPtr_;
}

The interesting thing is, in all the constructor of fvMesh, lduPtr_ are all defined as a nullptr. Therefore, the if loop has to run first. After the first run of the if loop, it will not be called again cause lduPtr_ is not a nullptr anymore. The *this here represents the fvMesh class objects. The lduPtr_ is defined in fvMesh.H as:

1
mutable fvMeshLduAddressing* lduPtr_;

In fvMeshLduAddressing.H, there will be LDU addressing informations:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// Member Functions
 
    //- Return lower addressing (i.e. lower label = upper triangle)
    const labelUList& lowerAddr() const
    {
        return lowerAddr_;
    }
 
    //- Return upper addressing (i.e. upper label)
    const labelUList& upperAddr() const
    {
        return upperAddr_;
    }
 
    //- Return patch addressing
    const labelUList& patchAddr(const label i) const
    {
        return *patchAddr_[i];
    }
 
    // Return patch field evaluation schedule
    const lduSchedule& patchSchedule() const
    {
        return patchSchedule_;
    }

To conclude, TEqn which has a type of fvScalarMatrix (actually fvMatrix<Scalar>) will be used here. There are several cases can happen:

  1. If we call TEqn.dimensions(), function dimensions() in class fvMatrix will be called. This is the same for function psi(), source(), internalCoeffs(), boundaryCoeffs() and faceFluxCorrectionPtr().

  2. If we call TEqn.mesh(), function mesh() in class lduMatrix will be called. This is the same for lower(), diag() and upper(). The last three functions are always used in peeping the coefficients in matrix of the linear equation.

  3. If we call TEqn.lduAddr(), function lduAddr() in lduMatrix will be called first. Then, we will have a tunnel to check some mesh LDU addressing. As what has been mentioned, lduMesh_ will be initilized with a type of fvMesh. So, it will call the lduAddr() function in fvMesh class, which will return a objects of type of fvMeshLduAddressing. Then we can further call the access function lowerAddr(), upperAddr(), patchAddr(const label i) and patchSchedule() in fvMeshLduAddressing class. It should be noted that fvMeshLduAddressing is a subclass of lduAddressing. Therefore, in function lduAddr(), the return type is lduAddressing&.