[vinz]

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:

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:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
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]

Originally Posted by vinz

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:

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:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
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]

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:

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]

I got it working. At least I can track one vortex center with this 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]

Originally Posted by vinz

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

Originally Posted by vinz

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

Originally Posted by vinz

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]

Originally Posted by gschaider

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.

Originally Posted by gschaider

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]

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]

Originally Posted by vinz

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]

Ok I see, not straight forward but brilliant!

I manage to get something which looks interesting with the following lines:

    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 :

    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]

Originally Posted by vinz

Ok I see, not straight forward but brilliant!

I think it is the other way round: straight forward (brute force) but not brilliant/elegant

Originally Posted by vinz

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 :

    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)

    
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

globalScopes ( vortexLoc );

where you want to use that location. See the circulatingSplash-example

[vinz]

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:

[gschaider]

Originally Posted by vinz

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:

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]

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]

Originally Posted by vinz

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

Originally Posted by vinz

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