[Pascal_doran]

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

[henrik]

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

[Pascal_doran]

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

[henrik]

Dear Pascal,

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

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

[Pascal_doran]

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

[Pascal_doran]

Any suggestion for tracking that local maximum?
Thanks

[henrik]

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++:

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

[Pascal_doran]

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:

    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

[henrik]

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.

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.

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

[Pascal_doran]

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

[henrik]

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

[Pascal_doran]

Here is the entire code for the local maximum vorticity utility (based on vorticity utility):

\*---------------------------------------------------------------------------*/
#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

[henrik]

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.

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

[Luis Batista]

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

[Pascal_doran]

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

[reza_65]

Originally Posted by henrik

Dear Pascal,

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

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

Hi All,

I have done the same sort of tracking method for finding the Flame tip position, and everything is working fine, but the only issue, is the precision of this value (maximum X value), which it has just 4 digit precision!
I increased writing precision in controldict to 12, everything (timesteps, cournat num,..) are having 12 digit precision in the log file, except this flame position (which is just 4 digit like: 1.3013 m).
Do you have any idea how can I increase the writing precision of this scalar?

Thank you advanced.
Reza