# finding average values for a plane in the vertical direction

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

 March 10, 2013, 08:55 finding average values for a plane in the vertical direction #1 Member   Nick Gutschow Join Date: Jan 2013 Posts: 36 Rep Power: 5 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

March 13, 2013, 19:46
#2
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,931
Rep Power: 41
Quote:
 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:
Code:
```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
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request

 March 15, 2013, 10:16 #3 Member   Nick Gutschow Join Date: Jan 2013 Posts: 36 Rep Power: 5 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 ....) in columns. Does anyone know if and how the sample utility could be used to do this? -Nick

 March 15, 2013, 11:18 #4 Member   Nick Gutschow Join Date: Jan 2013 Posts: 36 Rep Power: 5 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

March 15, 2013, 15:56
#5
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,931
Rep Power: 41
Quote:
 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)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request

 March 15, 2013, 16:02 #6 Member   Nick Gutschow Join Date: Jan 2013 Posts: 36 Rep Power: 5 Indeed. It is really just me shying away from Python for lack of experience. Thanks again for the input. -NG

March 18, 2013, 23:06
#7
Member

Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 7
Quote:
 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

March 19, 2013, 06:31
#8
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,931
Rep Power: 41
Quote:
 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:

Code:
```the Forces {
type patchExpression;
patches (
wall1
wall2
);
variables (
"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)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request

March 19, 2013, 10:43
#9
Member

Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 7
Quote:
 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: Code: ```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)
Thank you, Bernhard, this is exactly what I want.

For pressure force, the expression is simple,

Code:
`"forces= p*sf()*mag(sf()) ;"`
But for viscous force,

Code:
`"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

March 19, 2013, 11:16
#10
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,931
Rep Power: 41
Quote:
 Originally Posted by tfuwa Thank you, Bernhard, this is exactly what I want. For pressure force, the expression is simple, Code: `"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)

Quote:
 Originally Posted by tfuwa Code: `"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)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request

March 25, 2013, 03:44
#11
Member

Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 7
Quote:
 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,

The final expression is

Code:
`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

March 25, 2013, 06:55
#12
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,931
Rep Power: 41
Quote:
 Originally Posted by tfuwa Hi Bernhard, Many thanks for your help. The final expression is Code: `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

Quote:
 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)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request

March 25, 2013, 11:24
#13
Member

Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 7
Quote:
 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

March 25, 2013, 12:22
#14
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,931
Rep Power: 41
Quote:
 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
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request

 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 nishant_hull Open Source Meshers: Gmsh, Netgen, CGNS, ... 23 August 5, 2015 02:09 sunilpatil CFX 1 February 3, 2013 06:03 ouafa Open Source Meshers: Gmsh, Netgen, CGNS, ... 7 May 21, 2010 12:43 atulverma FLOW-3D 1 May 14, 2009 12:39 Dpu FLUENT 2 April 11, 2006 22:32

All times are GMT -4. The time now is 05:43.