[dEntropy]

(*updated code below*)

(*I've included, near the top, the code that didn't work, and I have the final, working, code at the bottom.*)

[dEntropy]

...and my searching has gotten me to this point (it was all much worse when I started).

i appreciate the help. thank you.

[dEntropy]

well, I've played around with all sorts of things. The current conditions are as below:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    
    inlet1
        {
            type            groovyBC;
        valueExpression    #codeStream {
                        codeInclude
                           #{
                                        
                          #};

                        codeOptions
                           #{
                            -I$(LIB_SRC)/groovyBC/lnInclude
                            -I$(LIB_SRC)/swakFunctionObjects/lnInclude
                            -I$(LIB_SRC)/swakTopoSources/lnInclude
                            -I$(LIB_SRC)/simpleSearchableSurfaces/lnInclude
                            -I$(LIB_SRC)/swak4FoamParsers/lnInclude
                            -I$(LIB_SRC)/swakSourceFields/lnInclude                                
                        #};

                        code
                        #{
                            
                            //SI units, please
                            const scalar Reynolds=60;
                            const scalar a=0.0012;
                            const scalar b=0.0012;
                            const scalar rho=997.1;
                            const scalar eta=0.00091;
                            const scalar Pi=3.14159265359;
                            const scalar beta=3*Reynolds*pow(eta,2)*a*b/(4*b*pow(a,3)*rho)*pow(1-(192*0.003825*a)/(pow(Pi,5)*b),-1);                            
                            int s=0;
                            
                            
                            
                            for (int i=1;i<=99999;i+=2) 
                            {
                                s+=(pow((-1),(.5*i-.5))*(1-Re(cosh((i*Pi*fpos().z)/(2*a)))/Re(cosh((i*Pi*b)/(2*a))))*(cos(i*Pi*fpos().y/(2*a))/pow(i,3)));
                            }
                            test=-(16*pow(a,2)*beta/(eta*pow(Pi,3)))*s*normal();
                            os << test;

                        #};
                        };
            value            uniform (0 0 0);
        }
    
    
    outlet1
        {
            type            zeroGradient;
        }

    

        fixedWalls
        {
            type            fixedValue;
            value           uniform (0 0 0);
        }

}



// ************************************************************************* //

Error:

Reading field p

Reading field U

Using #codeStream at line 28 in file "/home/mboehm/OpenFOAM/mboehm-2.4.0/swak4Foam_2.x/Examples/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1"
Using #codeStream with "/home/mboehm/OpenFOAM/mboehm-2.4.0/swak4Foam_2.x/Examples/straightDuct_VelocityProfileInletBC/dynamicCode/platforms/linux64GccDPOpt/lib/libcodeStream_1578c42d908c22e7136dbfd27e503fcb4c57 f809.so"
Creating new library in "dynamicCode/_1578c42d908c22e7136dbfd27e503fcb4c57f809/platforms/linux64GccDPOpt/lib/libcodeStream_1578c42d908c22e7136dbfd27e503fcb4c57 f809.so"
Invoking "wmake -s libso /home/mboehm/OpenFOAM/mboehm-2.4.0/swak4Foam_2.x/Examples/straightDuct_VelocityProfileInletBC/dynamicCode/_1578c42d908c22e7136dbfd27e503fcb4c57f809"
linux64GccDPOpt/options:60: *** commands commence before first target. Stop.
wmake error: file 'Make/linux64GccDPOpt/objectFiles' could not be created in /home/mboehm/OpenFOAM/mboehm-2.4.0/swak4Foam_2.x/Examples/straightDuct_VelocityProfileInletBC/dynamicCode/_1578c42d908c22e7136dbfd27e503fcb4c57f809


--> FOAM FATAL IO ERROR:
Failed wmake "dynamicCode/_1578c42d908c22e7136dbfd27e503fcb4c57f809/platforms/linux64GccDPOpt/lib/libcodeStream_1578c42d908c22e7136dbfd27e503fcb4c57 f809.so"


