This entry is part 2 of 2 in the series Boundary conditions

Information on the coded fixed value BC is available in the official release notes, in various forum threads like here and here, etc. There still seems to be some confusion about the boundary field assignment operator==(const fvPatchField&) as well as the availability of other (registered) objects that can be used within the coded fixed value BC. This blog post is an attempt to make this additional (C++) stuff clear.

A codedFixedValue tutorial case

We start with a simulation case for the interFoam solver, where we want to use the coded fixed value BC to define an inlet for the velocity field U and the volume fraction field α.

codedFixedValueSpecs

The test case is quite simple, it is a box that gets filled with water via a square shaped inlet inlet defined on one of its sides.

Note
The BCs described in this tutorial are available in the sourceflux blog code repository, in the tutorials/multiphase/interFoam/laminar/roomInflow tutorial case.

A codedFixedValue square inlet

Basic information on how to work with this BC is available in the official C++ Doxygen documentation. In this example, an inlet is generated for both the volume fraction α and the velocity U field and non-uniform values are set. The code used for both fields is more/less exactly the same, so only the velocity field is described in this section. You can compare this to the implementation of the coded fixed value BC for the α field available in the sourceflux repository.

Here is a complete definition of the coded fixed value BC that sets the inflow velocity to (0.1, 0, 0) for a sub-set of boundary patch faces whose centers are within a square defined by points P1 = (0, 0.501, 0.501) and P2 = (0, 0.751, 0.751):

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
    inlet
    {
        type            codedFixedValue;
        value           uniform (0 0 0);
        redirectType    velocitySquareInlet;   

        code
        #{
            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 

            vectorField& field = *this; 
            scalar min = 0.501;  
            scalar max = 0.751; 

            forAll(Cf, faceI)
            {
               if (
                    (Cf[faceI].z() > min) &&
                    (Cf[faceI].z() < max) &&
                    (Cf[faceI].y() > min) &&
                    (Cf[faceI].y() < max) 
                  )
                {
                    field[faceI] = vector(0.1, 0, 0);
                }
            }
        #};
    }

Let’s go over the sets of input parameters. The first set is clear, apart from probably redirectType.

parametermeaning
typeUsed by the RTS mechanism to select the BC.
valueDefault value.
redirectTypeName of the generated BC class.

The code snippet that documents the codedFixedValue BC can also be found in the official documentation. The redirectType is a name of the generated BC. That might be difficult do decipher. Basically, the coded fixed value BC does what all other coded objects in OpenFOAM do:

  1. Create the dynamicCode directory in the simulation case directory.
  2. Copy the template files of the appropriate coded object into the dynamicCode.
    • Create a sub-directory named as the value of the redirectType parameter.
    • All listed coded objects are compiled into separate dynamic libraries.
    • A single class is generated from template files per coded object.
    • The template files also include the required build files.
  3. Modify the source code of the template using user input.
  4. Build the library.

Depending on which coded object type is used, different member functions will be altered by the code provided by the user. Be it coded source finite volume option, coded function object or a coded boundary condition, they all have different member functions. Coded objects will be composited by a class. In fact, they are composited in a one-to-many relationship, because the user can create many coded objects at runtime. For example, there is no limit in the number of coded function objects that can be defined in the system/controlDict file. For reasons that we cover in depth in our programmer training, coded function objects are composited in a list by the Foam::Time class, so many of them can be appended to this list.

Note
The generated libraries in dynamicCode are linked to the solver at run-time, not compile time. This makes sense, since the libraries are case-specific and generated automagically by a meta-programming system that works with template files.

Depending on the type of the coded object and its member function that we are coding with the coded object, different variables and functions will be available to us. Since in this example a coded fixed value boundary condition is used, it is a boundary condition, and in OpenFOAM, this means that it is an fvPatchField. So, we can use all public and protected member functions and data members of fvPatchField. This brings us to the code part of our coded fixed value inlet BC:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
    code
    #{
        const fvPatch& boundaryPatch = patch(); 
        const vectorField& Cf = boundaryPatch.Cf(); 

        vectorField& field = *this; 

        scalar min = 0.501;  
        scalar max = 0.751; 

        forAll(Cf, faceI)
        {
           if (
                (Cf[faceI].z() > min) &&
                (Cf[faceI].z() < max) &&
                (Cf[faceI].y() > min) &&
                (Cf[faceI].y() < max) 
              )
            {
                field[faceI] = vector(1, 0, 0);
            }
        }
    #};

