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] topoSet - swakTopoSources - expressionToFace - problem regarding "variables" (https://www.cfd-online.com/Forums/openfoam-community-contributions/145235-toposet-swaktoposources-expressiontoface-problem-regarding-variables.html)

saimat December 1, 2014 12:13
topoSet - swakTopoSources - expressionToFace - problem regarding "variables"
 
Hello foamers,

I'd like to calculate the momentum flow through a special plane in my domain. This plane is formed by the faces of a special set of cells. At the moment I try to define a faceSet containing these faces. The following topoSetDict works well:

Code:
libs (
    "libsampling.so"
    "libsimpleSwakFunctionObjects.so"
    "libswakFunctionObjects.so"
    "libgroovyBC.so"
    "libswakTopoSources.so"
    );

actions
(
    {
        name    x_equals_halfXMax;
        type    faceSet;
        action  new;
        source  expressionToFace;
        sourceInfo
        {
            expression "(fpos().x>interpolate(0.0025-0.04/10*0.0025)) && (fpos().x<interpolate(0.0025+0.04/10*0.0025))"; //plane at x=0.5xmax (real cell faces)
        }
    }

 );
Now I tried to use variables inside the expression to get a more general formulation concerning the length of the given domain:

Code:
libs (
    "libsampling.so"
    "libsimpleSwakFunctionObjects.so"
    "libswakFunctionObjects.so"
    "libgroovyBC.so"
    "libswakTopoSources.so"
    );

theVariables
(
        "xMaxHalf=max(fpos().x)/2;"  //half the length of the domain in x direction
        "deltaX=0.04*0.0025;"        //length of cells in x direction
        "eps=deltaX/10;"
);

actions
(
    {
        name    x_equals_halfXMax;
        type    faceSet;
        action  new;
        source  expressionToFace;
        variables $theVariables;
        sourceInfo
        {
            expression "(fpos().x>xMaxHalf-eps) && (xMaxHalf+eps)"; //plane at x=0.5xmax (real cell faces)
        }
    }

 );
This unfortunately gives an error when running topoSet:

Code:
Create polyMesh for time = 0

Reading topoSetDict

Time = 0
    mesh not changed.
Created faceSet x_equals_halfXMax
    Applying source expressionToFace
    Adding all elements of for which "(fpos().x>xMaxHalf-eps) && (xMaxHalf+eps)" evaluates to true ...
swak4Foam: Allocating new repository for sampledMeshes
swak4Foam: Allocating new repository for sampledGlobalVariables
--> FOAM Warning :
    From function ConcretePluginFunction<DriverType>::exists
    in file lnInclude/ConcretePluginFunction.C at line 121
    Constructor table of plugin functions for FieldValueExpressionDriver is not initialized


--> FOAM FATAL ERROR:
 Parser Error for driver FieldValueExpressionDriver at "1.11-18" :"field xMaxHalf not existing or of wrong type"
"(fpos().x>xMaxHalf-eps) && (xMaxHalf+eps)"
            ^^^^^^^^
------------|     

Context of the error:


- Driver constructed from scratch
  Evaluating expression "(fpos().x>xMaxHalf-eps) && (xMaxHalf+eps)"


    From function parsingValue
    in file lnInclude/CommonValueExpressionDriverI.H at line 1181.

FOAM exiting
The error appears also when xMaxHalf is defined as "xMaxHalf=0.0025", which should, in my opinion, produce the same expression as in the hard-coded topoSetDict.

In the last days I tried to get this running - but with no success. Does anyone know about how to get rid of this error? What am I doing wrong?

gschaider December 1, 2014 16:29
Quote:
Originally Posted by saimat (Post 521903)
Hello foamers,

I'd like to calculate the momentum flow through a special plane in my domain. This plane is formed by the faces of a special set of cells. At the moment I try to define a faceSet containing these faces. The following topoSetDict works well:



Now I tried to use variables inside the expression to get a more general formulation concerning the length of the given domain:



This unfortunately gives an error when running topoSet:



The error appears also when xMaxHalf is defined as "xMaxHalf=0.0025", which should, in my opinion, produce the same expression as in the hard-coded topoSetDict.

In the last days I tried to get this running - but with no success. Does anyone know about how to get rid of this error? What am I doing wrong?
As a rule of thumb the "variables" entry has to be on the same level as the "expression". Try moving it into the sourceInfo-subDict

BTW: for sources and output use the CODE-tag. Not QUOTE.

saimat December 2, 2014 04:44
Thanks a lot for this quick and valuable reply. Is this written down somewhere so that I could have avoided asking this stupid question?
This runs for now:

Code:
theVariables
(
        "xMaxHalf=0.5*max(pos().x);"  //half the length of the domain in x direction
        "deltaX=0.04*0.0025;"        //length of cells in x direction
        "eps=deltaX/10;"
);

actions
(
    {
        name    x_equals_halfXMax;
        type    faceSet;
        action  new;
        source  expressionToFace;
        sourceInfo
        {
          variables $theVariables;
          expression "(fpos().x>interpolate(xMaxHalf-eps)) && (fpos().x<interpolate(xMaxHalf+eps))"; //plane at x=0.5xmax (real cell faces)
        }
    }

 );
Please note that I had to change "xMaxHalf=0.5*max(fpos().x);" into "xMaxHalf=0.5*max(pos().x);". The previous one gave the following error:

Code:
--> FOAM FATAL ERROR:
 Parser Error for driver FieldValueExpressionDriver at "1.9-12" :"syntax error, unexpected TOKEN_fposition"
"0.5*max(fpos().x)"
          ^^^^
----------| 

Context of the error:


- Driver constructed from scratch
  Evaluating expression "0.5*max(fpos().x)"


    From function parsingValue
    in file lnInclude/CommonValueExpressionDriverI.H at line 1181.
I don't understand this error because fpos() is used at the same level one line downwards without any error. Unfortunately I do not get the full length of my domain when using pos() instead. Do you know how to fix this?

gschaider December 2, 2014 07:11
Quote:
Originally Posted by saimat (Post 522024)
Thanks a lot for this quick and valuable reply. Is this written down somewhere so that I could have avoided asking this stupid question?
This runs for now:

Code:
theVariables
(
        "xMaxHalf=0.5*max(pos().x);"  //half the length of the domain in x direction
        "deltaX=0.04*0.0025;"        //length of cells in x direction
        "eps=deltaX/10;"
);

actions
(
    {
        name    x_equals_halfXMax;
        type    faceSet;
        action  new;
        source  expressionToFace;
        sourceInfo
        {
          variables $theVariables;
          expression "(fpos().x>interpolate(xMaxHalf-eps)) && (fpos().x<interpolate(xMaxHalf+eps))"; //plane at x=0.5xmax (real cell faces)
        }
    }

 );
Please note that I had to change "xMaxHalf=0.5*max(fpos().x);" into "xMaxHalf=0.5*max(pos().x);". The previous one gave the following error:

Code:
--> FOAM FATAL ERROR:
 Parser Error for driver FieldValueExpressionDriver at "1.9-12" :"syntax error, unexpected TOKEN_fposition"
"0.5*max(fpos().x)"
          ^^^^
----------| 

Context of the error:


- Driver constructed from scratch
  Evaluating expression "0.5*max(fpos().x)"


    From function parsingValue
    in file lnInclude/CommonValueExpressionDriverI.H at line 1181.
I don't understand this error because fpos() is used at the same level one line downwards without any error. Unfortunately I do not get the full length of my domain when using pos() instead. Do you know how to fix this?
The problem is that "0.5" is a value defined at cells and "fpos().x" is a value at face. The solution here would be "interpolate(0.5)" (see http://sourceforge.net/p/openfoam-ex...amReference.md) so that both of them are on the same "page".

Anyway: your way of selecting the faces (those within a small layer by selecting their face-value) is a bit unstable for several reasons.

Expression to face works in two modes depending on the result of expression: if the result is a face-field (the "fpos()"-variant) then "true" means "take this face". If the result is a cell-field it is a bit different. Then a face is taken only if one cell is True and the other is False. So selecting a plane at halfHeight is simply "pos().x<halfHeight". This second way makes it easier to select a "closed" plane in any grid

saimat December 4, 2014 04:13
Quote:
Expression to face works in two modes depending on the result of expression: if the result is a face-field (the "fpos()"-variant) then "true" means "take this face". If the result is a cell-field it is a bit different. Then a face is taken only if one cell is True and the other is False. So selecting a plane at halfHeight is simply "pos().x<halfHeight". This second way makes it easier to select a "closed" plane in any grid
Tank you for this great explanation. I changed the code with respect to your advice. "topoSetDict" now looks like:

Code:
theVariables
(
        "xHalf=min(pos().x)+0.5*(max(pos().x)-min(pos().x));"  //center of the domain in x direction
);

actions
(
    {
        name    x_equals_halfXMax;
        type    faceSet;
        action  new;
        source  expressionToFace;
        sourceInfo
        {
          variables $theVariables;
          expression "pos().x<xHalf"; //plane at x=0.5xmax (real cell faces)
        }
    }

 );
The calculation of the momentum flux through the plane (faceSet) is not the main topic of this thread, but maybe some other users find it interesting. I did this via swakExpression in the controlDict.:

Code:
    meanMomentumFlux
    {
        type swakExpression;
        valueType faceSet;
        setName x_equals_halfXMax;            //constructed in topoSet (see corresponding dictionary)
        autoInterpolate true;                //interpolate from cell center to cell face
        accumulations (
            sum                              //sum over faceset (vertical plane in the middle of the channel)
        );
        outputControlMode timeStep;
        outputInterval 20;
        after 0.1;                            //start after begining to calculate the averages (timeStart - see above)
        expression "997.535*UMean.x*UMean.x*area()/sum(area())";  //calculate mean momentum flux per area
        verbose true;
    }
Please note that in my case all the faces are found inside a plane defined by x=const. I tried a solution with flip(), phi, normal() and U, but it didn't work. Flip() command for example resulted in an error telling me about some missing slaveCells and I didn't know where to get them from.

gschaider December 4, 2014 15:54
Quote:
Originally Posted by saimat (Post 522419)
Tank you for this great explanation. I changed the code with respect to your advice. "topoSetDict" now looks like:

Code:
theVariables
(
        "xHalf=min(pos().x)+0.5*(max(pos().x)-min(pos().x));"  //center of the domain in x direction
);

actions
(
    {
        name    x_equals_halfXMax;
        type    faceSet;
        action  new;
        source  expressionToFace;
        sourceInfo
        {
          variables $theVariables;
          expression "pos().x<xHalf"; //plane at x=0.5xmax (real cell faces)
        }
    }

 );
The calculation of the momentum flux through the plane (faceSet) is not the main topic of this thread, but maybe some other users find it interesting. I did this via swakExpression in the controlDict.:

Code:
    meanMomentumFlux
    {
        type swakExpression;
        valueType faceSet;
        setName x_equals_halfXMax;            //constructed in topoSet (see corresponding dictionary)
        autoInterpolate true;                //interpolate from cell center to cell face
        accumulations (
            sum                              //sum over faceset (vertical plane in the middle of the channel)
        );
        outputControlMode timeStep;
        outputInterval 20;
        after 0.1;                            //start after begining to calculate the averages (timeStart - see above)
        expression "997.535*UMean.x*UMean.x*area()/sum(area())";  //calculate mean momentum flux per area
        verbose true;
    }
Please note that in my case all the faces are found inside a plane defined by x=const. I tried a solution with flip(), phi, normal() and U, but it didn't work. Flip() command for example resulted in an error telling me about some missing slaveCells and I didn't know where to get them from.
faceSet does not store the flip-vector (faceZone does that) so swak4foam determines it (this is a convention borrowed from the usual OF-utilities) by looking for a cellSet with the same name but SlaveCells added and assume that these cells are on one side of the faceSet. Your fix: create a cellSet with the appropriate name. As an expression you can use exactly the same you used to create the faceSet

saimat December 5, 2014 04:17
Thank you very much for helping. I'll keep this in mind (actually I took notes).


All times are GMT -4. The time now is 00:59.