fvm::laplacian in OpenFOAM

fvm::laplacian in OpenFOAM

Start from the deeper look into the laplacianFoam solver. We are going to solve the:

1
2
3
4
5
6
7
8
9
10
11
12
13
while (simple.correctNonOrthogonal())
{
    fvScalarMatrix TEqn
    (
        fvm::ddt(T) - fvm::laplacian(DT, T)
     ==
        fvOptions(T)
    );

    fvOptions.constrain(TEqn);
    TEqn.solve();
    fvOptions.correct(T);
}

In TEqn, we have a object called fvm::laplacian(DT, T) in its constructor. We will have a deep look into this term. For fvm::laplacian(DT, T), function laplacian in fvmLaplacian.C will be called:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
template<class Type, class GType>
tmp<fvMatrix<Type>>
laplacian
(
    const dimensioned<GType>& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const GeometricField<GType, fvsPatchField, surfaceMesh> Gamma
    (
        IOobject
        (
            gamma.name(),
            vf.instance(),
            vf.mesh(),
            IOobject::NO_READ
        ),
        vf.mesh(),
        gamma
    );

    return fvm::laplacian(Gamma, vf);
}

Here the Type and GType will be defined by the type of gamma and vf. Here both of them will be scalar. The return object will call another overloaded laplacian function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
template<class Type, class GType>
tmp<fvMatrix<Type>>
laplacian
(
    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return fvm::laplacian
    (
        gamma,
        vf,
        "laplacian(" + gamma.name() + ',' + vf.name() + ')'
    );
}

Here with a name is defined as "laplacian(" + gamma.name() + ',' + vf.name() + ')', and then another overloaded laplacian function will be called:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
template<class Type, class GType>
tmp<fvMatrix<Type>>
laplacian
(
    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const word& name
)
{
    return fv::laplacianScheme<Type, GType>::New
    (
        vf.mesh(),
        vf.mesh().laplacianScheme(name)
    ).ref().fvmLaplacian(gamma, vf);
}

Here, the fv::laplacianScheme<Type, GType>::New will be called, which is defined in laplacianScheme.C:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 template<class Type, class GType>
 tmp<laplacianScheme<Type, GType>> laplacianScheme<Type, GType>::New
 (
     const fvMesh& mesh,
     Istream& schemeData
 )
 {
     if (fv::debug)
     {
         InfoInFunction << "Constructing laplacianScheme<Type, GType>" << endl;
     }
 
...
 
     const word schemeName(schemeData);
 
     typename IstreamConstructorTable::iterator cstrIter =
         IstreamConstructorTablePtr_->find(schemeName);
 
...
 
     return cstrIter()(mesh, schemeData);
 }

In laplacianScheme.H :

1
2
3
4
5
6
7
8
declareRunTimeSelectionTable
(
    tmp,
    laplacianScheme,
    Istream,
    (const fvMesh& mesh, Istream& schemeData),
    (mesh, schemeData)
);

The RunTimeSelectionTable is declared, therefore, IstreamConstructorTablePtr_ can be used. Cause we use Gauss scheme in the fvScheme file, therefor here the returned value from this New function will be gaussLaplacianScheme.

The ref() is a function used to dereference the tmp pointer, which will return the dereferenced pointer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
template<class T>
inline T& Foam::tmp<T>::ref() const
{
    if (isTmp())
    {
        if (!ptr_)
        {
            FatalErrorInFunction
                << typeName() << " deallocated"
                << abort(FatalError);
        }
    }
    else
    {
        FatalErrorInFunction
            << "Attempt to acquire non-const reference to const object"
            << " from a " << typeName()
            << abort(FatalError);
    }

    return *ptr_;
}

to be continued…