file: /home/mboehm/OpenFOAM/mboehm-2.4.0/swak4Foam_2.x/Examples/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1 from line 27 to line 27.

From function functionEntries::codeStream::execute(..)
in file db/dictionary/functionEntries/codeStream/codeStream.C at line 177.

FOAM exiting

[dEntropy]

my obvious problems are i don't know C++. seriously, i like openFoam, but i just cannot do the simple calculations i need.

in mathematica, easy. matlab, easy. c++, impossible.

[dEntropy]

Let me explain what I'm trying to do:

I have a quasi-analytical solution for the velocity profile in a duct (that's the set of equations). Unfortunately, this solution comprises an infinite sum of trigonometric functions that depend on the y & z positions (because the velocity profile in 3D has only one non-zero term, i.e., Ux, which is dependent on two spatial positions, i.e., y & z). I want openFoam to do the calculations for me.

Easy as would involve me using Mathematica to solve the equation, saving in a text file then using that text file as the inlet conditions. However, this take two programs and increases the likelyhood that I use the wrong constants in the equation.

So, I thought that pos().y & pos().z was my way to put each individual y & z position into the equation, solve and have a velocity at every location on the inlet face.

I thought normal() made my value a vector, with Uy & Uz components = 0.

[dEntropy]

woops, I guess pos() access internal regions (because it expects a vector). Since I'm applying a boundary condition, I've changed to fpos() so as to access the surface. Sound right?

[gschaider]

Originally Posted by dEntropy

woops, I guess pos() access internal regions (because it expects a vector). Since I'm applying a boundary condition, I've changed to fpos() so as to access the surface. Sound right?

The problem is that you're trying something that groovyBC is there to help you avoid: C++ and codeStreams. The valueExpression is a simple string in which you write an expression. Something like "(1-pos().x)*(pos().x+1)" to get a parabola.

Have a look at http://openfoamwiki.net/index.php/Co...er_information (the presentation from the 9th workshop might be a good start).

But codestreams don't work with groovyBC. If you want to use them use the coded-BC that comes with OpenFOAM

[dEntropy]

Thank you, Bernhard I read through the slides. I can tell you've put a substantial amount of time into this.

I'm afraid I wasn't able to solve my problem, however. It would seem that people take for granted that fpos() or pos() or normal() work in any given code. Unfortunately, I receive the following errors (associated code included below):



/home/mboehm/OpenFOAM/mboehm-2.4.0/run/MySimulations/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1:45:49: error: ‘fpos’ was not declared in this scope

/home/mboehm/OpenFOAM/mboehm-2.4.0/run/MySimulations/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1:47:4: error: ‘os’ was not declared in this scope

/home/mboehm/OpenFOAM/mboehm-2.4.0/run/MySimulations/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1:47:55: error: ‘normal’ was not declared in this scope

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    
    inlet1
        {
            type            codedFixedValue;
        value            uniform (0 0 0);
        redirectType    parabolicSquareProfile;

        code
        #{
        
            //SI units, please
            const scalar Reynolds=60;
            const scalar a=0.0012;
            const scalar b=0.0012;
            const scalar rho=997.1;
            const scalar eta=0.00091;
            const scalar Pi=3.14159265359;
            const scalar beta=3*Reynolds*pow(eta,2)*a*b/(4*b*pow(a,3)*rho)*pow(1-(192*0.003825*a)/(pow(Pi,5)*b),-1);                            
            scalar s=0;
                        
            for (int i=1;i<=99999;i+=2) 
            {
                s+=(pow((-1),(.5*i-.5))*(1-(cosh((i*Pi*fpos().z)/(2*a)))/(cosh((i*Pi*b)/(2*a))))*(cos(i*Pi*fpos().y/(2*a))/pow(i,3)));
            }
            os << -(16*pow(a,2)*beta/(eta*pow(Pi,3)))*s*normal();
        
        #};    
        
        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude                                
        #};
        
    };
}
    
    
    outlet1
        {
            type            zeroGradient;
        }

    

        fixedWalls
        {
            type            fixedValue;
            value           uniform (0 0 0);
        }

}



