Notes for OpenFOAM | Zifeng Weng

This page will receive periodic updates.

writePrecision

Use finer precision to reduce truncated error when re-run the simulation from a given time step.
In AMR, large truncated error might lead to wrong point coordinates and zero cell volume.

Load balanced adaptive mesh refinement

The library doesn’t compatible with openmpi-3.1.0.

Implementation of fvc::div

fvc::div functions were realized in the following ways

// in src/finiteVolume/finiteVolume/fvc/fvcDiv.C
// fvc::div(vf,name) --- vf: volume field; name: scheme
return fv::divScheme<Type>::New
(
    vf.mesh(), vf.mesh().divScheme(name)
).ref().fvcDiv(vf);

// fvc::div(flux,vf,name) --- flux: surfaceScalarField&
return fv::convectionScheme<Type>::New
(
    vf.mesh(),
    flux,
    vf.mesh().divScheme(name)
).ref().fvcDiv(flux, vf);

There are two steps, fisrt, construct an object [type: tmp<convectionScheme<Type>>] of divScheme or convectionScheme; second, call fvcDiv function.

divScheme or convectionScheme

 // in src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C
 // divScheme
     typename IstreamConstructorTable::iterator cstrIter = IstreamConstructorTablePtr_->find(schemeName);
     return cstrIter()(mesh, schemeData);
 // convection scheme
     typename IstreamConstructorTable::iterator cstrIter = IstreamConstructorTablePtr_->find(schemeName);
     return cstrIter()(mesh, faceFlux, schemeData);

fvcDiv
1. For divergence
  
  Only one type is available: gaussDivScheme
  
  The essential code is c++ // fvcDiv function in src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C fvc::surfaceIntegrate ( this->tinterpScheme_().dotInterpolate(this->mesh_.Sf(), vf) )
  
  tinterpScheme_ was initialized in father class, c++ tinterpScheme_(surfaceInterpolationScheme<Type>::New(mesh, is)) // is: Istream for divScheme
2. For convection
  
  Three types were available
  - gaussConvectionScheme
  - boundedConvectionScheme
  - multivariateGaussConvectionScheme
  We only check gaussConvectionScheme, the essential code is
```
 fvc::surfaceIntegrate(flux(faceFlux, vf))

 // flux(faceFlux, vf) is defined as
     faceFlux*interpolate(faceFlux, vf);   
 // interpolate(faceFlux, vf) is defined as      
     return tinterpScheme_().interpolate(vf);
```
  during the construction, tinterpScheme_ is a pointer for the specified surface interpolation scheme.
```
 tinterpScheme_
 (
     surfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
 )
```
  The exact interpolation methods, like linear, WENO, were realized in corresponding class.

Implementation of fvm::div

fvm::div calculates the matrix for the divergence of the given field and flux. It only provides convection scheme. The code is

// fvm::div(flux,vf,name) is
return fv::convectionScheme<Type>::New
(
    vf.mesh(),
    flux,
    vf.mesh().divScheme(name)
)().fvmDiv(flux, vf);

In gaussConvectionScheme, the function fvmDiv calls several functions via the surface interpolation object tinterpScheme_, including

weights(vf)
corrected()
correction(vf)

// return type: tmp<fvMatrix<Type>> fvmDiv
    tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);
    const surfaceScalarField& weights = tweights();

    ...

    fvm.lower() = -weights.primitiveField()*faceFlux.primitiveField();
    fvm.upper() = fvm.lower() + faceFlux.primitiveField();
    fvm.negSumDiag();

    ...

    if (tinterpScheme_().corrected())
    {
        fvm += fvc::surfaceIntegrate(faceFlux*tinterpScheme_().correction(vf));
    }

    return tfvm;

Implementation of fvc::grad

For volume mesh fvc::grad(vf,name)

// fvc::grad(vf,name)
return fv::gradScheme<Type>::New
(
    vf.mesh(),
    vf.mesh().gradScheme(name)
)().grad(vf, name);

Here, several gradient schemes are available. Again, we only focus on gaussGrad.

For surface mesh

// fvc::grad(ssf,name)
return fv::gaussGrad<Type>::gradf(ssf, "grad(" + ssf.name() + ')');

grad function of gaussGrad (realized in father class gradScheme)
- cache
- calcGrad(vsf, name)
  - First, call gradf(tinterpScheme_().interpolate(vsf), name), tinterpScheme_ is created without face flux, i.e. surfaceInterpolationScheme<Type>::New(mesh, is)
  - Then call correctBoundaryConditions(vsf, gGrad);
gradf function of gaussGrad to update

Implementation of surfaceInterpolationScheme

The implementation of fvc::div, fvm::div and fvc::grad all reply on interpolation. In a summary, following functions should be provided in each intepolation method (or class),

// fvc::div(vf)
dotInterpolate
// fvc::div(flux,vf), fvc::grad
interpolate
// fvm::div(flux,vf)
weights(vf)
corrected()
correction(vf)

Here we use linear interpolation as an example,

coming soon

Implementation of mapFields

The mapFields is realized with two steps,

Map mesh via meshToMesh0 class

 meshToMesh0 meshToMesh0Interp
 (
     meshSource,
     meshTarget,
     patchMap, // neglected for consistent map
     cuttingPatches
 );

Map fields found in the target case via template class MapConsistentVolFields

 MapVolFields<scalar> // type includes vector, sphericalTensor, symmTensor, tensor
 (
     objects,
     meshToMesh0Interp,
     mapOrder,
     CombineOp<scalar>()
 );

The first step takes much longer time over the whole mapFields execution especially for large mesh size.

meshToMesh0 class

The meshToMesh0 class finds mapping between source and target mesh via calcAddressing() function, which mainly has four steps

Visit all boundaries and mark the cell next to the boundary.

 // Set reference to boundary
 const polyPatchList& patchesFrom = fromMesh_.boundaryMesh();

 forAll(patchesFrom, patchi)
 {
     // Get reference to cells next to the boundary
     const labelUList& bCells = patchesFrom[patchi].faceCells();

     forAll(bCells, facei)
     {
         boundaryCell[bCells[facei]] = true;
     }
 }

construct octree

The octree is used to search cell that containing a given point

 indexedOctree<treeDataCell> oc
 (
     treeDataCell(false, fromMesh_, polyMesh::CELL_TETS),
     shiftedBb,      // overall bounding box
     8,              // maxLevel
     10,             // leafsize
     6.0             // duplicity
 );

calculate addressing for internal points and boundary points (Most time consuming step)
```
 // realized in calculateMeshToMesh0Addressing.C
 cellAddresses
 (
     cellAddressing_,
     toMesh_.cellCentres(),
     fromMesh_,
     boundaryCell,
     oc
 );
```
In this step, there are three levels of search. Both level 2 and 3 use octree for search which is the major computation overhead.

Level 1: Simple neighbour array search.

It starts from a cell zero, searches its neighbours and finds one which is nearer to the target point than the current position.

The location of the “current position” is reset to that cell and search through the neighbours continues. The search is finished when all the neighbours of the cell are farther from the target point than the current cell

For the nearest cell to the target point is found in the source mesh, check if it is target point is inside the cell (the finite volume element) using polyMesh::pointInCell

Level 2 Handle boundary points

Level 3 Handle points in other regions of the doamin
simplification for mapping translated mesh

In most of my preprocessing, the target mesh the translated source mesh in one direction, in this case, the Level 2 and 3 search can be neglected. The octree construction is also neglected.