CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   fvc::div for surfaceVectorField (https://www.cfd-online.com/Forums/openfoam-programming-development/86164-fvc-div-surfacevectorfield.html)

 ARTem March 16, 2011 05:18

fvc::div for surfaceVectorField

Greetings !

I want to calculate divergence of velocity.
Firstly, I have volVectorField U. After some manipulations with U, I obtain surfaceScalarField Usurface. Next, i try to calculate divergence: fvc::div(Usurface), but this action leads to unexpected result.
For example, when i compute div of mesh.Sf() (dimension [m^2]), i recieve volVectorField (dimension [m^-1]) instead of volScalarField (dimension [m^1]).
But computation grad of mag(mesg.Sf()) (surfaceScalarField, dimension [m^2]) leads to volVectorField with dimension [m^1]. This is right,

There are my questions:
-Why fvc::grad(surfaceScalarField) gives correct result, but fvc::div(surfaceVectorField) - incorrect? May be i make mistake in calculations?
-Can i calculate divergence of surfaceVectorField using fvc::div?

Thank you!
Best regards!

 toboto January 9, 2017 08:21

I've exactly the same problem. When I interpolate velocity field to the surface I got a surfaceVectorField. When I try to calculate the divergence of this field I got a volVectorField while I expect a volScalarField.

According to the programmers' manual:
Quote:
 The fvc::div function can take as its argument either a surfaceField, in which case φf is speciﬁed directly, or a volField which is interpolated to the face by central diﬀerencing
Update:
I found in the source code of fvcDiv.H that if the fvc::div is applied to a volField it will return a field of one rank lower
Code:

```    GeometricField <typename innerProduct<vector, Type>::type, fvPatchField, volMesh> > div     (         const GeometricField<Type, fvPatchField, volMesh>&,         const word& name     );```
While if fvc::div is applied to a scalarField it will return a field of the same rank!!!
Code:

```tmp<GeometricField<Type, fvPatchField, volMesh>> div ( const GeometricField<Type, fvsPatchField, surfaceMesh>& );```

Any explanation is highly appreciated.
Best regards.

 Zeppo January 14, 2017 17:47

If you check source code you'll see there're 14 overloaded functions fvc:div() for several combinations of input parameters. This set of functions can be broken down into three groups with one "primary" function within each group being called by other functions in the same group. Those primary functions are declared as:

Code:

```template<class Type > tmp< GeometricField< typename innerProduct< vector, Type >::type, fvPatchField, volMesh > > div (     const GeometricField< Type, fvPatchField, volMesh > & vf,     const word &name )```
Code:

```template<class Type > tmp< GeometricField< Type, fvPatchField, volMesh > > div (     const surfaceScalarField &flux,     const GeometricField< Type, fvPatchField, volMesh > &vf,     const word &name )```
Code:

```template<class Type > tmp< GeometricField< Type, fvPatchField, volMesh > > div (     const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf )```
The first function on the list is the most general out of three:

Code:

```75 template<class Type> 76 tmp 77 < 78    GeometricField 79    < 80        typename innerProduct<vector, Type>::type, fvPatchField, volMesh 81    > 82 > 83 div 84 ( 85    const GeometricField<Type, fvPatchField, volMesh>& vf, 86    const word& name 87 ) 88 { 89    return fv::divScheme<Type>::New 90    ( 91        vf.mesh(), vf.mesh().divScheme(name) 92    )().fvcDiv(vf); 93 }```
It selects (through the Runtime Selection Mechanism) one of the available div schemes with name "name" and calls its member function fvcDiv. Currently the only one available div scheme in OpenFOAM is Gauss scheme (see Divergence theorem). Its fvcDiv is:

Code:

```42 template<class Type> 43 tmp 44 < 45    GeometricField 46    <typename innerProduct<vector, Type>::type, fvPatchField, volMesh> 47 > 48 gaussDivScheme<Type>::fvcDiv 49 ( 50    const GeometricField<Type, fvPatchField, volMesh>& vf 51 ) 52 { 53    tmp 54    < 55        GeometricField 56        <typename innerProduct<vector, Type>::type, fvPatchField, volMesh> 57    > tDiv 58    ( 59        fvc::surfaceIntegrate 60        ( 61            this->mesh_.Sf() & this->tinterpScheme_().interpolate(vf) 62        ) 63    ); 64 65    tDiv().rename("div(" + vf.name() + ')'); 66 67    return tDiv; 68 }```
Firstly, it interpolates your volume field (defined in cell centres) onto a surface field (defined in cell face centres) :

Code:

`this->tinterpScheme_().interpolate(vf)`
Secondly, it makes a dot (inner) product (&) of the mesh face vectors and the interpolated field. It is that operation which reduces the rank! :

Code:

`this->mesh_.Sf() & this->tinterpScheme_().interpolate(vf)`
This dot product can be thought of as an operation of getting surface flux of the field through the cell faces.

Lastly, it loops over the cells and for each cell it sums the fluxes through the cell's faces. This total sum of fluxes for each cell is then divided by a cell's volume

Code:

```42 template<class Type> 43 void surfaceIntegrate 44 ( 45    Field<Type>& ivf, 46    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf 47 ) 48 { 49    const fvMesh& mesh = ssf.mesh(); 50 51    const labelUList& owner = mesh.owner(); 52    const labelUList& neighbour = mesh.neighbour(); 53 54    const Field<Type>& issf = ssf; 55 56    forAll(owner, facei) 57    { 58        ivf[owner[facei]] += issf[facei]; 59        ivf[neighbour[facei]] -= issf[facei]; 60    } 61 62    forAll(mesh.boundary(), patchi) 63    { 64        const labelUList& pFaceCells = 65            mesh.boundary()[patchi].faceCells(); 66 67        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi]; 68 69        forAll(mesh.boundary()[patchi], facei) 70        { 71            ivf[pFaceCells[facei]] += pssf[facei]; 72        } 73    } 74 75    ivf /= mesh.Vsc(); 76 }```
The type of an object returned from fvcDiv is

Code:

```45    GeometricField 46    <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>```
It is a GeometricField with elements of type innerProduct<vector, Type>::type :

Code:

```88 template<class arg1, class arg2> 89 class innerProduct 90 { 91 public: 92 93    typedef typename typeOfRank 94    < 95        typename pTraits<arg1>::cmptType, 96        int(pTraits<arg1>::rank) + int(pTraits<arg2>::rank) - 2 97    >::type type; 98 };```
innerProduct is a class template. type nested in innerProduct defines the type of an object which inner product (basically, dot operation) of two objects (of classes arg1 and arg2) results in. typeOfRank is an example of a well known design pattern called "Trait". It can be instantiated with a specific value parameter rank:

Code:

```44 template<class Cmpt, int rank> 45 class typeOfRank 46 {};```
You just provide it with a rank and it tells you of what type an object of this rank is. OpenFOAM defines three partial specialisations of class template typeOfRank with template parameters (ranks) 0, 1, 2 for types scalar, vector and tensor, accordingly

Code:

```41 template<class Cmpt> 42 class typeOfRank<Cmpt, 0> 43 { 44 public: 45 46    typedef Cmpt type; 47 };```
Code:

```129 template<class Cmpt> 130 class typeOfRank<Cmpt, 1> 131 { 132 public: 133 134    typedef Vector<Cmpt> type; 135 };```
Code:

```178 template<class Cmpt> 179 class typeOfRank<Cmpt, 2> 180 { 181 public: 182 183    typedef Tensor<Cmpt> type; 184 };```
pTraits is ... well another example of design pattern "Trait". If you specialise it for a specific type it can tell you the rank of any object of that type. Of course, in OpenFOMA it is allready specialised for scalar, vector and tensor.

The other two "primary" fvc:div functions mentioned in the very beginning of this very lengthy post can be dissected in the similar way.

If something is unclear please feel free to ask.

 All times are GMT -4. The time now is 12:11.