CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Community Contributions (https://www.cfd-online.com/Forums/openfoam-community-contributions/)
-   -   [swak4Foam] Track vortices using Swak4FOAM (https://www.cfd-online.com/Forums/openfoam-community-contributions/99782-track-vortices-using-swak4foam.html)

vinz April 11, 2012 11:16
Track vortices using Swak4FOAM
 
Hi all,

I am simulating the propagation in space and time of some vortices.
I would like to be able to trace the center of these vortices into my field and then be able to compute the circulation on a closed surface around this center.

So, first step I need to trace the center of my vortex. It looks like defining some functions in the controlDict and using the libsimpleSwakFunctionObjects could be an option but I am not sure how to do it.
To start easy I was planning to find the center of my vortex by looking for the maximum vorticity in my field.

I wrote those lines in my dictionary:
Code:
functions
(
    defineFollowVortex {
        type addGlobalVariable;
        outputControl timeStep;
        outputInterval 1;

        globalScope followVortex;
        globalVariables {
            yThres {
                valueType scalar;
                value 0.2;
            }
        }
    };
    createVortexCenter {
        type createSampledSet;
        outputControl timeStep;
        outputInterval 1;
        setName vortexCenter;
        set {
            type cloud;
            axis y;
            points( (0.5 2 0.0) );
        }
    }

    vorticityMax {
        type swakExpression;
        verbose true;
        valueType set;
        setName vorticityMax;
        expression "(pos.y()>yThres) ? magVorticity : 0";
        accumulations ( max );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }

    vortexHeight {
        type swakExpression;
        verbose true;
        valueType set;
        setName  fillingHeight;
        globalScopes ( "fillDam" );
        expression "(magVorticity >= 0.99*vorticityMax) ? pos().y : 0";
        accumulations (
            max
        );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }
);
However I am really not sure it is the way to go and I get the following error message:
Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

[0]
[0]
[0] --> FOAM FATAL IO ERROR:
[0] "ill defined primitiveEntry starting at keyword 'functions' on line 60 and ending at line 124"
[0]
[0] file: /CFDBB/vincent/OpenFOAM/vincent-2.1.0/run/tutorials/compressible/rhoPimpleDyMFoamWithSources/LES/AfterMeeting120330/testVortexHeight/processor0/../system/controlDict at line 124.
[0]
[0]    From function primitiveEntry::readEntry(const dictionary&, Istream&)
[0]    in file lnInclude/IOerror.C at line 132.
[0]
FOAM parallel run exiting
Any idea?

gschaider April 11, 2012 15:53
Quote:
Originally Posted by vinz (Post 354201)
Hi all,

I am simulating the propagation in space and time of some vortices.
I would like to be able to trace the center of these vortices into my field and then be able to compute the circulation on a closed surface around this center.

So, first step I need to trace the center of my vortex. It looks like defining some functions in the controlDict and using the libsimpleSwakFunctionObjects could be an option but I am not sure how to do it.
To start easy I was planning to find the center of my vortex by looking for the maximum vorticity in my field.

I wrote those lines in my dictionary:
Code:
functions
(
    defineFollowVortex {
        type addGlobalVariable;
        outputControl timeStep;
        outputInterval 1;

        globalScope followVortex;
        globalVariables {
            yThres {
                valueType scalar;
                value 0.2;
            }
        }
    };
    createVortexCenter {
        type createSampledSet;
        outputControl timeStep;
        outputInterval 1;
        setName vortexCenter;
        set {
            type cloud;
            axis y;
            points( (0.5 2 0.0) );
        }
    }

    vorticityMax {
        type swakExpression;
        verbose true;
        valueType set;
        setName vorticityMax;
        expression "(pos.y()>yThres) ? magVorticity : 0";
        accumulations ( max );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }

    vortexHeight {
        type swakExpression;
        verbose true;
        valueType set;
        setName  fillingHeight;
        globalScopes ( "fillDam" );
        expression "(magVorticity >= 0.99*vorticityMax) ? pos().y : 0";
        accumulations (
            max
        );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }
);
However I am really not sure it is the way to go and I get the following error message:
Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

[0]
[0]
[0] --> FOAM FATAL IO ERROR:
[0] "ill defined primitiveEntry starting at keyword 'functions' on line 60 and ending at line 124"
[0]
[0] file: /CFDBB/vincent/OpenFOAM/vincent-2.1.0/run/tutorials/compressible/rhoPimpleDyMFoamWithSources/LES/AfterMeeting120330/testVortexHeight/processor0/../system/controlDict at line 124.
[0]
[0]    From function primitiveEntry::readEntry(const dictionary&, Istream&)
[0]    in file lnInclude/IOerror.C at line 132.
[0]
FOAM parallel run exiting
Any idea?
I'm not sure about the syntax error. Maybe ; following defineFollowVortex.

Anyway. What might help you is minPosition/maxPosition (this is currently only available in the development-version of swak. You'll have to get that via Mercurial): they take an expression and return the position where that expression has its minimum/maximum (you'll have to guess from the names which does which)

vinz April 12, 2012 03:08
Dear Bernhard,

Thank you for your reply. Indeed it looks like minPosition and maxPosition could be a good starting point.

I went to the wiki in order to see how to download the development version. I got it and apparently managed to install in on my openfoam version 2.1.

I looked at the examples and did not find one including the min or max Position function. Could you tell me if something like this should work:
Code:
vorticityMax {
        type swakExpression;
        verbose true;
        valueType pointField;
        setName vorticityMax;
        expression "maxPosition(magVorticity)";
    }
What does the function return? is it a vector with x y z coordinate?
Or do you have an example?

vinz April 12, 2012 09:08
I got it working. At least I can track one vortex center with this code:
Code:
functions
(
    maxVortPosY
    {
        type swakExpression;
        valueType internalField;
        expression "maxPosition(mag(vorticity)).y";
        accumulations (
          max
        );
        verbose true;       
    }
    maxVortPosX
    {
        type swakExpression;
        valueType internalField;
        expression "maxPosition(mag(vorticity)).x";
        accumulations (
          max
        );
        verbose true;       
    }

);
Now I am trying to do an XY plot instead of X versus time using pyFoamPlotWatcher.

gschaider April 12, 2012 15:52
Quote:
Originally Posted by vinz (Post 354385)
I got it working.
Great. Sorry for a lack of an example. I tested that feature with another case (that I can't publish). But if you have a simple example I'll be happy to include it into swak

Quote:
Originally Posted by vinz (Post 354385)
At least I can track one vortex center with this code:
More will be difficult with the min/max-approach. Unless you "mask out" the region of the first vortex

Quote:
Originally Posted by vinz (Post 354385)
Now I am trying to do an XY plot instead of X versus time using pyFoamPlotWatcher.
Sorry. PyFoamPlotWatcher can currently only do plots over time. Thanks for the reminder on that missing feature ;) (but of course you can easily plot the data written by it easily with any program you like)

vinz April 13, 2012 01:56
Quote:
Originally Posted by gschaider (Post 354444)
More will be difficult with the min/max-approach. Unless you "mask out" the region of the first vortex
Indeed, this is my problem. I have two counter-rotating vortices propagating in my field. They move in time so they are not always at the same place :(
I am not sure how to handle that. I would need to say something like "track the maximum vorticity only where vorticity.z is positive/negative" but I'm not sure if there is a way to do it.

Quote:
Originally Posted by gschaider (Post 354444)
Sorry. PyFoamPlotWatcher can currently only do plots over time. Thanks for the reminder on that missing feature ;) (but of course you can easily plot the data written by it easily with any program you like)
Nothing to be sorry of, you already do a lot! When you say that I can plot the data written by PyFoamPlotWatcher, are you talking about the file in the different directories swakExpression_*? or is there another way to use the data analysed by PyFoamPlotWatcher in order to create the graph I want?

vinz April 16, 2012 02:38
Dear Bernhard,

do you think it is possible to create a subset of my mesh based on the vorticity value or on a coordinate for instance and then look for a maximum only in that subset?

gschaider April 16, 2012 04:50
Quote:
Originally Posted by vinz (Post 354817)
Dear Bernhard,

do you think it is possible to create a subset of my mesh based on the vorticity value or on a coordinate for instance and then look for a maximum only in that subset?
That's what I meant with "mask out". Something like (this is just a sketch):

vortPos1=maxPosition(vorticity)
maskField=(mag(pos()-vortPos1)<rCut) ? 0 : 1
vortPos2=maxPosition( maskField ? vorticity : -1000)
maskField=(mag(pos()-vortPos2)<rCut) ? 0 : maskField
vortPos3= ....

Of course the criteria for masking could be something more advanced (the problem here is that rCut being to small might give you "vortex centers" that are in the other vortex and being too large might "hide" another vortex center)

Alternative would be "dynamic" cellSets based on an expression, but that feature is not yet implemented (partly because I'm not sure whether the parts of the OF-code where this could be interesting can "cope" with cellSets that change their size)

vinz April 16, 2012 10:32
Ok I see, not straight forward but brilliant! :)

I manage to get something which looks interesting with the following lines:
Code:
    mask
    {
        type swakExpression;
        valueType internalField;
        variables (
            "maskField=(vorticity.z>0 || pos().y<0.3) ? 0 : 1;"
        );
        storedVariables (
            {
                name maskField;
                initialValue "0";
            }
        );
        expression "maskField";
        accumulations (
          max
        );
        verbose true;       
    }


    maxVortPos
    {
        type swakExpression;
        valueType internalField;
        variables (
            "vortPos2= maxPosition((maskField>0) ? magVorticity : -1000);"
        );
        storedVariables (
            {
                name vortPos2;
                initialValue "vector(-0.4,1,0)";
            }
        );
        expression "vortPos2";
        accumulations (
          max
        );
        verbose true;       
    }

    maxVortPosY
    {
        type swakExpression;
        valueType internalField;
        expression "vortPos2.y";
        accumulations (
          max
        );
        verbose true;       
    }

    maxVortPosX
    {
        type swakExpression;
        valueType internalField;
        expression "vortPos2.x";
        accumulations (
          max
        );
        verbose true;       
    }
Now I would like to compute the integral of U on a surface which would be a sphere centered on my vortex.
I tried the following to create such a sphere :
Code:
    distField
    {
        type expressionField;
        outputControl timeStep;
        outputInterval 1;
        fieldName distanceFromVortex;
        expression "mag(pos()-vortPos2)";
        autowrite true;       
    }

    surfaceIsoDist
    {
        type swakExpression;
        valueType surface;
        surfaceName isoVortexDist;
        surface {
            type isoSurface;
            isoField distanceFromVortex;
            isoValue 0.15;
            interpolate true;
        }
        verbose true;
        expression "mag(U)";
        accumulations (
            min
            max
            average
        );
    }
but it crashes complaining that vortPos2 does not exist or is of wrong type although I use it earlier to get vortPos2.x for instance.

Then I am wondering what type is expected at this place? Any idea?

gschaider April 16, 2012 14:04
Quote:
Originally Posted by vinz (Post 354904)
Ok I see, not straight forward but brilliant! :)
I think it is the other way round: straight forward (brute force) but not brilliant/elegant

Quote:
Originally Posted by vinz (Post 354904)

Now I would like to compute the integral of U on a surface which would be a sphere centered on my vortex.
I tried the following to create such a sphere :
Code:
    distField
    {
        type expressionField;
        outputControl timeStep;
        outputInterval 1;
        fieldName distanceFromVortex;
        expression "mag(pos()-vortPos2)";
        autowrite true;       
    }

    surfaceIsoDist
    {
        type swakExpression;
        valueType surface;
        surfaceName isoVortexDist;
        surface {
            type isoSurface;
            isoField distanceFromVortex;
            isoValue 0.15;
            interpolate true;
        }
        verbose true;
        expression "mag(U)";
        accumulations (
            min
            max
            average
        );
    }
but it crashes complaining that vortPos2 does not exist or is of wrong type although I use it earlier to get vortPos2.x for instance.

Then I am wondering what type is expected at this place? Any idea?
I think you misunderstood the "storedVariable": this means that the variable keeps its value between timesteps. What you want is either

a) putting vortPos2 in an expression field
b) using a global variable (see slide 35ff in http://openfoamwiki.net/images/d/da/...Leoben2011.pdf).

Something like (untested)

Code:
   
calculateVortex {
        type calculateGlobalVariables;
        outputControl timeStep;
        outputInterval 1;
        valueType internalField;
        toGlobalNamespace vortexLoc;
        toGlobalVariables (
          vortexPos2
        );
        variables (
            "vortPos2= maxPosition((maskField>0) ? magVorticity : -1000);"
        );
}
Then make the scope known with

Code:
globalScopes ( vortexLoc );
where you want to use that location. See the circulatingSplash-example

vinz April 19, 2012 02:47
Thanks again for your help Bernhard!

I managed to get it working, and to track my two vortices.
Here is an example of what I was trying to do:
http://i.imgur.com/oRo3n.gif

gschaider April 19, 2012 07:23
Quote:
Originally Posted by vinz (Post 355528)
Thanks again for your help Bernhard!

I managed to get it working, and to track my two vortices.
Here is an example of what I was trying to do:
http://i.imgur.com/oRo3n.gif
Nice. Could you include this as "swak4Foam was not designed to do this - searching for vortices" on the Wiki-page?

BTW: looking at the picture it occurred to me that by searching for the minimum and the maximum of "vorticity*pos().x" the "masking" could have been avoided

vinz April 19, 2012 10:52
Honestly I have never user the wikki so I have no idea how to include this :)

