CFD Online Discussion Forums - Reduce() and Size() functions

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Reduce() and Size() functions (https://www.cfd-online.com/Forums/openfoam-programming-development/176851-reduce-size-functions.html)

Reduce() and Size() functions

Dear Fellows

I have few questions regarding the following snippet:

Code:

scalar ULeft = 0.0;

label leftI = mesh().boundaryMesh().findPatchID("left");

const fvPatchVectorField& fvp = U.boundaryField()[leftI];

if (fvp.size())

{

ULeft = fvp[0].x();

}

reduce(ULeft, maxOp<scalar>());

After labeling leftI, the 3rd line puts the velocity of the "left" ptach into fvp "object".
I assume the size of fvp is non-zero, so the "if" loop executes.
My first question is on

Code:

ULeft = fvp[0].x();

.
here, fvp is an array of vectors, since it keeps the velocity vectors of a patch, but the above command just stores the first component ".x()" of first vector "[0]" in the array of vectors. Am I correct?

For example assume
fvp = { (1 2 3) (4 5 6) (7 8 9)}
then fvp[0].x means 1
then why using the

Code:

reduce(ULeft, maxOp<scalar>());

command?

I assume that reduce() function with maxOp<scalar> finds the biggest scalar in an array of scalars, but here ULeft just has one member!!! (in my example 1)

Am I missing sth?:confused::confused:

Regards

Hey,
i assume you have this snippet from the potentialFoam/cylinder tutorial ?!

In the 0/U file you can see that they use a uniformFixedValue velocity (1 0 0) on the left boundary condition. Because you have the same velocity on each boundary face of the left BC it´s ok to just ask for a single element of fvp. Because you just have the x-component of velocity it´s also ok to just ask for fvp[0].x() ! Of course this is not very generic or reusable for other cases.

When you decompose the domain (parallel computing with MPI) you just get the correct velocity from the process (or processes) which compute the subdomain containing the left BC. Other processes will return a null velocity. The reduce operation searches for the maximum value in your parallel processes on the left BC which will be a 1. By not using the reduce operation there is a risk that the master process is not computing the part which contains the left-BC and hence returns a null velocity.

Hey Traction,

Yeah, I was reading the code snippet available in OpenFOAM Programmers guide 3.0.1.
Thanks for your clarification. Since I did not go through the whole sample case.:)

I just read the code. But with your suggestions, my questions regarding the code snippet is responded.

Dear Fellows

I have another question regarding this snippet inside O.F. programmers guide version 3 (page 49) .

Code:

dimensionedScalar uInfX

(

"uInfx",

dimensionSet(0, 1, -1, 0, 0),

ULeft

);

Why the dimensionSet class has just 5 components? I didn't know that we could write scalars or vectors with reduced number of units in OF. :confused::confused:

Code:

volScalarField error("error", mag(U-UA)/mag(UA));



            Info<<"Writing relative error in U to " << error.objectPath()

                << endl;

I executed this sample and the code produced an error file in 0 folder. Can some one tell me how does this "objectPath()" member function work?

Here, just a 0 folder was present, However in case of more time steps what would happen?

Hey babakflame,

Quote:

Why the dimensionSet class has just 5 components? I didn't know that we could write scalars or vectors with reduced number of units in OF. :confused::confused:

A good practice in OpenFOAM is to have a look at the specific source files. In this case the answer is presented in src/OpenFOAM/dimensionSet/dimensionSet.H and dimensionSet.C: there is a constructor for all seven dimensions and another constructor for just using the first five dimensions (and setting the leftover dimensions to zero by default). You can find these files for example by searching for dimensionSet.H in your installation or the repository.

Quote:

OpenFOAM stores results in the time directories. Hence the error.write() function stores the volScalarField error in the current time directory.

The OpenFOAM C++ Documentation says that objectPath() is a function to return the complete object path and the name of the object (in this case of the object error which is located in the current time directory). In this case it just produces an output on your screen / log file to find the location of the written error file.

Keep in mind that potentialFoam is independent of the time. If there would be a suitable solver with time dependency OpenFOAM would write the calculation results for the fields "UA" and "error" for every timeStep.

Hey Traction,

Thanks bro for your comprehensive response.:cool:

Quote:

Originally Posted by babakflame (Post 616012)

Code:

volScalarField error("error", mag(U-UA)/mag(UA));



            Info<<"Writing relative error in U to " << error.objectPath()

                << endl;

The IOobject-properties of a newly created object error is supposed to be the same as those of the object mag(U-UA)/mag(UA). mag(U-UA)/mag(UA) is constructed as a result of perfoming successive function (mag) and operation (operator-, operator/) calls on volScalarFields.
For example U-UA is computed with this code:

Code:

  //-- src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C

  ...

  616 template                                                                       \

  617 <class Type1, class Type2, template<class> class PatchField, class GeoMesh>    \

  618 tmp                                                                            \

  619 <                                                                              \

  620     GeometricField<typename product<Type1, Type2>::type, PatchField, GeoMesh>  \

  621 >                                                                              \

  622 operator op                                                                    \

  623 (                                                                              \

  624     const GeometricField<Type1, PatchField, GeoMesh>& gf1,                     \

  625     const GeometricField<Type2, PatchField, GeoMesh>& gf2                      \

  626 )                                                                              \

  627 {                                                                              \

  628     typedef typename product<Type1, Type2>::type productType;                  \

  629     tmp<GeometricField<productType, PatchField, GeoMesh>> tRes                 \

  630     (                                                                          \

  631         new GeometricField<productType, PatchField, GeoMesh>                   \

  632         (                                                                      \

  633             IOobject                                                           \

  634             (                                                                  \

  635                 '(' + gf1.name() + #op + gf2.name() + ')',                     \

  636                 gf1.instance(),                                                \

  637                 gf1.db(),                                                      \

  638                 IOobject::NO_READ,                                             \

  639                 IOobject::NO_WRITE                                             \

  640             ),                                                                 \

  641             gf1.mesh(),                                                        \

  642             gf1.dimensions() op gf2.dimensions()                               \

  643         )                                                                      \

  644     );                                                                         \

  645                                                                                \

  646     Foam::opFunc(tRes.ref(), gf1, gf2);                                        \

  647                                                                                \

  648     return tRes;                                                               \

  649 }

  ...

Thanks sergei.

I read this part that you send me.:)

It seems that both add and subtract operations are done by this snippet. (overloading operators)

Am I correct?

Yes, that code is a part of a macro which generates both

Code:

operator +                                                                

(                                                                              

const GeometricField<Type1, PatchField, GeoMesh>& gf1,                     

const GeometricField<Type2, PatchField, GeoMesh>& gf2

)

and

Code:

operator -                                                                

(                                                                              

const GeometricField<Type1, PatchField, GeoMesh>& gf1,                     

const GeometricField<Type2, PatchField, GeoMesh>& gf2

)