[swak4Foam] topoSet - swakTopoSources - expressionToFace - problem regarding "variables"

saimat · December 1, 2014, 13:13

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, 17:29

Quote:

Originally Posted by saimat

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, 05: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, 08:11

Quote:

Originally Posted by saimat

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, 05: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, 16:54

Quote:

Originally Posted by saimat

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, 05:17

Thank you very much for helping. I'll keep this in mind (actually I took notes).

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Scheme Programming: problem writing "set!" for multiple variables	nsaeidi	Fluent UDF and Scheme Programming	1	June 16, 2016 15:58
Sod Shock Tube problem	cfd	Main CFD Forum	4	July 7, 2009 14:24
Fortran 77 problem	CFD_nerd	Main CFD Forum	5	February 23, 2009 08:44
Adiabatic and Rotating wall (Convection problem)	ParodDav	CFX	5	April 29, 2007 20:13
extremely simple problem... can you solve it properly?	Mikhail	Main CFD Forum	40	September 9, 1999 10:11

December 5, 2014, 05:17		#7
saimat New Member Peter im Hof Join Date: Nov 2013 Posts: 9 Rep Power: 13	Thank you very much for helping. I'll keep this in mind (actually I took notes).