Initialization of U by code
 

Hi all!. I'm using scalarTransportFoam to solve de transport equation. It's running, but i have a problem. I need to set the initial field of U. I know that there are several utilities to do similar things like funkySetFields, but in this case, i think that i need to write some code.
The mesh is simply a tube, and i want to set up a rigid rotation as a velocity field. I wrote a code using the vector operations provided by OpenFOAM but it crashes with i try to compile it.
Then i tryed putting a portion of H. Jasak code in a .H file, as follows:

USetUp_Jasak.H

// Do cells
const volVectorField& centres = mesh.C();
point origin(1, 1, 0.05);
vector axis(0, 0, -1);
U = axis ˆ (centres.internalField() - origin);
U.write();


This code have no sense to my purposes, but it's only for trying to compile.

Then i added this line in scalarTransportFoam.C file:

int main(int argc, char *argv[])
{

# include "setRootCase.H"

# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "USetUp_Jasak.H" <<<<<<<<<<----------------------

...

It implies that, next U was created and initialized by reading U in /0 directory, the code overwrite the values with the new ones.

When i try to compile, i obtain this:

santiago@elduque:~/OpenFOAM/santiago-1.5/applications/myScalarTransportFoam$ wmake
SOURCE=myScalarTransportFoam.C ; g++ -m32 -Dlinux -DDP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -O3 -DNoRepository -ftemplate-depth-40 -I/home/santiago/OpenFOAM/OpenFOAM-1.5/src/finiteVolume/lnInclude -IlnInclude -I. -I/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude -I/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OSspecific/Unix/lnInclude -fPIC -pthread -c $SOURCE -o Make/linuxGccDPOpt/myScalarTransportFoam.o
In file included from myScalarTransportFoam.C:47:
USetUp_Jasak.H: In function ‘int main(int, char**)’:
USetUp_Jasak.H:5: error: no match for ‘operator=’ in ‘U = Foam::operator^(const Foam::VectorSpace<Form, Cmpt, nCmpt>&, const Foam::tmp<Foam::Field<Type> >&) [with Form = Foam::Vector<double>, Cmpt = double, int nCmpt = 3, Type = Foam::Vector<double>](((const Foam::tmp<Foam::Field<Foam::Vector<double> > >&)((const Foam::tmp<Foam::Field<Foam::Vector<double> > >*)(& Foam::operator-(const Foam::UList<T>&, const Foam::VectorSpace<Form, Cmpt, nCmpt>&) [with Type = Foam::Vector<double>, Form = Foam::Vector<double>, Cmpt = double, int nCmpt = 3](((const Foam::VectorSpace<Foam::Vector<double>, double, 3>&)((const Foam::VectorSpace<Foam::Vector<double>, double, 3>*)(& origin.Foam::Vector<double>::<anonymous>))))))))’
/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricField.C:1032: note: candidates are: void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=(const Foam::GeometricField<Type, PatchField, GeoMesh>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricField.C:1057: note: void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=(const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricField.C:1093: note: void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=(const Foam::dimensioned<Type>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
make: *** [Make/linuxGccDPOpt/myScalarTransportFoam.o] Error 1

Any clues?, thanks in advance. Santiago.

In the line:

U = axis ˆ (centres.internalField() - origin);

You have U as a volVectorField, axis as a vector, the internalField of centres and a single point, origin.

Have you tried:

U.internalField() = axis ˆ (centres.internalField() - origin);

Always think about what the fields mean.

Luarence, your insight was very useful, thank you. I've tryed that with other code but this code was buggy, then i didn't realized that i have to access the U data by this method.
The code of this post was a test code, i wrote the final version and it compiled, but in runtime i have some problems. This is the code:

// Velocity field initialization

// Angular velocity
scalar omega=10;

// Cell centroid coordinates
const volVectorField& centres = mesh.C();

// Coordinate system origin
point origin(0, 0, 0);

// Axis of rotation
vector axis(0, 0, 1);

// Force arm calculation
// Memory allocation
volVectorField palanca
(
IOobject
(
"pal",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh
);


palanca.internalField() = centres - ((centres.internalField()-origin) & axis)*axis;

// Velocity field calculation
U.internalField()=omega*mag(palanca.internalField( ))*axis ^ palanca.internalField();


in runtime i obtain this:

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Exec : myScalarTransportFoam
Date : Aug 18 2009
Time : 11:52:36
Host : elduque
PID : 6414
Case : /home/santiago/OpenFOAM/santiago-1.5/run/tuboMyTransport
nProcs : 1

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading field T

Reading field U

Reading transportProperties

Reading diffusivity D

Reading/calculating face flux field phi



NO_READ specified for read-constructor of object pal of class IOobject#0 Foam::error::printStack(Foam::Ostream&) in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Foam::regIOobject::readStream() in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#3 Foam::regIOobject::readStream(Foam::word const&) in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::GeometricField(Foam::IOobject const&, Foam::fvMesh const&) in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#5 main in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#6 __libc_start_main in "/lib/i686/cmov/libc.so.6"
#7 Foam::regIOobject::readIfModified() in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"


From function regIOobject::readStream(const word&)
in file db/regIOobject/regIOobjectRead.C at line 51.

FOAM aborting

Aborted

Really the only reason because i define the volVectorField palanca is to use it as a temporal variable. Because that i set NO_READ and NO_WRITE. It seems that this settings are incorrect. How I can define a temporal vector field to store the intermediate values (with the smallest useful structure)?. If volVectorField is the unique structure that i can use, how i can set it properly?.

Thanks.

You could create a temporary object:

tmp<volVectorField> variableName;

In your case your temporary variable is just a field of vectors, so you could just have a tmp<vectorField>. The vectorField is like a volVectorField but without the boundary conditions.
So something like:

tmp<vectorField> palanca = centres - ((centres.internalField()-origin) & axis)*axis;

Thanks Laurence!, slowly I'm understanding the structure of FOAM classes. tmp<vectorField> did the trick.

Now i could impose a rigid rotation velocity field with this lines:

tmp<vectorField> palanca;

palanca = centres.internalField() - ((centres.internalField()-origin) & axis)*axis;

// Velocity field calculation
U.internalField()=omega*axis ^ palanca;
U.write();

It compiled and ran properly, i verified the field using Paraview. Now i want to scale the field by the palanca's magnitude (the norm of arm of force vector), then i tryed the following:

U.internalField()=mag(palanca)*omega*axis ^ palanca;

It's good at compilation but don't run, this is the output:

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Exec : myScalarTransportFoam
Date : Aug 18 2009
Time : 19:38:28
Host : elduque
PID : 9150
Case : /home/santiago/OpenFOAM/santiago-1.5/run/tuboMyTransport
nProcs : 1

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading field T

Reading field U

Reading transportProperties

Reading diffusivity D

Reading/calculating face flux field phi



temporary deallocated#0 Foam::error::printStack(Foam::Ostream&) in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/santiago/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#3 Foam::tmp<Foam::Field<Foam::crossProduct<Foam::Vec tor<double>, Foam::Vector<double> >::type> > Foam::operator^<Foam::Vector<double>, Foam::Vector<double> >(Foam::tmp<Foam::Field<Foam::Vector<double> > > const&, Foam::tmp<Foam::Field<Foam::Vector<double> > > const&) in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#4 main in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"
#5 __libc_start_main in "/lib/i686/cmov/libc.so.6"
#6 Foam::regIOobject::readIfModified() in "/home/santiago/OpenFOAM/santiago-1.5/applications/bin/linuxGccDPOpt/myScalarTransportFoam"


From function const T& tmp<T>::operator()() const
in file /home/santiago/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/tmpI.H at line 190.

FOAM aborting

Aborted

I think something is wrong with the magnitude operator. Regards.

Maybe forget about the use of tmp, so just define palanca as a vectorField:

Someone can correct me, but I don't think you can take the mag of a tmp. You could try:

But I think for your purposes the first option would make life easier. Unless you're worried about memory usage.

Laurence, thank you again, your tips were very helpful, I'm imposing the velocity field as I wanted. Now, I'm not completely sure about freeing the temporal memory used, i tryed to invoke the vectorField destructor without success. What is the way in FOAM?. Regards.

I did a quick search on the forums about tmp, you should be able to find out more about using it. tmp variables can be cleared by using .clear(). You can avoid dereferencing them by putting brackets after the variable when you use it.

Thanks Laurance for all your support. I'm dealing now with the appropriate algebra to set the vector field, i want to compare the effects of high rotation velocities in advecition-diffusion solvers, specifically between Fluent and OpenFOAM. This is a work in colaboration with a team of chemical engineers with the aim of model polimerization reactors. We aren't completely sure about the influence of false diffusion in RTD determination (using the tracer method). Regards.

Sounds interesting, I do population balance modelling and turbulent, reactive, multiphase flow modelling. Numerical diffusion is something I've had to be careful about too.
Do you use population balances when modelling the reaction? Or are you just trying to look at the flow?
We've done a little bit of basic RTD analysis using Lagrangian particle tracking for a stirred tank in both experiments and simulations. The comparison wasn't great, and I think the main problems arose due to:

1) boundary conditions
2) inadequate turbulence modelling

I would "imagine" that numerical diffusion would have less effect on the mean residence time but a greater effect on the shape/width of the distribution?

Keep everyone posted on your progress.

Laurence, actually we aren't do chemical simulations really. We use a polimer-equivalent materials, most of them with a very high dynamic viscosity (like 1 Pa-seg). All of this lead us to use laminar-regime solvers. Once the velocity field is obtained we run a tracer to obtain the RTD.
In industry, because of the large time of simulation needed, all the people uses very long time steps, usually the RTD obtained looks like a Poisson distribution a lot.
My co-advisor suggested that behavior is due mainly to numerical diffusion in stead of real diffusion, in fact we use small diffusion factors in the solution of adv-diff equation. I'm suspecting now that this sounds more accurate in FEM solvers, because of the streamline diffusion methods (SUPG).
I'm finishing the coding now, when i have the results i'll post them here, so as to retain the topic.

Hi, Laurence & all the community, I've implemented the initialization for the U field. My new question is how change only one component of a field. Digging in the marvelous forums i found the following method, it compiled without warnings,

// Parabolic distribution in z component
scalarField zComp=Vmax*(1-sqr(mag(palanca)/tubeR));

// Velocity field calculation
U.internalField()=omega*(axis^palanca)+vector(0,0,1)*zComp;

basically it consists in using a mask, but... Is there a more "pure" way to do this, for example by a class' method?

Thanks in advance.

You may use Doxygen to find this out. If you search for scalarField in it, you'll find that it's a type definition: Field<type>. The following link gives the class reference for Field:

http://foam.sourceforge.net/doc/Doxy..._1_1Field.html

In the public member functions you'll find component. So to access one component of zComp would be (0=x axis, 1 = yaxis etc.):

Doxygen is very useful!

Yes, it's true, Doxygen is a powerful tool. I forgot to tell it, but I was doing research in Doxygen too. I've found the component method, but probably i used it in a wrong way, so your message push me to try again and now it works.
Actually I think that my previous message was bad written, i wanted to overwrite the z component of the U field (this was really the main topic). Anyway, the final code is the following:

// Velocity field calculation
//Solid rotation
U.internalField()=omega*(axis^palanca);

// Parabolic distribution in z component
U.internalField().component(2)=Vmax*(1-sqr(mag(palanca)/tubeR));

Regards.

Hi folks, as a continuation of my last post I'm copying all the code:

// Velocity field initialization

// Angular velocity
scalar omega=100;

// Maximum velocity en parabolic profile
scalar Vmax=2.0;

// Tube radius
scalar tubeR=0.1;

// Cell centroid coordinates
const volVectorField& centres = mesh.C();

// Coordinate system origin
point origin(0, 0, 0);

// Axis of rotation
vector axis(0, 0, 1);

// Force arm calculation
vectorField palanca;

palanca = centres.internalField() - ((centres.internalField()-origin) & axis)*axis;

// Velocity field calculation
//U.internalField()=omega*(axis^palanca);

// Parabolic distribution in z component
U.internalField().component(2)=Vmax*(1-sqr(mag(palanca)/tubeR));

U.write();

It's OK in compilation & running time but the results are bad. I obtain a zero field in the z component of U in stead of the parabolic field that i wanted. palanca has the distance of each cell center to the axis of rotation, then the equation for component 2 may give me a laminar parabolic profile in all the tube.
I don't know if the problem is in the mag or sqr methods or something else. All answers are welcome. Regards.

Some other important data:

1. I initialized U field to (0, 0, 0).
2. If we try to the following line:

U.internalField().component(2)=mag(palanca);

U field continue having (0, 0, 0)...

I guess that the component method only give us a copy of the data not a reference to overwrite the original.

Double check your axis vector. It should be the axis that your flow rotates around. Let's say that 'y' is up, 'x' is right and 'z' is out of the page. You want flow in the 'z' direction, so your axis vector should be your x or y axis. i.e. (0 1 0) or (1 0 0).
Does that make sense?

Edit: Think I misunderstood. Doesn't look like you can overwrite component. Try 'replace' (see Doxygen for GeometriField).

Maybe my no so pretty drawing can help us to understand. I want to impose a rotational field with a non-zero tangential velocity.

// Velocity field calculation
//U.internalField()=omega*(axis^palanca);

Due I impose (0,0,0) as the initial condition when I set the rotational velocity the axial velocity remains zero. Then I want to use:

// Parabolic distribution in z component
U.internalField().component(2)=Vmax*(1-sqr(mag(palanca)/tubeR));

to set the axial velocity. Component(2) implies setting the z component of the field, i.e. the axial one.
Actually the problem is that I've checked out with Paraview that this last line don't change the field. I guess that the component method only give us a copy of the component field not a reference, then the assignation is useless, but, I'm only guessing.
When I use this:

// Parabolic distribution in z component
scalarField zComp=Vmax*(1-sqr(mag(palanca)/tubeR));

// Velocity field calculation
U.internalField()=omega*(axis^palanca)+vector(0,0, 1)*zComp;things goes very well, trying to use component() was only to be more "pure" but it seems not to work. Will be another method or correct form of using component() to code this initialization a bit more elegant (and also to learn something more about OpenFOAM coding)?
Regards.

Yes, I see, try;

scalarField zComp=Vmax*(1-sqr(mag(palanca)/tubeR))
U.internalField().replace(Vector::Z,zComp);

I got replace from the Doxygen entry on GeometricField.

Laurence thanks for your answer, the code is running properly, now a detail the line:

must be

or

in order to run. The first line is calling directly the class Vector without giving all the information to the template, in this case the type of Vector that you want.
Finished this matter, and have another question, this one stopped my from your last post. When I run the case from a ./0/U with all the internalField set in 0 I can't obtain transport within the tube. But once I've ran the case and delete all the timesteps but the first one (0) and then run again things go well. I've checked that OpenFOAM writes the calculated (with our routine) U field in ./0/U field. I guess that FOAM uses the read data in stead of calculated or something like this, then things only run the second time.
Actually I'm not completely sure about what is the function of .write method (use of write is like in my first post, next of the initialization). Firstly I thought that it was to overwrite the values effectively, but now I think that this is the step when U field is being written in ./0 directory. But it remains in memory next to the .write calling?

Regards.