[NJG]

Hi Folks,

I have a finished run that looks good, paraView is showing expected behavior, but now I wish to create some plots for my data that paraView won't do.

Specifically, I want to take some scalar, let's say T, and average it on a given horizontal plane. So let's say z is vertical and x and y go from 0 to 10; I want to average the T values for all 100 volumes (spread through my x-y plane) at a given z value. Then, I want to end up with results that I can use to plot z vs T in a conventional Cartesian manner.

When I investigate my T file for some time step however, there is just a huge list of values. A given value is not explicitly tagged with a given x,y,z value in that file.

So, my questions are as follows I guess:

1) Is anyone aware of good way to accomplish my go to produce data for my z vs T plot?

2) If not, looking ahead to making my own function, where is it referenced what index in a given scalar file at a given timestep corresponds to a given spatial location? Or, is there a trick to output or export OpenFOAM data in a more readable and manipulable format to my end?

Thanks everyone!

-NG

[gschaider]

Originally Posted by NJG

Hi Folks,

I have a finished run that looks good, paraView is showing expected behavior, but now I wish to create some plots for my data that paraView won't do.

Specifically, I want to take some scalar, let's say T, and average it on a given horizontal plane. So let's say z is vertical and x and y go from 0 to 10; I want to average the T values for all 100 volumes (spread through my x-y plane) at a given z value. Then, I want to end up with results that I can use to plot z vs T in a conventional Cartesian manner.

When I investigate my T file for some time step however, there is just a huge list of values. A given value is not explicitly tagged with a given x,y,z value in that file.

So, my questions are as follows I guess:

1) Is anyone aware of good way to accomplish my go to produce data for my z vs T plot?

2) If not, looking ahead to making my own function, where is it referenced what index in a given scalar file at a given timestep corresponds to a given spatial location? Or, is there a trick to output or export OpenFOAM data in a more readable and manipulable format to my end?

Thanks everyone!

-NG

Well you can have the cake and eat it to: do it in paraview although paraview doesn't know it. The thing is that Paraview ... if you have a Paraview with Python support compiled in ... and the version is not totally old .... comes with numpy a nice library for matrix operations. Instead of cutting planes you'd divide the range of z into N "bins" have a look which cell falls into which bin (there is a function in numpy for that) and calculate a weighted average. Approximately like this:

bins=numpy.linspace(min(minY),max(maxY),25)
vol,edges=numpy.histogram(heightField,bins=bins,weights=volume)
alpha,edges=numpy.histogram(heightField,bins=bins,weights=volume*alphaField)

xVals=0.5*(edges[1:]+edges[0:-1])
vals=alpha/vol

data=numpy.asarray([xVals,vals]).T
numpy.savetxt(dataFile("alphaByHeight"),data)

The things why you can't immediately use this is:
- This is for alpha1 not for T
- for height Y not Z
- You'll probably have to import numpy
- This was lifted from an example of swak4Foam (fillingTheDam) and the variables minY, maxY, heightField (which is Y), volume (the cell volumes), and alphaField reach the python-namespace through some swak-Stuff (this is from a function-object that evaluates this python-code at each write-time).

I must admit I haven't done much python in Paraview recently so I can't tell you how to best get to the data (with the "Python Calculator", "Programmable Filter" or by working on the Python Shell) but if you get to a solution I think it would be much appreciated if you shared it here

[NJG]

gschaider,

Thanks for the reply! I have got to say however, there MUST be an easier way. Actually, this is driving me a bit nuts now.

Using the sample utility I am able to produce a 1d line of data using the "raw format for sets as follows:

1 1 1 300
2 1 1 301
3 1 1 302

Where that is x, y, z, and T for temperature say. But I can only vary one of the x, y or, z values I guess? I want an output just like this, but I want it to just pump out the T value for every cell, not just the cells in a single line.

I find it hard to believe that this functionality does not exist in this utility. It seems so basic... The most basic of outputs; just a text file giving (x y z <field 1> <field 2> <field 3>....) in columns.

Does anyone know if and how the sample utility could be used to do this?

-Nick

[NJG]

Ok, I think I found my own work around.

From paraView I can just save my data as a csv, which outputs all my fields and points.

From here I am half way through a MATLAB code to basically area weight average the scalar in question for a given plane and then set these to vertical displacement values equal to the center point of that plane.

-NG

[gschaider]

Originally Posted by NJG

Ok, I think I found my own work around.

From paraView I can just save my data as a csv, which outputs all my fields and points.

From here I am half way through a MATLAB code to basically area weight average the scalar in question for a given plane and then set these to vertical displacement values equal to the center point of that plane.

Which basically boils down to my solution (or are you writing out the plane data?): binning all cells at the same height together.

That is the nice thing about people: that they have different views what is "easy" and each of these ideas is valid: for you the idea of clicking through two different programs and letting them "communicate" through files, for me to invest some work in the beginning and then just starting the run and have all the data ready for me when it is finished (of course your way is more time efficient if you only have to do this once for one run)

[NJG]

Indeed. It is really just me shying away from Python for lack of experience. Thanks again for the input.

-NG

[tfuwa]

Originally Posted by gschaider

Which basically boils down to my solution (or are you writing out the plane data?): binning all cells at the same height together.

That is the nice thing about people: that they have different views what is "easy" and each of these ideas is valid: for you the idea of clicking through two different programs and letting them "communicate" through files, for me to invest some work in the beginning and then just starting the run and have all the data ready for me when it is finished (of course your way is more time efficient if you only have to do this once for one run)

Hi there,

I would like to choose part of a patch (for instance, cylinder faces in a specific domain), and then calculated the forces on this selected part. Is there a way to do this in OpenFoam at run-time? Thanks a lot, Albert

