How to implement non-uniform initial conditions to a field in folder 0

rucky96 · September 26, 2019, 06:29

Hi! I am new in OF so sorry if what I am asking is too basic!

I want to implement non-uniform initial conditions to a field in folder 0. Let's say, for example, U. My system is a cube with cyclic boundaries so I am just worried about the internal Field.

My aim is to define a function dependent on the space U(x,y,z,t=0)=(Ux(x,y,z,0),Uy(x,y,z,0),Uz(x,y,z,0)). I want the input to be the nodes of the mesh (hexaedra, 40cellsx40cellsx40cells) and return the function U(x,y,z,t=0).

My plan was to define a function in a file for Ux, Uy and Uz; call it in the U file in 0 folder and apply it over the nodes of the mesh where internalField is written. The problem is that I do not know how to perform this since I did not find any tutorial in OF doing it. Could any of you help me by telling me where I can find information on how to define a function in OF, take it to the U file at 0 folder and / or apply it on the nodes of my mesh?

PS:I am a Python user and this is the most logical way to do it that I thought, but maybe that is not the right way to do it in OF...

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

FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField #Here the output list;

boundaryField
{

left
{
type cyclic;

}

[..]

back
{
type cyclic;

}

}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Edit I tried to use #codeStream but I get a FATAL ERROR. The code:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.x                                 |
|   \\  /    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   #codeStream //Use codeStream to set the value of the initial conditions
{
    codeInclude
    #{
        #include "fvCFD.H"
    #};
    codeOptions
    #{
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude
    #};
    codeLibs
    #{
        -lmeshTools \
        -lfiniteVolume
    #};
    //Depending of what are you trying to do, you will need to add new files, options and libraries.
    //For most of the cases, this part is always the same.

    code //Insert your code here. At this point, you need to know how to access internal mesh information
    #{
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());
        //Access internal mesh information
       
        vectorField iniU(mesh.nCells(), uniform(0.,0.,0.));//Initialize vector field to zero for all components
       
        forAll(iniU, i) //forAll loop to access cell centers and to assign iniU values. Notice the iniU was previously initialized. The size of the loop is defined by iniU and the iterator is i.
        {
            IOobject Sheader("iniU",runTime.timeName(),mesh,IOobject::MUST_READ);

            const scalar x = mesh.C()[i][0];
            const scalar y = mesh.C()[i][1];
            const scalar z = mesh.C()[i][2];
            //Access cell centers coordinates

            iniU[i][0] = -2*sin(y);
            iniU[i][1] =  2*sin(x);
            iniU[i][2] = 0;
        }
        iniU.write();
    #};
};

boundaryField
{
    left
    {
        type            cyclic;
    }

    right
    {
        type            cyclic;
    }

    down
    {
        type            cyclic;
    }

    up
    {
        type            cyclic;
    }

    front
    {
        type            cyclic;
    }
    back
    {
        type            cyclic;
    }
}
// ************************************************************************* //

The ERROR:

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 7-6ef561967074
Exec   : mhdFoam
Date   : Sep 26 2019
Time   : 11:33:48
Host   : "b483d6a79080"
PID    : 479
I/O    : uncollated
Case   : /home/foamusr/localfolder/tfg/cube2
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field U

Using #codeStream at line 19 in file "/home/foamusr/localfolder/tfg/cube2/0/U"
Using #codeStream with "/home/foamusr/localfolder/tfg/cube2/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libcodeStream_6204b24330a8236a303c9e9f30a6fbc5014e11ef.so"
Creating new library in "dynamicCode/_6204b24330a8236a303c9e9f30a6fbc5014e11ef/platforms/linux64GccDPInt32Opt/lib/libcodeStream_6204b24330a8236a303c9e9f30a6fbc5014e11ef.so"
Invoking "wmake -s libso /home/foamusr/localfolder/tfg/cube2/dynamicCode/_6204b24330a8236a303c9e9f30a6fbc5014e11ef"
wmake libso /home/foamusr/localfolder/tfg/cube2/dynamicCode/_6204b24330a8236a303c9e9f30a6fbc5014e11ef
    ln: ./lnInclude
    wmkdep: codeStreamTemplate.C
    Ctoo: codeStreamTemplate.C
