[ARTem]

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]

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:

The fvc::div function can take as its argument either a surface<Type>Field, in which case
φf is specified directly, or a vol<Type>Field which is interpolated to the face by central differencing

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

    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!!!

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

Any explanation is highly appreciated.
Best regards.

[Zeppo]

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:

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

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

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:

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:

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) :

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! :

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

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

45     GeometricField
46     <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>

It is a GeometricField with elements of type innerProduct<vector, Type>::type :

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:

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

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

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

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.

[alexvaleije]

Hi,

I'm stucked with this kind of problem also, and I did not understand really well the answer (maybe because my knowledge about the programming is not very good).

I want to do the divergence of a surfaceVectorField (it actually is a volVector, but i can convert it) to add it to a ddt term, which is a volScalarField, but the divergence is not converting the vector into a scalar.

Can somebody explain to me what I'm doing wrong??

Thank u very much
Alex

[alexvaleije]

Originally Posted by alexvaleije

Hi,

I'm stucked with this kind of problem also, and I did not understand really well the answer (maybe because my knowledge about the programming is not very good).

I want to do the divergence of a surfaceVectorField (it actually is a volVector, but i can convert it) to add it to a ddt term, which is a volScalarField, but the divergence is not converting the vector into a scalar.

Can somebody explain to me what I'm doing wrong??

Thank u very much
Alex

Never mind, I realized I was doing another mistake I was unable to see, and once corrected, I was able to permorm the operation

Regards,
Alex