[gschaider]

Originally Posted by tfuwa

Hi there,

I would like to choose part of a patch (for instance, cylinder faces in a specific domain), and then calculated the forces on this selected part. Is there a way to do this in OpenFoam at run-time? Thanks a lot, Albert

I'm not clear on your definition of "cylinder faces" and "calculate forces".

If you're going for complete patches you might look into the forces-functionObject that comes with OF. If you only want to select parts of the patches then you can use for instance a swak-patch expression that looks like this:

the Forces {
  type patchExpression;
  patches (
    wall1
    wall2 
  );
  variables (
     "radius=0.1;"
     "forces= .... ;" // that is the part you have to specify
  );
  expression "(sqrt(pos().x*pos().x+pos().y*pos().y)<radius) ? forces : vector(0,0,0)";
  accumulations (
     sum
  );
}

that would sum up the force for all faces on the patches which are within 0.1 to the z-axis. Specification of force is you responsibility but you can use any field that is available and lots of operators (see the groovyBC-documentation)

[tfuwa]

Originally Posted by gschaider

I'm not clear on your definition of "cylinder faces" and "calculate forces".

If you're going for complete patches you might look into the forces-functionObject that comes with OF. If you only want to select parts of the patches then you can use for instance a swak-patch expression that looks like this:

the Forces {
  type patchExpression;
  patches (
    wall1
    wall2 
  );
  variables (
     "radius=0.1;"
     "forces= .... ;" // that is the part you have to specify
  );
  expression "(sqrt(pos().x*pos().x+pos().y*pos().y)<radius) ? forces : vector(0,0,0)";
  accumulations (
     sum
  );
}

that would sum up the force for all faces on the patches which are within 0.1 to the z-axis. Specification of force is you responsibility but you can use any field that is available and lots of operators (see the groovyBC-documentation)

Thank you, Bernhard, this is exactly what I want.

For pressure force, the expression is simple,

"forces= p*sf()*mag(sf()) ;"

But for viscous force,

"forces= rho()*nu*dev(twoSymm(grad(U))) *sf()*mag(sf());"

It seems to me that dev, grad, curl and laplacian are not available in groovyBC. How to apply those operators? Many thanks. Also, please correct me if there are any errors in the force expressions.

Regards,
Albert

[gschaider]

Originally Posted by tfuwa

Thank you, Bernhard, this is exactly what I want.

For pressure force, the expression is simple,

"forces= p*sf()*mag(sf()) ;"

But for viscous force,

mag(sf()) is called area in OF. So alternative might be "p*normal()*area()" (not sure about the sign)

Originally Posted by tfuwa

"forces= rho()*nu*dev(twoSymm(grad(U))) *sf()*mag(sf());"

It seems to me that dev, grad, curl and laplacian are not available in groovyBC. How to apply those operators? Many thanks. Also, please correct me if there are any errors in the force expressions.

The reason is that to calculate these you would have to calculate them in the whole field and I don't want a lowly BC that "only" works on a patch to trigger that.

For your expression: would have to look it up myself how the original functionObjects calculate that. But what might be sufficient to reformulate it in terms of snGrad(U) (that is available in groovyBC). Otherwise you'd have to use an expressionField-functionObject that creates a field named devGradU with the expression "dev(twoSymm(grad(U)))" and use that field in your expression (not sure but you might have to use internalField(devGradU) if that field is 0 on the boundary).

Also: rho is a field. So no (). nu depends on the solver whether this is available as a field (only then swak can access it. There are ways to get it anyway)

[tfuwa]

Originally Posted by gschaider

mag(sf()) is called area in OF. So alternative might be "p*normal()*area()" (not sure about the sign)


The reason is that to calculate these you would have to calculate them in the whole field and I don't want a lowly BC that "only" works on a patch to trigger that.

For your expression: would have to look it up myself how the original functionObjects calculate that. But what might be sufficient to reformulate it in terms of snGrad(U) (that is available in groovyBC). Otherwise you'd have to use an expressionField-functionObject that creates a field named devGradU with the expression "dev(twoSymm(grad(U)))" and use that field in your expression (not sure but you might have to use internalField(devGradU) if that field is 0 on the boundary).

Also: rho is a field. So no (). nu depends on the solver whether this is available as a field (only then swak can access it. There are ways to get it anyway)

Hi Bernhard,

Many thanks for your help.

The final expression is

p*rho*normal()*area()-rho*nu*snGrad(U)*area()

where the first term is due to pressure and the second represents the contribution of viscous.

BTW, what is the difference between sf() and area() ?

Regards,
Albert

[gschaider]

Originally Posted by tfuwa

Hi Bernhard,

Many thanks for your help.

The final expression is

p*rho*normal()*area()-rho*nu*snGrad(U)*area()

where the first term is due to pressure and the second represents the contribution of viscous.

Cool. Even cooler if you could add a small recipe to
http://openfoamwiki.net/index.php/Co...Usage_examples

Originally Posted by tfuwa

BTW, what is the difference between sf() and area() ?

Basically "Sf()" is "area()*normal()" (for OpenFOAM that describes all it needs to know about that face in a single value)

[tfuwa]

Originally Posted by gschaider

Cool. Even cooler if you could add a small recipe to
http://openfoamwiki.net/index.php/Co...Usage_examples



Basically "Sf()" is "area()*normal()" (for OpenFOAM that describes all it needs to know about that face in a single value)

Hi Bernhard,

Thanks. Just let you know the wiki-page has been updated to include this usage.

Kind regards,
Albert

[gschaider]

Originally Posted by tfuwa

Hi Bernhard,

Thanks. Just let you know the wiki-page has been updated to include this usage.

Just saw it. Great