/home/foamusr/localfolder/tfg/cube2/0/U.#codeStream: In function 'void Foam::codeStream_6204b24330a8236a303c9e9f30a6fbc5014e11ef(Foam::Ostream&, const Foam::dictionary&)':
/home/foamusr/localfolder/tfg/cube2/0/U.#codeStream:44:41: error: 'uniform' was not declared in this scope
/home/foamusr/localfolder/tfg/cube2/0/U.#codeStream:44:41: note: suggested alternative: 'union'
/home/foamusr/localfolder/tfg/cube2/0/U.#codeStream:48:37: error: 'runTime' was not declared in this scope
/home/foamusr/localfolder/tfg/cube2/0/U.#codeStream:48:37: note: suggested alternative: 'dimTime'
/home/foamusr/localfolder/tfg/cube2/0/U.#codeStream:52:26: warning: unused variable 'z' [-Wunused-variable]
/home/foamusr/localfolder/tfg/cube2/0/U.#codeStream:59:14: error: 'Foam::vectorField {aka class Foam::Field<Foam::Vector<double> >}' has no member named 'write'
/opt/openfoam7/wmake/rules/General/transform:25: recipe for target 'Make/linux64GccDPInt32Opt/codeStreamTemplate.o' failed
make: *** [Make/linux64GccDPInt32Opt/codeStreamTemplate.o] Error 1


--> FOAM FATAL IO ERROR:
Failed wmake "dynamicCode/_6204b24330a8236a303c9e9f30a6fbc5014e11ef/platforms/linux64GccDPInt32Opt/lib/libcodeStream_6204b24330a8236a303c9e9f30a6fbc5014e11ef.so"


file: /home/foamusr/localfolder/tfg/cube2/0/U from line 17 to line 17.

    From function static void (* Foam::functionEntries::codeStream::getFunction(const Foam::dictionary&, const Foam::dictionary&))(Foam::Ostream&, const Foam::dictionary&)
    in file db/dictionary/functionEntries/codeStream/codeStream.C at line 218.

FOAM exiting

Does anyone see the mistake(s)?

Thanks!

rucky96 · September 26, 2019, 08:15

As an alternative, I tried to implement a #codeStream in the U internalField. I have modified an example of this page http://www.cfd-china.com/assets/uplo...ing_of_ics.pdf. The code for U is:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.x                                 |
|   \\  /    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   #codeStream //Use codeStream to set the value of the initial conditions
{
    {
        codeInclude
        #{
            #include "fvCFD.H"
        #};
        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};
        codeLibs
        #{
            -lmeshTools \
            -lfiniteVolume
        #};
        //Depending of what are you trying to do, you will need to add new files, options and libraries.
        //For most of the cases, this part is always the same.
       
        code //Insert your code here. At this point, you need to know how to access internal mesh information
        #{
            const IOdictionary& d = static_cast<const IOdictionary&>(dict);
            const fvMesh& mesh = refCast<const fvMesh>(d.db());
            //Access internal mesh information
           
            vectorField iniU(mesh.nCells(), uniform(0.,0.,0.));//Initialize vector field to zero for all components
           
            forAll(iniU, i) //forAll loop to access cell centers and to assign iniU values. Notice the iniU was previously initialized. The size of the loop is defined by iniU and the iterator is i.
            {
                IOobject Sheader("iniU",runTime.timeName(),mesh,IOobject::MUST_READ);

                const scalar x = mesh.C()[i][0];
                const scalar y = mesh.C()[i][1];
                const scalar z = mesh.C()[i][2];
                //Access cell centers coordinates

                iniU[i][0] = -2*sin(y);
                iniU[i][1] =  2*sin(x);
                iniU[i][2] = 0;
            }
            iniU.write();
        #};
    };
}

boundaryField
{
    left
    {
        type            cyclic;
    }

    right
    {
        type            cyclic;
    }

    down
    {
        type            cyclic;
    }

    up
    {
        type            cyclic;
    }

    front
    {
        type            cyclic;
    }
    back
    {
        type            cyclic;
    }
}
// ************************************************************************* //

But I get a FATAL ERROR:

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 7-6ef561967074
Exec   : mhdFoam
Date   : Sep 26 2019
Time   : 10:39:22
Host   : "b483d6a79080"
PID    : 473
I/O    : uncollated
Case   : /home/foamusr/localfolder/tfg/cube2
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field U

