CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Tracking position of maximum vorticity

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

Reply
 
LinkBack Thread Tools Display Modes
Old   July 4, 2009, 20:20
Default Tracking position of maximum vorticity
  #1
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hi all,

I'm simulating vortex in ground effect using OpenFOAM 1-5 and I would like to know how I could track the position (x,y) of the maximum vorticity for each time step.

I try to modify the vorticity utility but no success.

Thanks a lot,

Pascal
Pascal_doran is offline   Reply With Quote

Old   July 6, 2009, 04:52
Default
  #2
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Dear Pascal,

My recipie:

1) Modify vorticity utility to report cell with maximum vorticity
2) Report position in real coordinates
3) Add code to solver
4) Write function object

Guess you are stuck on 1. Please post your modifications so that we get in position where we can help. An error message would be nice, too.

Henrik
henrik is offline   Reply With Quote

Old   July 6, 2009, 17:17
Default
  #3
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Thanks for your reply,

All this OpenFOAM programming is new for me... (It is the first time a use OpenFOAM)
As you said I'm stuck at step 1.

Here is what I tried to do in the vorticity utility (I renamed it trace_max_vor.C) for writing the x and y position of the maximum vorticity :

#include "calc.H"
#include "fvc.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{
bool writeResults = !args.options().found("noWrite");

IOobject Uheader
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ
);

