This entry is part of 5 in the series Run-Time Type Selection

Featured Image: Math Wall by João Trindade

A question popped up during the Programmer training on calculating gradients using alternative schemes to those automatically chosen by the gradient operator from the system/fvSchemes file. In OpenFOAM, the operators rely on the RTS mechanism to select a specific differencing scheme. A subset of discrete operators support named scheme selection that allows the user of the operator to name the scheme to be used by the operator. This post contains a small example that shows how to select a scheme using the name argument of the operator, as well as a description of the scheme selection implementation for the discrete (Finite Volume) operators in OpenFOAM.

Scheme selection example

  1. Add the scheme name parameter to system/fvSchemes.
  2. Call the operator with two arguments: a field and a scheme name.

 

The following example shows how to do that for a gradient calculation. The system/fvSchemes file is modified like this

    gradSchemes
    {
        default         Gauss linear;
        alphaGrad       pointCellsLeastSquares;
    }

This small application loads the alpha.water field, computes two gradient fields and saves them:

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
    #include "fvCFD.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    // Main program:

    int main(int argc, char *argv[])
    {
        #include "setRootCase.H"
        #include "createTime.H"
        #include "createMesh.H"

        volScalarField alpha 
        (
            IOobject
            (
                "alpha.water", 
                runTime.timeName(), 
                mesh, 
                IOobject::MUST_READ, 
                IOobject::AUTO_WRITE
            ),
            mesh
        );

        tmp<volVectorField> alphaGradTmp = fvc::grad(alpha, "alphaGrad"); 
        const volVectorField& alphaGrad = alphaGradTmp(); 

        tmp<volVectorField> linearGradTmp = fvc::grad(alpha); 
        volVectorField& linearGrad = linearGradTmp(); 
        linearGrad.rename("linearGrad"); 

        linearGrad.write(); 
        alphaGrad.write(); 

        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

        Info<< "\nEnd\n" << endl;
        return 0;
    }

The application is available in the sourceflux blog code repository, in applications/postProcessing/calcGrad. You can compile it and run it in any simulation case directory that contains a scalar field named alpha.water in the 0 directory. I have used the tutorials/multiphase/laminar/rising-bubble-2D in the sourceflux blog code repository. If you want to understand why the line linearGrad.rename("linearGrad"); is used to save the gradient computed with the default scheme (Gauss linear), read the following sections that cover the implementation of the scheme selection process for operators in OF.

Note
The gradient calculation usually involves steps: operator approximation and field value interpolation. The default gradient uses Gauss approximation of a cell-centered gradient (second order accuracy), and an interpolation scheme that is required for the face-centered values. We cover the FVM discretization as an optional topic in our trainings.

Scheme selection in transport equations

Scheme selection for an operator that is a part of a transport equation is done using a prescribed naming rule. Wherever the operator is called in the source code

   tmp<volVectorField> gradAlpha = fvc::grad(alpha); 

The operator only forwards the work to the scheme class (template), after selecting the scheme based on the entry in system/fvSchemes:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
    template<class Type>
    tmp
    <
        GeometricField
        <
            typename outerProduct<vector,Type>::type, fvPatchField, volMesh
        >
    >
    grad
    (
        const GeometricField<Type, fvPatchField, volMesh>& vf
    )
    {
        return fvc::grad(vf, "grad(" + vf.name() + ')');
    }

From the return call site, we can see that the name used for scheme selection will be grad(nameOfField), where nameOfField is the name given to the field constructor (IOobject) argument somewhere in the application. For solver applications, the fields are initialized in the createFields.H file, so they are named there as well. If you are programming a new application, you can choose where to initialize the fields and set their names. You can fix the names, let the application user define the names in the command line, or make the application read the field name from a configuration file. The operator above forwards the call to the corresponding operator that works with two arguments (a field and a word):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
    tmp
    <
        GeometricField
        <
            typename outerProduct<vector,Type>::type, fvPatchField, volMesh
        >
    >
    grad
    (
        const GeometricField<Type, fvPatchField, volMesh>& vf,
        const word& name
    )
    {
        return fv::gradScheme<Type>::New
        (
            vf.mesh(),
            vf.mesh().gradScheme(name)
        )().grad(vf, name);
    }

Here we can see that the gradient operator makes the scheme class (template) do all the work. Kind of. Let’s go over the important lines to make sure we get everything right.

1
2
3
4
5
    return fv::gradScheme<Type>::New
    (
        vf.mesh(),
        vf.mesh().gradScheme(name)
    )().grad(vf, name);

The return line returns the type declared by the operator. This will be a tmp<GeometricField<typename outerProduct... Long story short – a gradient of a scalar field is a vector field, so the gradient in the fvc::grad(alpha) call will return a tmp<volVectorField>. Short story long : visit our Programmer Training and learn stuff like traits, policies, template specialization, smart pointers, and other C++ template stuff important for efficient programming in OF.

The return forwards the call to fv::gradScheme<Type>::New – a static class function from the gradScheme class template used for RTS – Runtime Type Selection. A fancy way of saying that the New function creates an object of a specific scheme that we want to be used in our gradient calculation. If you want to understand the workings of RTS, follow the link to our blog series on RTS above.