Using #codeStream at line 19 in file "/home/foamusr/localfolder/tfg/cube2/0/U"
--> FOAM Warning :
    From function entry::New(dictionary&, Istream&)
    in file db/dictionary/entry/entryIO.C at line 144
    Reading /home/foamusr/localfolder/tfg/cube2/0/U
    found on line 21 the punctuation token '{'
    expected either } or EOF
Using #codeStream with "/home/foamusr/localfolder/tfg/cube2/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libcodeStream_fa6bcb9352550c3007ffd513d48777a6a88042f2.so"


--> FOAM FATAL IO ERROR:
expected keyword 'uniform' or 'nonuniform', found codeInclude

file: /home/foamusr/localfolder/tfg/cube2/0/U from line 17 to line 40.

    From function Foam::Field<Type>::Field(const Foam::word&, const Foam::dictionary&, Foam::label) [with Type = Foam::Vector<double>; Foam::label = int]
    in file /home/ubuntu/OpenFOAM/OpenFOAM-7/src/OpenFOAM/lnInclude/Field.C at line 223.

FOAM exiting

Can anyone see the mistake(s) in the code?

Thanks!

piotr.mecht · September 26, 2019, 08:36

It is reasonably easy with swak4Foam installed. You can apply groovyBC type of boundary condition:
https://openfoamwiki.net/index.php/Contrib/groovyBC

You can add this line to your controlDict:

Code:

libs ( "libgroovyBC.so" );

and then your velocity BC can look like this:

Code:

Inlet
    {
        type            groovyBC;
        valueExpression "vector(0.15*pow((1-(sqrt(pow(pos().z,2)+pow(pos().y,2))/0.0115)),0.14),0,0)";
        value           uniform (0.12 0 0); //dummy value for compatibility
    }

missios · February 2, 2022, 12:07

Hi Rucky96,

Perhaps you need to remove the double curly brakets after the internalField line

Code:

