|
[Sponsors] |
July 4, 2009, 21:20 |
Tracking position of maximum vorticity
|
#1 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
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, 05:52 |
|
#2 |
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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 6, 2009, 18:17 |
|
#3 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
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 |
|
July 7, 2009, 09:45 |
|
#4 |
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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; |
|
July 7, 2009, 13:14 |
|
#5 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
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 17:29. |
|
July 9, 2009, 01:11 |
|
#6 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
Any suggestion for tracking that local maximum?
Thanks |
|
July 9, 2009, 19:59 |
|
#7 |
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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; Henrik |
|
July 11, 2009, 18:54 |
|
#8 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
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; Regards, Pascal |
|
July 12, 2009, 14:03 |
|
#9 | ||
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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:
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:
http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam Then add your code. Henrik |
|||
July 13, 2009, 02:54 |
|
#10 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
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, 15:36 |
|
#11 |
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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, 16:26 |
|
#12 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
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; } } // ************************************************************************* // Pascal |
|
July 13, 2009, 17:05 |
|
#13 |
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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! } Henrik |
|
May 25, 2013, 12:48 |
Tracking Vorticity
|
#14 |
New Member
Luis Batista
Join Date: Mar 2013
Location: Lisboa / Setúbal
Posts: 17
Rep Power: 13 |
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, 10:58 |
|
#15 |
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 17 |
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 |
|
February 1, 2017, 19:09 |
|
#16 | |
Member
|
Quote:
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 |
||
Tags |
tracking vorticity |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Position Tracking | mvoss | CFX | 4 | May 19, 2019 11:08 |
Tracking droplet position | Monika | FLUENT | 1 | September 4, 2014 11:30 |
DPM UDF particle position using the macro P_POS(p)[i] | dm2747 | FLUENT | 0 | April 17, 2009 02:29 |
transform navier-stokes eq. to euler-eq. | pxyz | Main CFD Forum | 37 | July 7, 2006 09:42 |
Combustion Convergence problems | Art Stretton | Phoenics | 5 | April 2, 2002 06:59 |