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.

Sorry for the late reply, I've been really busy!
I'm not sure why you're getting the behaviour that you observe. Possibly U has gone out of scope or something? The write() function just writes the current U field to the 0/U file.
I'd probably need to see it for myself, feel free to post your files in a tarball and I can take a quick look.

Thanks Laurence for your attention even with all your work. Well, this is the code:

The sequence of includes is that follows:

(mine is the last one, added in scalarTransportFoam.C)

As you can see, after rewriting in memory the U vector field I call the .write() method to write the field to the hard disk. The problem is that if you run the code once, the initializated U field doesn't remain in memory, i.e. the advection equation is solved with another field, maybe zero, I don't
know. But, in the other hand, if you run the code some steps, then stop it and run it again with the /0/U recently written things go well.

That's all. Regards.

This worried me slightly, so I took a look, and I think I've figured it out.

The scalar transport equation takes phi as the coefficient in fvm::div, not U. So you have overwritten U, but now you have to recreate the flux.

So just put below Hope that's the problem.

Laurence, I've been working in another things last months but I've returned to this problem. Your suggestion was right, FOAM was assembling the flux with the original U not with the one calculated by my code. I had to put the

below

Nevertheless it is necessary to comment the original invocation of createPhi.H file in createFields.H because the

clause within createPhi.H. As you know, this line impedes that definition can be repeated.

Regards.

Dear Santiago,

I am struggling with getting a parabolic internal field for pipe flow in openfoam and I've found your post very useful.
Could you please give me some hints on the line "U.internalField()=omega*(axis^palanca);" ?
Another silly question, where shall I put the code, in 0/U or somewhere else? Do I need to derive my own solver for it?
I really appreciate your reply. Thank you very much.

Best regards,
Tony

Greetings to all!

FYI for future readers: Tony's question has been in essence answered here: http://www.cfd-online.com/Forums/ope...tml#post468380


@Tony: Regarding your question about Santiago's implementation: from what I can figure out from the first few posts, it seems that he created a variant of the solver scalarTransportFoam and added a new header file that includes the code discussed in this thread.

It takes a while to explain it in more detail, so I suggest that you have a look at the following tutorials, in order to get a bit more perspective into the details:

• http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam
• http://openfoamwiki.net/index.php/ScalarTransportFoam

In fact, you'll see that those two tutorials are similar to each other, since the first one implements the explanation given in the second one.


Best regards,
Bruno

Hi, Bruno

Thank you very much for the information. I will have a look at that and get back to you if I have further questions.

Best regards,
Tony