// ************************************************************************* //

[gschaider]

Originally Posted by dEntropy

Thank you, Bernhard I read through the slides. I can tell you've put a substantial amount of time into this.

I'm afraid I wasn't able to solve my problem, however. It would seem that people take for granted that fpos() or pos() or normal() work in any given code. Unfortunately, I receive the following errors (associated code included below):



/home/mboehm/OpenFOAM/mboehm-2.4.0/run/MySimulations/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1:45:49: error: ‘fpos’ was not declared in this scope

/home/mboehm/OpenFOAM/mboehm-2.4.0/run/MySimulations/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1:47:4: error: ‘os’ was not declared in this scope

/home/mboehm/OpenFOAM/mboehm-2.4.0/run/MySimulations/straightDuct_VelocityProfileInletBC/0/U.boundaryField.inlet1:47:55: error: ‘normal’ was not declared in this scope

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    
    inlet1
        {
            type            codedFixedValue;
        value            uniform (0 0 0);
        redirectType    parabolicSquareProfile;

        code
        #{
        
            //SI units, please
            const scalar Reynolds=60;
            const scalar a=0.0012;
            const scalar b=0.0012;
            const scalar rho=997.1;
            const scalar eta=0.00091;
            const scalar Pi=3.14159265359;
            const scalar beta=3*Reynolds*pow(eta,2)*a*b/(4*b*pow(a,3)*rho)*pow(1-(192*0.003825*a)/(pow(Pi,5)*b),-1);                            
            scalar s=0;
                        
            for (int i=1;i<=99999;i+=2) 
            {
                s+=(pow((-1),(.5*i-.5))*(1-(cosh((i*Pi*fpos().z)/(2*a)))/(cosh((i*Pi*b)/(2*a))))*(cos(i*Pi*fpos().y/(2*a))/pow(i,3)));
            }
            os << -(16*pow(a,2)*beta/(eta*pow(Pi,3)))*s*normal();
        
        #};    
        
        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude                                
        #};
        
    };
}
    
    
    outlet1
        {
            type            zeroGradient;
        }

    

        fixedWalls
        {
            type            fixedValue;
            value           uniform (0 0 0);
        }

}



// ************************************************************************* //

