Tracking position of maximum vorticity

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

 July 4, 2009, 20:20 Tracking position of maximum vorticity #1 Member   Pascal Join Date: Jun 2009 Location: Montreal Posts: 65 Rep Power: 8 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

 July 6, 2009, 04:52 #2 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 9 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

 July 7, 2009, 08:45 #4 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 9 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

 July 7, 2009, 12:14 #5 Member   Pascal Join Date: Jun 2009 Location: Montreal Posts: 65 Rep Power: 8 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.

 July 9, 2009, 00:11 #6 Member   Pascal Join Date: Jun 2009 Location: Montreal Posts: 65 Rep Power: 8 Any suggestion for tracking that local maximum? Thanks

 July 9, 2009, 18:59 #7 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 9 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

 July 11, 2009, 17:54 #8 Member   Pascal Join Date: Jun 2009 Location: Montreal Posts: 65 Rep Power: 8 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

July 12, 2009, 13:03
#9
Senior Member

Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
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

Henrik

 July 13, 2009, 01:54 #10 Member   Pascal Join Date: Jun 2009 Location: Montreal Posts: 65 Rep Power: 8 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

 July 13, 2009, 14:36 #11 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 9 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

 July 13, 2009, 15:26 #12 Member   Pascal Join Date: Jun 2009 Location: Montreal Posts: 65 Rep Power: 8 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

 July 13, 2009, 16:05 #13 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 9 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

 May 25, 2013, 11:48 Tracking Vorticity #14 New Member   Luis Batista Join Date: Mar 2013 Location: Lisboa / Setúbal Posts: 17 Rep Power: 4 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

 May 27, 2013, 09:58 #15 Member   Pascal Join Date: Jun 2009 Location: Montreal Posts: 65 Rep Power: 8 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

 Tags tracking vorticity

 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 Monika FLUENT 1 September 4, 2014 10:30 mvoss CFX 3 June 23, 2009 08:07 dm2747 FLUENT 0 April 17, 2009 01:29 pxyz Main CFD Forum 37 July 7, 2006 08:42 Art Stretton Phoenics 5 April 2, 2002 05:59

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