# fvc::div for surfaceVectorField

 Register Blogs Members List Search Today's Posts Mark Forums Read

 March 16, 2011, 05:18 fvc::div for surfaceVectorField #1 Member   Artem Shaklein Join Date: Feb 2010 Location: Russia, Izhevsk Posts: 43 Rep Power: 9 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! Last edited by ARTem; March 16, 2011 at 11:10.

January 9, 2017, 08:21
#2
New Member

toboto
Join Date: Jun 2016
Posts: 13
Rep Power: 3
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.

Last edited by toboto; January 9, 2017 at 09:52. Reason: Add additional information

 January 14, 2017, 17:47 #3 Senior Member     Sergei Join Date: Dec 2009 Posts: 195 Rep Power: 10 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 tmp< GeometricField< typename innerProduct< vector, Type >::type, fvPatchField, volMesh > > div ( const GeometricField< Type, fvPatchField, volMesh > & vf, const word &name ) Code: template tmp< GeometricField< Type, fvPatchField, volMesh > > div ( const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name ) Code: template 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 76 tmp 77 < 78 GeometricField 79 < 80 typename innerProduct::type, fvPatchField, volMesh 81 > 82 > 83 div 84 ( 85 const GeometricField& vf, 86 const word& name 87 ) 88 { 89 return fv::divScheme::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 43 tmp 44 < 45 GeometricField 46 ::type, fvPatchField, volMesh> 47 > 48 gaussDivScheme::fvcDiv 49 ( 50 const GeometricField& vf 51 ) 52 { 53 tmp 54 < 55 GeometricField 56 ::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 43 void surfaceIntegrate 44 ( 45 Field& ivf, 46 const GeometricField& 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& 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& 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 ::type, fvPatchField, volMesh> It is a GeometricField with elements of type innerProduct::type : Code: 88 template 89 class innerProduct 90 { 91 public: 92 93 typedef typename typeOfRank 94 < 95 typename pTraits::cmptType, 96 int(pTraits::rank) + int(pTraits::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 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 42 class typeOfRank 43 { 44 public: 45 46 typedef Cmpt type; 47 }; Code: 129 template 130 class typeOfRank 131 { 132 public: 133 134 typedef Vector type; 135 }; Code: 178 template 179 class typeOfRank 180 { 181 public: 182 183 typedef Tensor 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. qjh888 and mengweilm425 like this. Last edited by Zeppo; January 15, 2017 at 06:21.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post ivan_cozza OpenFOAM Running, Solving & CFD 2 February 6, 2010 07:09 johanna OpenFOAM Programming & Development 2 August 31, 2009 07:12 su_junwei OpenFOAM Running, Solving & CFD 6 October 13, 2008 07:09

All times are GMT -4. The time now is 21:10.