Take a look at the official C++ documentation of fvPatchField. It has a pretty fat interface, and we can use all of it. In order to generate a rectangle inlet for the velocity U, we need to find a way to access the geometry of the patch. We can work with face points or face centers. Since the boundary field of a vol*Field maps to face centers (covered in the Primer and the Programmers Guide), face centers are chosen to create the rectangle. Getting access to the face centers of the patch is done like this

1
2
3
4
    code
    #{
        const fvPatch& boundaryPatch = patch(); 
        const vectorField& Cf = boundaryPatch.Cf(); 

The fvPatchField is a field but it has a patch, namely a constant reference to the boundary mesh patch:

1
2
3
4
5
6
7
8
9
    template<class Type>
    class fvPatchField
    :
        public Field<Type>
    {
        // Private data

            //- Reference to patch
            const fvPatch& patch_;

It offers a public member function fvPatchField<Type>::patch() that retrieves the necessary constant reference. Since codedFixedValueFvPatchField publicly inherits from fvPatchField<Type>, we can call patch() without problems.

We want to be able to set the individual face centered values differently depending on some spatial and/or temporal conditions. This requires a non-constant access to the field stored in our coded boundary condition. Well, as noted before, the coded boundary condition is a fvPatchField and fvPatchField is a Field<Type>. Since we are working with the volVectorField U, we can upcast the coded boundary condition to its grandparent class Field<vector> a.k.a. vectorField, and work with that.

Note
Upcasting is safe in C++, because it is a basic object orientation principle, called Liskov subsitution principle. The Field<Type> defines a sub-set of the codedFixedValueFvPatchField type, which means that any object instantiated from codedFixedValueFvPatchField will also behave like an object of Field<Type>.

The rest should be self-explanatory. To define the rectangle, I’ve used min and max for the z and y coordinates and set the inlet velocity to (1, 0, 0) for all faces whose centers are within the rectangle:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
    scalar min = 0.501;  
    scalar max = 0.751; 

    forAll(Cf, faceI)
    {
       if (
            (Cf[faceI].z() > min) &&
            (Cf[faceI].z() < max) &&
            (Cf[faceI].y() > min) &&
            (Cf[faceI].y() < max) 
          )
        {
            field[faceI] = vector(1, 0, 0);
        }
    }

For the volume fraction field, the inlet rectangle is positioned in exactly the same way as for the velocity field:

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
    inlet
    {
        type            codedFixedValue;
        value           uniform 0;
        redirectType    alphaSquareInlet;   

        code
        #{
            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 

            scalarField& field = *this; 
            field = patchInternalField(); 

            scalar min = 0.501;  
            scalar max = 0.751; 

            forAll(Cf, faceI)
            {
               if (
                    (Cf[faceI].z() > min) &&
                    (Cf[faceI].z() < max) &&
                    (Cf[faceI].y() > min) &&
                    (Cf[faceI].y() < max) 
                  )
                {
                    field[faceI] = 1; 
                }
            }
        #};
    
    }

The main difference is in the type of the tensor we are using: for the volume fraction we are dealing with scalars. You can also examine the definition of the codedFixedValue BC applied to the volume fraction field in the sourceflux blog code repo. The official documentation makes use of the fvPatchField assignment operator

1
2
3
4
    code
    #{
        operator==(min(10,0.1,this->db().time().value()));
    #};

This line basically tells the compiler to look up the name of the function operator== and call the appropriate function with the argument resulting from min. The process of looking up a name of the function and/or variable is called name lookup and detail information on that topic can be found here and here. In this example, the compiler finds a member function operator of the parent class fixedValueFvPatchField<vector>, and calls it to set the uniform boundary field value resulting from the min function.

Note
To work with coded objects in OpenFOAM, it is important to understanding two things: which class is being implemented by the coded object and which member function is metaprogrammed by which code block. The names of the functions and objects used in the code block will be looked up by the compiler following a set of name lookup rules.

The codedFixedValue dynamic library

Once the codedFixedValue BC has been applied to both the velocity and the member function, a switch needs to be turned on in order to allow OpenFOAM to build dynamic code. The switch is located in the $WM_PROJECT_DIR/etc/controlDict file:

1
2
3
4
5
6
7
8
9
10
    InfoSwitches
    {
        writePrecision  6;
        writeJobInfo    0;
        writeDictionaries 0;
        writeOptionalEntries 0;

        // Allow case-supplied C++ code (#codeStream, codedFixedValue)
        allowSystemOperations   1;
    }

Executing the solver causes the dynamic code system to build the appropriate libraries following the steps covered in the previous section. For example, for the velocity field, in the dynamicCode directory, we find the velocitySquareInlet library directory:

1
2
3
4
    ?> ls dynamicCode
    ?> ls dynamicCode/velocitySquareInlet/
    fixedValueFvPatchFieldTemplate.C  fixedValueFvPatchFieldTemplate.dep  
    fixedValueFvPatchFieldTemplate.H  lnInclude  Make

This directory contains the template files with the automatically generated fixedalueFvPatchField<vector> for our velocity field. The library naming is using SHA1 hashes to for library names:

1
2
3
4
5
6
    /* dynamicCode:
     * SHA1 = 91231d9d7ebcda548c39dc7d6ed624a2b3222e33
     */
    fixedValueFvPatchFieldTemplate.C

    LIB = $(PWD)/../platforms/$(WM_OPTIONS)/lib/libvelocitySquareInlet_91231d9d7ebcda548c39dc7d6ed624a2b3222e33

This is done because for a simulation case, there is no limit to the number of dynamic libraries that can be defined by the user. For the boundary conditions, the upper limit is $N_{patch} \cdot \N_{fields}$, since we could hypothetically applied a coded BC for each patch of each field. Those of you that work with technical applications that involve a large number of outlets/inlets know how large this number can get. That’s why using hash system to generate unique names makes perfect sense.

If we examine the declaration (.H) file in the velocitySquareInlet directory, we’ll see that the dynamic code mechanism meta-programs concrete types:

1
2
3
4
5
6
7
8
9
10
11
    class velocitySquareInletFixedValueFvPatchVectorField
    :
        public fixedValueFvPatchField<vector>
    {
    public:

        //- Information about the SHA1 of the code itself
        static const char* const SHA1sum;

        //- Runtime type information
        TypeName("velocitySquareInlet");

A class template is not generated, but a concrete class that inherits from specialized fixedValueFvPatchField<vector>. This has a serious consequence for the code we program into the code block, because the name lookup mentioned in the previous section will allways be independent : it will not depend on the template parameters, since there are none. Using fixedVlaueFvPatchField<vector> causes the compiler to explicitly instantiate the template with the vector argument. That is why we don’t have to qualify the name with the this operator when calling the patch member function

    const fvPatch& boundaryPatch = this->patch(); 

These minor details become important when one starts working with the dynamic code system, and the time comes to read the compilation error. After all, the dynamic code system is a system that allows us to program C++ code that gets built and linked against the solver at runtime. We are still programming C++, so things like this can and will happen:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
    Making dependency list for source file fixedValueFvPatchFieldTemplate.C
    path/to/your/simulationCase/0/U.boundaryField.inlet: In member function ‘virtual void Foam::velocitySquareInletFixedValueFvPatchVectorField::updateCoeffs()’:
    path/to/your/simulationCase/0/U.boundaryField.inlet:36:39: error: ‘mesh’ was not declared in this scope
    path/to/your/simulationCase/0/U.boundaryField.inlet:36:35: warning: unused variable ‘C’ [-Wunused-variable]
    fixedValueFvPatchFieldTemplate.dep:492: recipe for target 'Make/linux64GccDPOpt/fixedValueFvPatchFieldTemplate.o' failed
    make: *** [Make/linux64GccDPOpt/fixedValueFvPatchFieldTemplate.o] Error 1

    --> FOAM FATAL IO ERROR: 
    Failed wmake "dynamicCode/velocitySquareInlet/platforms/linux64GccDPOpt/lib/libvelocitySquareInlet_76a01ae3375b66215148138e199da3c9b61d9466.so"

    file: path/to/your/simulationCase/0/U.boundaryField.inlet from line 25 to line 30.

        From function codedBase::createLibrary(..)
        in file db/dynamicLibrary/codedBase/codedBase.C at line 213.

    FOAM exiting

This is the compiler letting us know that we can’t simply use the global mesh object within BC code like this: const volVectorField& C = mesh.C();. However, fixing this issue by removing the troublesome line, doesn’t change the compiler error report. There are two options:

  1. Delete the dynamicCode/velocitySquareInlet directory each time a modification is done in the 0/U file.
    • Takes a lot of time if the case is even a bit larger.
  2. Simply continue programming the automatically generated library.
    • You can use git to track the files and that is just perfect.
    • You can get the compiler errors faster – no need to start a simulation.
    • Remember to port the final code back to 0/U or whichever field requires the codedFixedValue BC.

I opt for 2, for obvious reasons. So basically, once your first prototype of the coded object is finished or you have had your share of typing and starting the solver still results with an error, it’s time to start working with the automatically generated dynamicCode library. The implementation file fixedValueFvPatchFieldTemplate.C contains the updateCoeffs member function.

Note
In case you don’t know which member function you are coding with the coded object, search for the code string in the generated code.

Here is what the updateCoeffs member function looks like:

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
    void velocitySquareInletFixedValueFvPatchVectorField::updateCoeffs()
    {
        if (this->updated())
        {
            return;
        }

        if (false)
        {
            Info<<"updateCoeffs velocitySquareInlet sha1: 76a01ae3375b66215148138e199da3c9b61d9466\n";
        }

    //{{{ begin code
        #line 31 "/path/to/your/simulationCase/0/U.boundaryField.inlet"

            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 

            vectorField& field = *this; 

            const volVectorField& C = mesh.C();  // Messed up on purpose here. 

            scalar min = 0.501;  
            scalar max = 0.751; 

            forAll(Cf, faceI)
            {
               if (
                    (Cf[faceI].z() > min) &&
                    (Cf[faceI].z() < max) &&
                    (Cf[faceI].y() > min) &&
                    (Cf[faceI].y() < max) 
                  )
                {
                    field[faceI] = vector(1, 0, 0);
                }
            }
    //}}} end code

        this->fixedValueFvPatchField<vector>::updateCoeffs();
    }

You can see where I have messed the code up on purpose. Sure, all solvers and applications that use the fvMesh/dynamicFvMesh classes build it by including the createMesh.H/crateDynamicFvMesh.H boilerplate files, that instantiate a global mesh objects. However we do not want to make a library aware of a named global object. This will cause compiler errors for all applications that will have named the global mesh object differently. That is why I’m not describing that this can be done by writing extern fvMesh mesh line above the mesh.C() call site, oh, wait… Following the rules covered in the previous section, we dig through the fvPatchField, it’s data members and parent classes to find ways to get the constant reference to the mesh. Here is one way to do that:

1
2
3
    const DimensionedField<vector,volMesh> df = dimensionedInternalField(); 
    const fvMesh& mesh = df.mesh(); 
    const volVectorField& C = mesh.C();  

Another way would look like this

1
2
3
    const fvBoundaryMesh& boundaryMesh = boundaryPatch.boundaryMesh(); 
    const fvMesh& mesh = boundaryMesh.mesh();
    const volVectorField& C = mesh.C(); 

And the code compiles, with a warning of course: warning: unused variable ‘C’ [-Wunused-variable]. That’s because we’re not using the vector field of all cell centers anywhere in our codedFixedValue BC. Remember: use git and port the implemented (updateCoeffs) member function into the code block of the codedFixedValue BC, in the 0/U file – or to whichever field you chose for applying a codedFixedValue BC.

Accessing registered objects

We have covered registering objects in OpenFOAM in the OpenFOAM Primer. The mesh acts as registry of objects with respect to the fields. The object registry allows non-const access to its registered objects. Actually, the approach is an implementation of the Observer/Subject OOD pattern. So the answer to the question is a lookup of any field that is registered to the mesh, using its name. For example, if we were to access a pressure field within the velocity BC, we could code something like this

1
2
    const volScalarField& p = mesh.lookupObject<volScalarField>("p"); 
    // Working with the pressure field.

The rest is straightforward: loop through the boundary fields, access the boundary field with the same patch ID as the codedFixedValue inlet boundary field, and make the velocity values dependent somehow on the pressure. Leave a comment in the comment section below if you have an example of such a BC dependency that you would like to see implemented. If we find the time, we might code it and write about it.

The codedFixedValue square inlet in action

This time we have even made an animation of the codedFixedValue inlet BC 🙂

In this very simple setup, the (min,max) pair is used to define a rectangular inlet shape. This can be extended easily so that temporal variation is included – the motion of the inlet can be added such that the location of the inlet varies with time, the geometrical shape and dimensions of the inlet can be changed (also in time), and so on. This simplest example is only to show how to work with non-uniform boundary field values on a case-by-case basis. Once you get one of these up and running for a simulation case, you can extract the library from the dynamicCode directory and use it as a general purpose (case-independent) boundary condition.

Conclusions

The coded boundary condition makes it possible to prototype case-specific boundary conditions quickly in OpenFOAM. However, as soon as the same condition is applied to two fields and/or multiple simulation cases, we should consider implementing a new boundary condition library to avoid unnecessary source code duplication. If you re-examine the codedFixedValue BC for both the volume fraction field α and the velocity field U, you will notice that they are exactly the same, apart from the tensor rank. When same algorithms work with different types, not associated in any way, the usual way to reduce code duplication would be to code a class template.

This is an extremely simple example – If you have a more application driven / complex / interesting idea for such a codedFixedValue BC, let us know in the comment section below – if we find time we’ll code it and release the code in the blog post repository.

Series Navigation

<< Modify Boundary Fields