if (Uheader.headerOk())
{
Info<< " Reading U" << endl;
volVectorField U(Uheader, mesh);

Info<< " Calculating vorticity" << endl;
volVectorField vorticity
(
IOobject
(
"vorticity",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
fvc::curl(U)
);

volScalarField magVorticity
(
IOobject
(
"magVorticity",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
mag(vorticity)
);


volScalarField pos_x_Vorticity
(
IOobject
(
"pos_x_Vorticity",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
max(magVorticity).value().component(0)
);

volScalarField pos_y_Vorticity
(
IOobject
(
"pos_y_Vorticity",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
max(magVorticity).value().component(1)
);

Info<< "vorticity max/min : "
<< max(magVorticity).value() << " "
<< min(magVorticity).value() << endl;

if (writeResults)
{
vorticity.write();
magVorticity.write();
pos_x_Vorticity.write();
pos_x_Vorticity.write();
}
}
else
{
Info<< " No U" << endl;
}
}

// ************************************************** *********************** //

And here is the error message :

SOURCE=trace_max_vor.C ; g++ -m64 -Dlinux64 -DDP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -march=opteron -O3 -DNoRepository -ftemplate-depth-40 -I/home/pdoran/OpenFOAM/OpenFOAM-1.5/src/postProcessing/postCalc -I/home/pdoran/OpenFOAM/OpenFOAM-1.5/src/finiteVolume/lnInclude -IlnInclude -I. -I/home/pdoran/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude -I/home/pdoran/OpenFOAM/OpenFOAM-1.5/src/OSspecific/Unix/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/trace_max_vor.o
trace_max_vor.C: In function ‘void Foam::calc(const Foam::argList&, const Foam::Time&, const Foam::fvMesh&)’:
trace_max_vor.C:92: error: request for member ‘component’ in ‘Foam::max(const Foam::GeometricField<Type, PatchField, GeoMesh>&) [with Type = double, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]().Foam::dimensioned<Type>::value [with Type = double]()’, which is of non-class type ‘double’
trace_max_vor.C:104: error: request for member ‘component’ in ‘Foam::max(const Foam::GeometricField<Type, PatchField, GeoMesh>&) [with Type = double, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]().Foam::dimensioned<Type>::value [with Type = double]()’, which is of non-class type ‘double’
make: *** [Make/linux64GccDPOpt/trace_max_vor.o] Error 1
Pascal_doran is offline   Reply With Quote

Old   July 7, 2009, 08:45
Default
  #4
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Dear Pascal,

you cannot use field algebra to determine the position where the maximum vorticity is. Try something like this:

Code:
scalar maxVorticity = -GREAT;
label maxCellI;

forAll(magVorticity, cellI)
{
    if ( maxVorticity < magVorticity[cellI] )
    {
        maxVorticity = magVorticity[cellI];
        maxCellI = cellI;
    }
}
Info << "max @ cell : " << maxCellI << endl;
Info << "pos @ max : " << mesh.C()[maxCellI] << endl;
Henrik
henrik is offline   Reply With Quote

Old   July 7, 2009, 12:14
Default
  #5
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Dear Dr Rusche,

I really appreciate your very good advice. It works perfectly!
Thank you very much

But now I would like to create a new utility that will track a local maximum. The previous utility was used for tracking the maximum vorticity (max. in all the field).

What I would like to do is:

first: Find the max vorticity in all the field at time 0.0 (which is already done thanks to you)
second : For the next time step (0.1, 0.2, 0.3, ..., n), find the nearest maximum vorticity from the previous max. at time n-1.

And this is the issue, I'm not able to write the appropriate line of code to say: based on the previous position of maximum vorticity, find the nearest max. vorticity.

Here's what I Tried:

scalar maxVorticity = -GREAT;
label maxCellI;
scalar distance = 0.0;
scalar prev_distance = 0.0;
scalar pos_x = 0.0;
scalar pos_y = 0.0;
scalar prev_pos_x = 0.0;
scalar prev_pos_y = 0.0;

forAll(magVorticity, cellI)
{
pos_x = mesh.C()[cellI].component(0);
pos_y = mesh.C()[cellI].component(1);
distance = :ow((pos_x-prev_pos_x)*(pos_x-prev_pos_x) + (pos_y-prev_pos_y)*(pos_y-prev_pos_y),0.5);

if ((runTime.timeName() == 0 || distance < 0.5) && maxVorticity < magVorticity[cellI])
//if (maxVorticity < magVorticity[cellI])
{
maxVorticity = magVorticity[cellI];
maxCellI = cellI;

prev_pos_x = pos_x;
prev_pos_y = pos_y;
prev_distance = distance;
}
}
Info<< "X,Y " << mesh.C()[maxCellI].component(0) << " "<< mesh.C()[maxCellI].component(1) << endl;
Info<< "Dist. from prev. max. : " << prev_distance << endl;
Info<< "Local max. vorticity : "<< maxVorticity << endl;

I have no error message, it just gives wrong x,y positions. I think the problem is that I overwrite prev_pos_x and prev_pos_y for every time increment, but I don't know how fixed that since the time loop is not in that file.

Thank you for your help,

Pascal

Last edited by Pascal_doran; July 7, 2009 at 16:29.
Pascal_doran is offline   Reply With Quote

Old   July 9, 2009, 00:11
Default
  #6
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Any suggestion for tracking that local maximum?
Thanks
Pascal_doran is offline   Reply With Quote

Old   July 9, 2009, 18:59
Default
  #7
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Dear Pascal,

the biggest problem I can spot is

runTime.timeName() == 0

First of all, have a look at $FOAM_SRC/OpenFOAM/lnInclude/Time.H

runTime.timeName() is a word!

I think you want

runTime.value() < SMALL

and not even that will not work if startTime is not zero. But let's not worry about this right now.

Some more words about programming in OF's C++:

Code:
scalar maxVorticity = -GREAT;
label maxCellI;
scalar distance = 0.0;
scalar prev_distance = 0.0;
vector pos (0,0,0);
vector prev_pos (0,0,0);
    
forAll(magVorticity, cellI)
{
     distance = mag(mesh.C()[cellI]-prev_pos);

     if ((runTime.timeName() == 0 || distance < 0.5) && maxVorticity < magVorticity[cellI])
     //if (maxVorticity < magVorticity[cellI])
     {
            maxVorticity = magVorticity[cellI];
            maxCellI = cellI;

            prev_pos = mesh.C()[cellI];
            prev_distance = distance;
    }
}
Info<< "pos " << mesh.C()[maxCellI] << endl;
Info<< "Dist. from prev. max. : " << prev_distance << endl;
Info<< "Local max. vorticity : "<< maxVorticity << endl;
Get my drift?

Henrik
henrik is offline   Reply With Quote

Old   July 11, 2009, 17:54
Question
  #8
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Thank you for advise Dr. Rusche,

But I'm still overwriting the vector prev_pos for each loop of time. Should I change something in the application file?

here's what the code look alike:
Code:
    scalar maxVorticity = -GREAT;
    label maxCellI;
    scalar distance = 0.0;
    scalar prev_distance = 0.0;
    vector prev_pos (0,0,0);

    forAll(magVorticity, cellI)
    {
        distance = mag(mesh.C()[cellI]-prev_pos);

           if ((runTime.value() < SMALL || distance < 0.5) && (maxVorticity < magVorticity[cellI]))
        {
            //Info<< "runTime.value() "<< runTime.value() << endl;
            Info<< "mesh.C()[cellI] "<< mesh.C()[cellI] << endl; //This output show me that I overwrite prev_pos every time loop.
            //Info<< "distance "<< distance << endl;
            //Info<< "maxVorticity "<< maxVorticity << endl;
            //Info<< "magVorticity[cellI] "<< magVorticity[cellI] << endl;
            //Info<< "---------------------------------------------------"<< endl;

            maxVorticity = magVorticity[cellI];
            maxCellI = cellI;

            prev_pos = mesh.C()[cellI];
            prev_distance = distance;
            }
    }
    Info<< "pos " << mesh.C()[maxCellI] << endl;
I have no error message, just the wrong position of the local maximum vorticity.
Regards,

Pascal
Pascal_doran is offline   Reply With Quote

Old   July 12, 2009, 13:03
Default
  #9
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Dear Pascal,

obviously there is an error in the logic of your code. Just do what everybody has to do: Add Info statements, analyse where it goes wrong and fix it - This is too hard to do remotely.

Quote:
But I'm still overwriting the vector prev_pos for each loop of time.
Why do you think this is a problem?

One thing I can spot is that you overwrite prev_pos within the loop over all cells. This has an unwanted effect on the distance calculation!

I think, prev_pos should be reset after the loop over all cell and exactly once per time step.

Quote:
Should I change something in the application file?
Yes, I would hardcode it into the application because then you do not to rely on the dump frequency. If you don't know how to change an existing application - Have a look at step 2 of this "How to"

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

Then add your code.

Henrik
henrik is offline   Reply With Quote

Old   July 13, 2009, 01:54
Default
  #10
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Dr. Rusche,

I already added info statement analyze where it goes wrong (it is the uncommented info statement, see my previous post).

Firstly, I know that this code, how it is written right now, give me the right maximum vorticity at initial time (which is 0,0 s) and prev_pos is set to the right location of that maximum vorticity, let say prev_pos = (x1,y1,z1) (after the loop over all cells). But, at the next time loop (which is 0,1 s), the prev_pos is re-initialized to (0,0,0) which is incorrect. I would like that prev_pos keep its last value which should be (x1,y1,z1). I don’t know where the time loop in that utility is and I would like to find it.

Secondly, I want to keep the previous value of prev_pos for the next time step because I need that information to find the closest local maximum vorticity of the previous one. Remember that local maximum is moving, but I know that it will be in a distance less than 0,5 from the previous local maximum vorticity.


That is the reason why I need to know the previous and the actual position location of the local maximum vorticity. I want to track the location of the same vortex for all the duration of the simulation.

Thank you, I really appreciate your help.

Pascal
Pascal_doran is offline   Reply With Quote

Old   July 13, 2009, 14:36
Default
  #11
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Dear Pascal,

it's difficult to help because you did not post the full source code including the time or pseudo time loop which reads the fields.

My feeling is that since, prev_pos keeps track of the "old" vortex position, it should be initalised BEFORE the time loop and updated with the new position AT THE END of the time loop AFTER the loop over all cells.

Henrik
henrik is offline   Reply With Quote

Old   July 13, 2009, 15:26
Default
  #12
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Here is the entire code for the local maximum vorticity utility (based on vorticity utility):

HTML Code:
\*---------------------------------------------------------------------------*/
#include "calc.H"
#include "fvc.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{
bool writeResults = !args.options().found("noWrite");
IOobject Uheader
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ
);
if (Uheader.headerOk())
{
Info<< " Reading U" << endl;
volVectorField U(Uheader, mesh);
Info<< " Calculating vorticity" << endl;
volVectorField vorticity
(
IOobject
(
"vorticity",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
fvc::curl(U)
);
volScalarField magVorticity
(
IOobject
(
"magVorticity",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
mag(vorticity)
);
scalar maxVorticity = -GREAT;
label maxCellI;
scalar distance = 0.0;
scalar prev_distance = 0.0;
scalar pos_x = 0.0;
scalar pos_y = 0.0;
scalar prev_pos_x = 0.0;
scalar prev_pos_y = 0.0;
forAll(magVorticity, cellI)
{
pos_x = mesh.C()[cellI].component(0);
pos_y = mesh.C()[cellI].component(1);
distance = ::pow((pos_x-prev_pos_x)*(pos_x-prev_pos_x) + (pos_y-prev_pos_y)*(pos_y-prev_pos_y),0.5);
//if ((runTime.timeName() == 0 || distance < 0.5) && (maxVorticity < magVorticity[cellI]))
if (maxVorticity < magVorticity[cellI])
{
maxVorticity = magVorticity[cellI];
maxCellI = cellI;
prev_pos_x = pos_x;
prev_pos_y = pos_y;
prev_distance = distance;
}
}
Info<< "X,Y " << mesh.C()[maxCellI].component(0) << " "<< mesh.C()[maxCellI].component(1) << endl;
// Info<< "Dist. from prev. max. : " << prev_distance << endl;
// Info<< "Local max. vorticity : "<< maxVorticity << endl;
if (writeResults)
{
// vorticity.write();
// magVorticity.write();
}
}
else
{
Info<< " No U" << endl;
}
}
// ************************************************************************* //
Thank you,

Pascal
Pascal_doran is offline   Reply With Quote

Old   July 13, 2009, 16:05
Default
  #13
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Dear Pascal,

now I understand. Haven't looked at the source of the utilities for a while - so it didn't click.

Your problem is that the utility is actually a global function Foam::calc which linked against a generic main compiled into postCalc.o. The consequence of this is that you where missing the "time loop" (Yes, you said and I did't listen carefully) and prev_pos get's deleted at the end of the function.

The easiest (and pretty ugly) option is to define prev_pos outside of the function which makes it a global variable.

Code:
Foam::vector prev_pos(0,0,0);

void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{
// Do your stuff
prev_pos = mesh.C()[cellI]; // No preceding class name (vector) here!
}
Make sure not to instantiate another prev_pos in the function!

Henrik
henrik is offline   Reply With Quote

Old   May 25, 2013, 11:48
Default Tracking Vorticity
  #14
New Member
 
Luis Batista
Join Date: Mar 2013
Location: Lisboa / Setúbal
Posts: 17
Rep Power: 4
Luis Batista is on a distinguished road
Hello,

I understand this is quite an old post, however, I would like to know how did you manage to sort out this issue.

Can you please detail?

regards,
Luis
Luis Batista is offline   Reply With Quote

Old   May 27, 2013, 09:58
Default
  #15
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hi Luis,

Take a look here and you will learn how to code the time loop as a post-process utility :

/home/.../OpenFOAM/OpenFOAM-2.1.0/applications/utilities/postProcessing/scalarField/pPrime2

Regards,
Pascal
Pascal_doran is offline   Reply With Quote

Reply

Tags
tracking vorticity

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Tracking droplet position Monika FLUENT 1 September 4, 2014 10:30
Position Tracking mvoss CFX 3 June 23, 2009 08:07
DPM UDF particle position using the macro P_POS(p)[i] dm2747 FLUENT 0 April 17, 2009 01:29
transform navier-stokes eq. to euler-eq. pxyz Main CFD Forum 37 July 7, 2006 08:42
Combustion Convergence problems Art Stretton Phoenics 5 April 2, 2002 05:59


All times are GMT -4. The time now is 14:16.