• CFD, Fluid Flow, FEA, Heat/Mass Transfer
  • +91 987-11-19-383
  • amod.cfd@gmail.com

New boundary conditions, utility and solver programming in OpenFOAM

OpenFOAM: Programming


This content is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OpenFOAM and OpenCFD trademarks.
One of the strengths of OpenFOAM is that new solvers and utilities can be created by its users with some prior knowledge of the underlying method, physics and programming techniques involved. This tutorial by Tommaso Lucchini from Department of Energy Politecnico di Milano serves to generate this understanding.

OpenFOAM code for the N-S equation:

Note that the specifications fvm and fvc are selected by the user from the fvSchemes dictionary in the system dictionary. Here fvm::laplacian means an implicit finite volume discretization for the Laplacian operator, and similarly for fvm::div for the divergence operator. On the other hand, fvc::grad means an explicit finite volume discretization for the gradient operator. The parentheses ( , ) means a product of the enclosed quantities, including tensor products. Similarly, & denotes dot product of two vectors in OpenFOAM. The code to solve Navier-Stokes equation described below is written for demonstration of fvm and fvc class.

Naier-Stokes Equation

Solve (		
  fvm::ddt(rho,U) + fvm::div(U,U) - fvm::laplacian(mu,U) == f - fvc::grad(p)

Similarly, the statement to solve k-equation in OpenFOAM is given by
 solve (
   fvm::ddt(k) +
   fvm::div(phi, k) -
   fvm::laplacian(nu() + nut, k) ==
   nut*magSqr(symm(fvc)::grad(U)) - fvm::Sp(epsilon/k, k)); 
Here k-equation is expressed as

k-equation for turbulence

Su:Implicit evaluation of the source term in the equation. SuSp(phi, A): Implicit/Explicit source term in equation, depending on sign of phi.
OpenFOAM can be used to solve equations governed by theory of elasticity (solid mechanics) as well.

Solid mechanics equation

This equation can be solved in OpenFOAM by:
fvVectorMatrix UEqn (	
rho * fvm::d2dt2(U) == fvm::laplacian(2*mu + lambda, U)
    + fvc::div
    - (mu + lambda) * gradU
    + lambda * tr(gradU)*I
    + (sigma & gradU)

C++ Syntax and Conventions

C++ is a statically typed [that error checking is done only during compilation], compiled [a separate application called compiler is needed and executable needs to be created], case-sensitive, free-form [unlike FORTRAN or Python where indentation specifies a block of code it does not recognize end of the line as a terminator . Thus, it does not matter where a statement in a line is put] programming language that supports procedural, object-oriented programming. An interpreted language such as Python may also be object-oriented but error checking is done during run-time and no compilation is required.

The :: operator: It is called the scope-resolution operator and by prefixing a type-name with this, it tells compiler to look in the global namespace for the type. In other words, :: links either of a variable, function, class, typedef etc... to a namespace or to a class. If there is no left hand side before :: then it refers to the global namespace. Thus, in 'fvm::laplacian', fvm is a namespace and laplacian is a function (it can also be another namespace of a variable or an onject). Foam: namespace for OpenFOAM, constant: namespace for collection of constants, fv: namespace for finite-volume, fvc (Finite Volume Calculus): namespace of functions to calculate explicit derivatives - given a field creats a new field, fvm (Finite Volume Method): namespace of functions to calculate implicit derivatives returning a matrix.
  • namespace is a method to distinguish two variables or functions having same names. It is analogous to roll numbers to differentiate two students having exactly same name.
  • using namespace Foam::constant; → The 'using' directive here is used to refer to a particular item within a namespace.
  • Foam::scalar Foam::noiseFFT::p0 = 1e-5; → Here the variable p0 in scopes 'Foam::scalar' and 'Foam::noiseFFT' are assigned a value of 1e-5.
  • typedef: create a new name for an existing type. e.g. typedef int single; → single is another name for int, typedef float real; → real is another name for float. Note that the term 'void' represents the absence of type. In case of a 'function' it tells that the function does not return anything.
  • Constants refer to fixed values that the program cannot alter, they are also called literals. Constants can be defined either using #define preprocessor or 'const' keyword.
  • x = &myObj; defines x as pointer to object myObj. x -> read(); accesses the member function read(). e.g. turbulence -> correct(); in $FOAM_SOLVERS/incompressible/pimpleFoam/pimpleFoam.C
  • 'forAll' defined in UList.H is an alias to 'for-loop' defined to march through all entries of a list of objects of any class.
Concepts in OOP (Object-oriented Programming)

Object oriented progamming

Classes in OpenFOAM:
  • 'polyMesh' class holds mesh definition objects whereas polyBoundaryMesh is a list of polyPatches

OpenFOAM root folder structure

  1. src: the core OpenFOAM libraries with *.C and *H files, including include files. Typing 'foam' on the command line will take you to installation folder. Similarly, app = $WM_PROJECT_DIR/applications = cd $FOAM_APP and util = $WM_PROJECT_DIR/applications/utilities. Command alias | grep "FOAM" retrieves following aliases
    • alias app='cd $FOAM_APP'
    • alias lib='cd $FOAM_LIBBIN'
    • alias run='cd $FOAM_RUN'
    • alias sol='cd $FOAM_SOLVERS'
    • alias src='cd $FOAM_SRC'
    • alias tut='cd $FOAM_TUTORIALS'
    • alias util='cd $FOAM_UTILITIES'
  2. $FOAM_APP/solvers/: solvers (such as icoFoam, pimpleFoam, sonicFoam ...) and utilities (such as blockMesh, setFields, topoSet ...). Solvers are grouped into 12 categoryies: basic, compressible, DNS, financial, incompressible, multiphase, combustion, discreteMethods, electromagnetics heatTransfer, lagrangian and stressAnalysis. Move into the folder solver using command 'sol'
  3. $FOAM_SRC/finiteVolume/finiteVolume: this folder contains discretization schemes namely convectionSchemes, ddtSchemes, fv, fvm, d2dt2Schemes, divSchemes, fvc, fvSchemes, fvSolution, gradSchemes, laplacianSchemes, snGradSchemes. Use alias 'foamfv' to move into $FOAM_SRC/finiteVolume directory.
  4. $FOAM_TUTORIALS: sample cases (simulation data) to demonstrate OpenFOAM functionalities, tutorial folder by typing 'tut' on command line
  5. test: cases that test and validate individual models or aspects of OpenFOAM functionality
  6. etc: scripts -foamyHexMeshDict, mesh, meshQualityDict, postProcessing, setConstraintTypes, solvers, surface
  7. doc: documentation such as user guides in PDF format, Doxygen
  8. In addition to these folders, the root folder also contains a file call 'Allwmake' and a directory called 'wmake'. These are 'make' utility for compilations.
  9. $WM_PROJECT_DIR: Full path to locate binary executables of OpenFOAM release
  10. Use command <find -iname "sampleDict"> to find a file named sampleDict. To search a folder named 'pimpleFoam' in folder '$WM_PROJECT_DIR' use command: find $WM_PROJECT_DIR –type d -name "*pimpleFoam*" where *is known as wildcard character. To find a file replace 'd' with 'f'. Similarly, command "find $FOAM_SRC -iname \*fixedValue\*.[CH]" lists all the occurrences of fixedValue.C and fixedValue.H files in folder $FOAM_SRC.
  11. $WM_PROJECT_DIR/src/OpenFOAM/lnInclude: list of *.H and *.c files linked to specific utility

C++ File Structure
  • There needs to be a *.C file with the solver name that calls all of the necessary libraries as well as initializing the time, meshes, and material properties. This file contains the time loop as well as the calls to the individual solvers.
  • The other .H files in the solver directory relate to the initialization and the calls within the time loop.
  • The Make directory includes a file titled options which "contains the full directory paths to locate the header files … and library names" while the file titled "files", lists the entire path and name of the solver once it has been compiled.
  • The difference between compilation of solver and boundary condition lies in the fact that solvers are compiled into executables and boundary conditions are compiled into shared object libraries which can be linked to any solver without need to recompile for every solver.
  • All the source codes related to boundary conditions can be found in the file /src/finiteVolume/fields/fvPatchField/derived. In this folder, you can also check whicl all boundary conditions types are available with the installed version of OpenFOAM.
  • OpenFOAM usage 'wmake' utility for compilation - a reference copy is available in folder '$FOAM_SRC'.
The file structure is shown below.

OpenFOAM solver file heirarchy

FVM solution method

Before further moving into the OpenFOAM program structure, let's look at the governing equations and methods adopted for their numerical solutions.

FVM discretization scheme

The integration over a cell volume need to be converted into surface intergral over the face of the cells using Gauss theorem as explained below.

FVM discretization Gauss Theorem

Let's look at some of the settings in fvScheme dictionary. Note that the variable φ defined above is a generic representation of u, v, w ... and does not represent 'phi' defined below.
ddtSchemes   {               //ddtScheme is a 'class' under namespace 'fv'
    default         Euler;   

gradSchemes  {               //gradScheme is a 'class' under namespace 'fv'
    default         Gauss linear;

divSchemes   {
    default          none;
    div(phi,U)       bounded Gauss linearUpwind grad(U);
    div(phi,k)       bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;       //phi = U.A = Flux field
    div(phi,R)       bounded Gauss upwind;
    div(R)           Gauss linear;
    div(phi,nuTilda) bounded Gauss upwind;       //nuTilda = Spalart-Allmaras model
    div((nuEff*dev2(T(grad(U))))) Gauss linear;  //nuEff = nuLam + nuTurb
dev() calculates the deviatoric part of the matrix: ADEV = A - AHYD = A - 1/3 × tr(A)I. OpenFOAMr versions like 4.x onwards, the incompressible solvers will use dev2() [introduced in version 3 for better stability and convergence] otherwise used for compressible flows. The difference between dev2 and dev is the factor two at some position of the code. dev2() calculates the deviatoric part of the matrix: ADEV = A - 2AHYD = A - 2/3 × tr(A)I.
laplacianSchemes {
    default         Gauss linear corrected;

interpolationSchemes {
    default         linear;
Solvers operate on the lduMatrix class located in folder FOAM_SRC/OpenFOAM/matrices/lduMatrix/ - lduMatrix is a matrix class in which the coefficients of matrix [A] in linear equation [A]{x}={b} are stored as three arrays. One for the lower triangle (l), one for the the diagonal of the matrix (d) and the third array for uppper triangle (u) - hence known as lduMatrix. Inside of lduMatrix.H and lduMatrix.C, the member functions - const scalarField & lower () const, const scalarField & diag () const and const scalarField & upper () const are defined which returns the list of cofficients for the lower, diagonal and upper part of the matrix.

fvSolution: pimpleFoam This dictionary defines the selection of linear solver, smoothing algorithm and convergence tolerance for inner iterative process.

  1. BICCG - Diagonal Incomplete LU preconditioned BiCG solver
  2. diagonalSolver - diagonal solver for both symmetric and asymmetric problems
  3. GAMG - generally known as Geometric-agglomerated Algebraic MultiGrid solver (named Generalised geometric-Algebraic MultiGrid in the user manual)
  4. ICC - Incomplete Cholesky preconditioned Conjugate Gradients solver
  5. PBiCG - Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices using a run-time selectable preconditioner
  6. PCG - Preconditioned conjugate gradient solver for symmetric lduMatrices using a run-time selectable preconditiioner
  7. smoothSolver - Iterative solver using smoother for symmetric and asymmetric matrices which uses a run-time selected smoother
Preconditioner and smoother: Preconditioners improves the convergence of the origina system [A]{x}={b} using a matrix [M] and resulting in a new system [M]-1[A]{x} = [M]-1{b}. Smoother on the other hand is to reduce the dependency of number of iterataions on mesh type.

solvers {
    p     {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        smoother         DICGaussSeidel;

    pFinal     {
        relTol          0;

    "(U|k|epsilon)"     {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;

    "(U|k|epsilon)Final"     {
        relTol          0;

    nNonOrthogonalCorrectors 0;  
    nOuterCorrectors    100; //Outer Loops (Pressure - Momentum Correction )
    nCorrectors         2;   //Inner Loops (Pressure Correction )

Utility statements

  • Centre of faces: Cf[faceI].x() / patchName.Cf()[i][0], Cf[faceI].y() / patchName.Cf()[i][1], Cf[faceI].z() / patchName.Cf()[i][2]
  • Access time: this->db().time().value()
  • Loop over all faces of a patch: forAll(Cf, faceI) { field[faceI] = ... }

Source code: pimpleFoam Large time-step transient solver for incompressible, flow using the PIMPLE (merged PISO-SIMPLE) algorithm. Sub-models include: turbulence modelling, i.e. laminar, RAS or LES, run-time selectable finite volume options, e.g. MRF, explicit porosity. The lines starting with .. are the content of 'included' header files. Some changes in code structure hs been made to reduce number of lines.

/* Make all the FV machinery or templates available to this code: this header
file is located at $FOAM_SRC/finiteVolume/cfdTools/general/include           */
#include "fvCFD.H"

#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "IObasicSourceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])   {
// Check case folder: $FOAM_SRC/OpenFOAM/include
   #include "setRootCase.H"
      ..  #include "listOptions.H"
      ..  Foam::argList args(argc, argv);
      ..   if (!args.checkRootCase()) {
      ..   Foam::FatalError.exit();  //exit() is a member function of class 'FatalError' in namespace 'Foam'
      ..  }
      ..  #include "listOutput.H"

//  Create time directories: $FOAM_SRC/OpenFOAM/include
    #include "createTime.H"
      ..  Foam::Info<< "Create time\n" << Foam::endl;
      ..  Foam::Time runTime(Foam::Time::controlDictName, args);

//  Read mesh: $FOAM_SRC/OpenFOAM/include
    #include "createMesh.H"
      ..  Foam::Info
      ..    << "Create mesh for time = "
      ..    << runTime.timeName() << Foam::nl << Foam::endl;
      ..  Foam::fvMesh mesh 
      ..    (Foam::IOobject 
      ..      (Foam::fvMesh::defaultRegion,
      ..       runTime.timeName(),
      ..       runTime,
      ..       Foam::IOobject::MUST_READ  //Must read mesh every timestep
      ..      )
      ..    );

//  Create read U, p, k, epsilon ... dictionaries and initalize fields
//  File is solver specific: located in $FOAM_APP/solvers/incompressible/pimpleFoam	  
    #include "createFields.H"
      ..  Info<< "Reading field p\n" << endl;
      ..  volScalarField p  (      //define p volume scalar field
      ..    IOobject  (
      ..      "p",
      ..      runTime.timeName(),
      ..      mesh,
      ..      IOobject::MUST_READ,
      ..      IOobject::AUTO_WRITE
      ..    ),
      ..   mesh
      ..  );

      ..  Info<< "Reading field U\n" << endl;
      ..  volVectorField U  (        //define U volume vector field
      ..    IOobject (
      ..      "U",
      ..      runTime.timeName(),    // look for dictionary "U"
      ..      mesh,
      ..      IOobject::MUST_READ,
      ..      IOobject::AUTO_WRITE
      ..    ),
      ..   mesh
      ..  );

      ..  #include "createPhi.H"    // create flux field phi at cell faces

      ..  label pRefCell = 0;
      ..  scalar pRefValue = 0.0;
      ..  setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);

      ..  singlePhaseTransportModel laminarTransport(U, phi);

      ..  autoPtr turbulence  (
      ..    incompressible::turbulenceModel::New(U, phi, laminarTransport)
      ..  );
    //Declare and initialise the cumulative continuity error.
    #include "initContinuityErrs.H"
      ..  #ifndef initContinuityErrs_H
      ..  #define initContinuityErrs_H
      ..  scalar cumulativeContErr = 0;
      ..  #endif

    pimpleControl pimple(mesh);
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    Info<< "\nStarting time loop\n" << endl;  //Time loop as per controlDict
    while (runTime.run())  {                  //While time < endTime
        #include "readTimeControls.H"  // Read adjustTimeStep, maxCo, maxDeltaT
        #include "CourantNo.H"         // Get the Mean and max Courant Numbers 

        /*$FOAM_SRC/finiteVolume/cfdTools/general/include:  Reset the timestep 
          to maintain a constant maximum courant Number. Reduction of time-step 
          is immediate, but increase is damped to avoid unstable oscillations */
        #include "setDeltaT.H"        // Get the deltaT defined in controlDict
          ..  if (adjustTimeStep) {
          ..    scalar maxDeltaTFact = maxCo/(CoNum + SMALL);
          ..    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
          ..    runTime.setDeltaT (
          ..      min (
          ..        deltaTFact*runTime.deltaTValue(),
          ..        maxDeltaT
          ..      )
          ..    );
          ..  }


        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity PIMPLE corrector loop:  as many outer correctors 
        // as nOuterCorrectors  and as many inner correctors as nCorrectors defined
        // in fvSolution dictionary.
        while (pimple.loop())   {    
            #include "UEqn.H"      //Solution of the momentum equation

            // --- Pressure corrector loop
            while (pimple.correct()) {
                #include "pEqn.H"  //Include the pressure Equation
            if (pimple.turbCorr())  {
        //Print on the screen information - computational and clock time
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    Info<< "End\n" << endl;
    return 0;
// ************************************************************************* //


This content is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OpenFOAM and OpenCFD trademarks.
  1. OpenFOAM's basic solvers for linear systems of equations Solvers, preconditioners, smoothers by Tim Behrens, February 18, 2009
  2. OpenFOAM Programmers Guide
  3. Mathematics, Numerics, Derivations and OpenFOAM by Tobias Holzmann at Holzmann CFD
  4. Dictionaries in tutorial folders of OpenFOAM
  5. https://www.youtube.com/watch?v=L94iYGvZr9Q
  6. https://www.youtube.com/watch?v=2FKoMY91SdM&t=3231s
Contact us
Disclaimers and Policies

The content on CFDyna.com is being constantly refined and improvised with on-the-job experience, testing, and training. Examples might be simplified to improve insight into the physics and basic understanding. Linked pages, articles, references, and examples are constantly reviewed to reduce errors, but we cannot warrant full correctness of all content.