However, it is not urgent, it's an ongoing work so I cannot share everything right now. The contract ends on july and I should have other interesting (hopefully) thing to show afterwards and place it on the wikki.

Regarding the other possibility to track the vortices, I think it is less general. And although it might work in this specific case, I have other cases with wind comming from the left or where one vortex might turn around the other and I think your last proposal would not work anymore.

Thanks again for your help.

gschaider April 23, 2012 19:18
Quote:
Originally Posted by vinz (Post 355650)
Honestly I have never user the wikki so I have no idea how to include this :)
That is Wiki with one K. Unless you're talking about the company

It's quite easy. Really

Quote:
Originally Posted by vinz (Post 355650)
However, it is not urgent, it's an ongoing work so I cannot share everything right now. The contract ends on july and I should have other interesting (hopefully) thing to show afterwards and place it on the wikki.

Regarding the other possibility to track the vortices, I think it is less general. And although it might work in this specific case, I have other cases with wind comming from the left or where one vortex might turn around the other and I think your last proposal would not work anymore.

Thanks again for your help.
Well. The problem with your solution is that if you're unlucky the vortices can oszillate (exchange positions): if the vorticities for both of them are almost the same then for some time vortex A has the "first maximum", then B, then A again .... and in your graph you'll see the vortices exchanging place.

But I guess there is no general solution that will always work


All times are GMT -4. The time now is 22:03.