Now we come to a part that’s kinda confusing to most people: New(some, arguments)().grad(some, arguments). The question that we usually hear on the training is what the () brackets are doing. Well, the New operator returns an OpenFOAM smart pointer (which kind depends on the RTS table declaration/definition), and the () operator is dereferencing this pointer. It becomes clear when written like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    // Select the scheme: build a scheme object and wrap it in a smart
    // pointer.
    tmp<gradScheme<Type> > selectedSchemePtr = fv::gradScheme<Type>::New
    (
        vf.mesh(),
        vf.mesh().gradScheme(name)
    );

    // Get the reference to a scheme: dereference the smart pointer with
    // tmp<gradScheme<Type> >::operator()() call operator. 
    gradScheme<Type>& selectedScheme = selectedSchemePtr();  
    
    // Make the scheme do the work and return the tmp<GeometricField<.. 
    return selectedScheme.grad(vf, name);

We know which smart pointer is returned by the scheme New selector because it is declared in the gradScheme.H file:

1
2
3
4
5
6
    //- Return a pointer to a new gradScheme created on freestore
    static tmp<gradScheme<Type> > New
    (
        const fvMesh& mesh,
        Istream& schemeData
    );

That is all nice and now we understand how the scheme selection is fixed with respect to the naming structure. Well, there is one question left: how do we know that the file system/fvSchemes is used for scheme selection?

Named scheme selection

The question on selecting schemes can be extended from the gradient to selecting an alternative scheme for any other operator in the OF toolbox, because they all work in exactly the same way when it comes to scheme selection. In this post, we stick to computing gradients. However, exactly same thing can be done for most discrete (differential) operators in [fFoOaAmM]. In order to allow an alternative, named scheme selection, the operator has to provide an implementation that takes both the field and the name as arguments. The grad operator provides this implementation, it is posted in the code snippet above. An operator with this signature will do two things

  1. Select the scheme.
  2. Call the scheme.
  3. Return the result from the scheme.

 

Even though some operators like the gradient do 1-3 in the same line, semantically, these 3 operations are distinct. If you haven’t figured out so far how the scheme lookup is done from system/fvSchemes and you don’t have time to read the blog posts on RTS, the answer lies in the following line from the above gradient operator: vf.mesh().gradScheme(name). Each geometric field holds a handle (in this case it is a const reference) to the mesh by inheritance from a dimensioned field:

1
2
3
4
5
6
   template<class Type, class GeoMesh>
   inline const typename GeoMesh::Mesh&
   Foam::DimensionedField<Type, GeoMesh>::mesh() const
   {
       return mesh_;
   }

That accounts for the vf.mesh() part. Now comes the interesting part: The mesh is the fvSchemes dictionary:

1
2
3
4
5
6
7
8
9
10
    class fvMesh
    :
        public polyMesh,
        public lduMesh,
        public surfaceInterpolation,
        public fvSchemes,
        public fvSolution,
        public data
    {
        // Private data

The finite volume fv operators work on surface* and vol* (dimensioned) geometrical fields used by the finite volume method, and fvMesh implements the mesh that supports finite volume numerics. From this point, everything becomes straightforward. The fvSchemes class delivers data from system/fvSchemes directory, precisely – the name of the scheme used by the operator:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
    Foam::ITstream& Foam::fvSchemes::gradScheme(const word& name) const
    {
        if (debug)
        {
            Info<< "Lookup gradScheme for " << name << endl;
        }

        if (gradSchemes_.found(name) || defaultGradScheme_.empty())
        {
            return gradSchemes_.lookup(name);
        }
        else
        {
            const_cast<ITstream&>(defaultGradScheme_).rewind();
            return const_cast<ITstream&>(defaultGradScheme_);
        }
    }

The gradSchemes_ data member is as expected a dictionary. This means that the location of the scheme name parameter is fixed, it is prescribed by the fv operators, since they forward the search request for the scheme name to the mesh (fvSchemes). As a consequence, in the first section, we can write alphaGrad pointCellsLeastSquares in the system/fvSchemes dictionary. The fvSchemes class will look up the gradScheme for the name (alphaGrad), and pass the scheme type information to the scheme New selector. In this tutorial, it was the gradScheme<Type>::New selector.

Conclusions

Nothing much to conclude here. The scheme selection is prescribed by the operators in fv. Mesh is used to read the scheme parameters from a prescribed place (system/fvOptions). You can add your name to select a scheme in any line of your code by using a syntax like operator(field, name). The table below lists some of the operators with and without allowing the operator caller the ability to select schemes using names other than operator(fieldName)' insystem/fvSchemes`. The table may not be complete, if you want to contribute and have found another operator that can be added to the table, leave a comment below and I’ll check the operator and extend the table.

OperatorNamed scheme selection
average[o] Fixed linear interpolation.
curl[x] Used for the gradient.
div[x]
flux[x]
grad[x]
laplacian[x]
magSqrGradGrad[o]
snGrad[x]
Series Navigation

Run-time type selection in OpenFOAM – the type name system >>