The problem is that this idioms are not standard OpenFOAM-C++ but shorthands that only work in swak to make life easier. At first they have to be called with the current patch prefixed (something like this->patch().foo() but I'd have to look up how to access the current patch in coded myself). Then pos() is a shorthand for the faceCentres()-method of the patch (fpos() only makes sense on the internalField). And accessing components with .x won't work in C++ in that form either

[dEntropy]

First, thanks to those commenting. Also, thanks to Christian at OIST--he wrote the pivotal pieces of the codestream.

So, I went with standard openFoam (i.e., no swak4Foam) and used the #codeStream functionality (code shown below). I also had to be certain about the formatting near the codeOptions library includes because my compiler settings wasn't properly tabbing in the output inside the DynamicCode folder (which is made when running pisoFoam). The code below solves the analytical solution of the Poisson equation, applies that to the inlets ...and bob's your uncle. I hope others find this useful.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    
inlet1
{
  type   codedFixedValue;
  value  uniform (0 0 0);
  redirectType    parabolicSquareProfile1;

  code
  #{
    const double Reynolds = 30;
    const double width = 0.0012;
    const double height = 0.0012;
    const double a=0.5*width;
    const double b=0.5*height;
    const double maxZ=.0012;
    const double minZ=0;
    const double maxX=0.0132;
    const double minX=0.012;
    const double Dh=4*width*height/(2*width + 2*height);
    const double Area=(width*height);
    const double rho=997.1;
    const double eta=0.00091;
    const double Pi=3.14159265359;
    const double beta=(3*Area*Reynolds*std::pow(eta,2))/(4*std::pow(a,3)*b*rho*Dh*0.421731);

    // We start with a fixedValueFvPatchField.
    const fvPatch& patch = this->patch();
    const vectorField &centres = patch.Cf();
    const double zCentral=maxZ-0.5*(maxZ-minZ);
    const double xCentral=maxX-0.5*(maxX-minX);
    
    Field<vector> newValues(this->patchInternalField());

    forAll(newValues, c)
    {
      vector centre = centres[c];                
      const double z = centre.z();
      const double x = centre.x();
      double s = 0;
      for (int i=1;i<=99;i+=2) 
      {
        s+= (std::pow((-1),(.5*i-.5))*(1-(cosh((i*Pi*std::abs(z-zCentral))/(2*a)))/(cosh((i*Pi*b)/(2*a))))*(cos(i*Pi*std::abs(x-xCentral)/(2*a))/std::pow(i,3)));
      }
      
      double value = (16*std::pow(a,2)*beta*s) / (eta*std::pow(Pi,3));
      newValues[c] = vector(0, value, 0);
    }
            
      this->operator==(newValues);
   #};


   codeOptions
   #{
       -I$(LIB_SRC)/finiteVolume/lnInclude \
       -I$(LIB_SRC)/meshTools/lnInclude
   #};

   codeInclude
   #{
      #include "fvCFD.H"
      #include <cmath>
      #include <iostream>
     
   #};

}

inlet2
{
  type   codedFixedValue;
  value  uniform (0 0 0);
  redirectType    parabolicSquareProfile2;

  code
  #{
    const double Reynolds = 30;
    const double width = 0.0012;
    const double height = 0.0012;
    const double a=0.5*width;
    const double b=0.5*height;
    const double maxZ=.0012;
    const double minZ=0;
    const double maxX=0.0132;
    const double minX=0.012;
    const double Dh=4*width*height/(2*width + 2*height);
    const double Area=(width*height);
    const double rho=997.1;
    const double eta=0.00091;
    const double Pi=3.14159265359;
    const double beta=(3*Area*Reynolds*std::pow(eta,2))/(4*std::pow(a,3)*b*rho*Dh*0.421731);

    // We start with a fixedValueFvPatchField.
    const fvPatch& patch = this->patch();
    const vectorField &centres = patch.Cf();
    const double zCentral=maxZ-0.5*(maxZ-minZ);
    const double xCentral=maxX-0.5*(maxX-minX);
    
    Field<vector> newValues(this->patchInternalField());

    forAll(newValues, c)
    {
      vector centre = centres[c];                
      const double z = centre.z();
      const double x = centre.x();
      double s = 0;
      for (int i=1;i<=99;i+=2) 
      {
        s+= (std::pow((-1),(.5*i-.5))*(1-(cosh((i*Pi*std::abs(z-zCentral))/(2*a)))/(cosh((i*Pi*b)/(2*a))))*(cos(i*Pi*std::abs(x-xCentral)/(2*a))/std::pow(i,3)));
      }
      
      double value = -1*(16*std::pow(a,2)*beta*s) / (eta*std::pow(Pi,3));
      newValues[c] = vector(0, value, 0);
    }
            
      this->operator==(newValues);
   #};


   codeOptions
   #{
       -I$(LIB_SRC)/finiteVolume/lnInclude \
       -I$(LIB_SRC)/meshTools/lnInclude
   #};

   codeInclude
   #{
      #include "fvCFD.H"
      #include <cmath>
      #include <iostream>
     
   #};

}


        outlet1
        {
          type  zeroGradient;
        }

        outlet2
        {
          type  zeroGradient;
        }
        
        fixedWalls
        {
          type  fixedValue;
          value uniform (0 0 0);
        }

}



// ************************************************************************* //

[dEntropy]

One note: when you name the codeStream (e.g., parabolicflowprofile1 or parabolicflowprofile2) you need to be certain to name each instance something different. if you do not, the solver will still solve but all boundaries sharing that codestream will also share values for, e.g., velocity. So, the code solves...but solves wrong.