internalField   #codeStream //Use codeStream to set the value of the  initial conditions 
{ 

{

After that I would suggest to replace this line

Code:

vectorField iniU(mesh.nCells(), uniform(0.,0.,0.));//Initialize vector field to zero for all components

with this

Code:

vectorField U_0(mesh.nCells(), Foam::Vector<double>(0.,0.,0.));

Best,
K

joshwilliams · February 6, 2022, 12:22

In the past, I have used python to assign non-uniform initial fields (i.e. Taylor-Green flow). I have code for parsing and writing openfoam fields to and from python here (https://github.com/jvwilliams23/OpenFOAMparser).

Code:

import numpy as np
from numpy import pi

from Ofpp import mesh_parser, field_parser


cells = field_parser.parse_internal_field("0/C")
vols = field_parser.parse_internal_field("0/V")
# print("centres", c, c.shape)

nx = np.unique(cells[:,0].round(decimals=6)).size
ny = np.unique(cells[:,1].round(decimals=6)).size
nz = np.unique(cells[:,2].round(decimals=6)).size

xSorted = np.unique(cells[:,0].round(decimals=6))
xSorted.sort()
delta2D = xSorted[1] - xSorted[0]

lx = 1.0
ly = 1.0
lz = 0.4

# create pressure and velocity fields
uf = np.zeros(cells.shape)
p = np.zeros(cells.shape[0])
for i, coords in enumerate(cells):
  x,y,z = coords
  #u = np.cos(2.*pi*coords[0]/lx)*np.cos(2.*pi*coords[1]/ly)
  u = v_scale*np.sin(x/l_scale)*np.cos(y/l_scale)*np.cos(z/l_scale)
  v = 0
  w = -v_scale*np.cos(x/l_scale)*np.cos(y/l_scale)*np.sin(z/l_scale)
  uf[i][0] = u
  uf[i][1] = v
  uf[i][2] = w
  p[i] = v_scale*v_scale/16. * (np.cos(2.*x/l_scale) + np.cos(2.*y/l_scale)) * (np.cos(2.*z/l_scale)+2.)

# setup files for writing

basefiles = listdir("0/")
for bf in basefiles:
    if "air" in bf:
        phaseName = ".air"
        break

copy("0/U.orig","0/U"+phaseName)
copy("0/p.orig","0/p")


field_parser.write_field_data(uf,"0/U"+phaseName)
field_parser.write_field_data(p, "0/p")

It works very well for coarse grids such as yours. For fine grids, there is probably a better way of doing it using some native OpenFOAM tool but I went with the quick and easy option here.

jherb · September 6, 2022, 07:27

Hello Josh,

the code on github does not contain the write method. Do you plan to publish it there? It would be really useful.

Joachim

Quote:

Originally Posted by joshwilliams

In the past, I have used python to assign non-uniform initial fields (i.e. Taylor-Green flow). I have code for parsing and writing openfoam fields to and from python here (https://github.com/jvwilliams23/OpenFOAMparser).

joshwilliams · September 6, 2022, 07:37

Hi Joachim,

I have just checked the repository, and the write_field_data method is there (see line 235 on link to class).
Maybe the issue is that my version is not yet published on pypi, so if using via "pip install openfoamparser", it will not be available. I submitted a pull request to try get my changes to the upstream fork and published on PyPi but no response yet. You could try copying field_parser.py to your working directory and using it directly as a module (i.e. "import field_parser as fp")

Josh

Quote:

Originally Posted by jherb

the code on github does not contain the write method. Do you plan to publish it there? It would be really useful.

jherb · September 6, 2022, 08:53

Thank you for your really fast answer. I installed it by cloning the github repo and then:

Code:

python setup.py sdist
pip install dist/openfoamparser-0.13.tar.gz

Just another question: Can I also write the boundary field or "only" the internal field?

Quote:

Originally Posted by joshwilliams

Hi Joachim,

I have just checked the repository, and the write_field_data method is there (see line 235 on link to class).
Maybe the issue is that my version is not yet published on pypi, so if using via "pip install openfoamparser", it will not be available. I submitted a pull request to try get my changes to the upstream fork and published on PyPi but no response yet. You could try copying field_parser.py to your working directory and using it directly as a module (i.e. "import field_parser as fp")

Josh

joshwilliams · September 6, 2022, 10:27

Ah good question - that is not something I had thought about doing. Currently, it does not support that, but there are methods there to parse the boundary field, so I think it would be simple to implement something to write the boundary field. Feel free to submit an issue to the GitHub repo (or a pull request if you implement your own solution).

Quote:

Originally Posted by jherb

]Just another question: Can I also write the boundary field or "only" the internal field?

September 26, 2019, 08:36		#3
piotr.mecht Member Piotr Ładyński Join Date: Apr 2017 Posts: 55 Rep Power: 9	It is reasonably easy with swak4Foam installed. You can apply groovyBC type of boundary condition: https://openfoamwiki.net/index.php/Contrib/groovyBC You can add this line to your controlDict: Code: libs ( "libgroovyBC.so" ); and then your velocity BC can look like this: Code: Inlet { type groovyBC; valueExpression "vector(0.15*pow((1-(sqrt(pow(pos().z,2)+pow(pos().y,2))/0.0115)),0.14),0,0)"; value uniform (0.12 0 0); //dummy value for compatibility } rucky96 likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
chtMultiRegionSimpleFoam turbulent case	Aditya Patil	OpenFOAM Running, Solving & CFD	6	April 24, 2017 23:13
pressure in incompressible solvers e.g. simpleFoam	chrizzl	OpenFOAM Running, Solving & CFD	13	March 28, 2017 06:49
High Courant Number @ icoFoam	Artex85	OpenFOAM Running, Solving & CFD	11	February 16, 2017 14:40
Upgraded from Karmic Koala 9.10 to Lucid Lynx10.04.3	bookie56	OpenFOAM Installation	8	August 13, 2011 05:03
Differences between serial and parallel runs	carsten	OpenFOAM Bugs	11	September 12, 2008 12:16

February 2, 2022, 12:07		#4
missios New Member Konstantinos Missios Join Date: Mar 2017 Location: Copenhagen, Denmark. Posts: 12 Rep Power: 9	Hi Rucky96, Perhaps you need to remove the double curly brakets after the internalField line Code: internalField #codeStream //Use codeStream to set the value of the initial conditions { { After that I would suggest to replace this line Code: vectorField iniU(mesh.nCells(), uniform(0.,0.,0.));//Initialize vector field to zero for all components with this Code: vectorField U_0(mesh.nCells(), Foam::Vector<double>(0.,0.,0.)); Best, K