patchInternalField() returning a random Field

Benben · July 12, 2018, 08:18

Hello,

I am using a “codedFixedValue” boundary condition with reactingFoam, on OpenFoam 5.x.
I am trying to set the composition of “H2” at the boundary as a function of some other properties in the cell next to the boundary.
My “Case/0/H2” file contains the following boundary condition :

Code:

fixedWalls
    {
        type            codedFixedValue;
        value           uniform 1.;
        redirectType    EvaporatingSurface;
        code
        #{
            //Getting the list of faces that are part of the boundary patch
            const fvPatch& boundaryPatch = patch();//boundary patch
            const vectorField& Cf = boundaryPatch.Cf();//List of face centers
            
            //read properties at the boundary patch
            scalarField& w_H2 = *this;
            const scalarField& w_N2 = this->patch().lookupPatchField<volScalarField, scalar>("N2");
            const scalarField& P = this->patch().lookupPatchField<volScalarField, scalar>("p");
            const scalarField& T = this->patch().lookupPatchField<volScalarField, scalar>("T");
            
            //Reading properties in the cells next to the boundary patch
            const scalarField& w_H2_inf = this->patchInternalField();
            const scalarField& w_N2_inf = this->patch().lookupPatchField<volScalarField, scalar>("N2").patchInternalField();
            const scalarField& P_inf = this->patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();
            const scalarField& T_inf = this->patch().lookupPatchField<volScalarField, scalar>("T").patchInternalField();
            
            //Looping on the id of every faces of the patch
            forAll(Cf, faceID)
            {
              //printing the properties that were read, for every face of the patch
              std::cout << std::endl;
              std::cout << std::endl;
              std::cout << "w_H2 = " << w_H2[faceID] << std::endl;
              std::cout << "w_N2 = " << w_N2[faceID] << std::endl;
              std::cout << "T = " << T[faceID] << std::endl;
              std::cout << "P = " << P[faceID] << std::endl;
              std::cout << std::endl;
              std::cout << "w_H2_inf = " << w_H2_inf[faceID] << std::endl;
              std::cout << "w_N2_inf = " << w_N2_inf[faceID] << std::endl;
              std::cout << "T_inf = " << T_inf[faceID] << std::endl;
              std::cout << "P_inf = " << P_inf[faceID] << std::endl;
            }
        #};
    }

This allows me to print to screen the values i want to use for my calculations. However, the values returned are not the one i expected, and change from one execution of "reactingFoam" to another. For example, one of the results was :

Code:

w_H2 = 1
w_N2 = 0.5
T = 600
P = 101325
　
w_H2_inf = 101325
w_N2_inf = 1
T_inf = 600
P_inf = 101325

It seems that the mass fraction of H2 and N2 get randomly mixed up properties. Even weirder, these values change between two executions of reactingFoam.

here is another try :

Code:

w_H2 = 1
w_N2 = 0.5
T = 600
P = 101325
　
w_H2_inf = 600
w_N2_inf = 101325
T_inf = 600
P_inf = 101325

Does someone has any idea of what is going on ?
Thank you for your help.

Benben · July 12, 2018, 08:28

Here is a test case to reproduce :
reactingFoam_TestBug.zip

simrego · July 12, 2018, 09:08

Hi!

I think this boundary condition is not a correct one. You set nothing, just read data. Honestly i don't know what's happening in the solver in this case, maybe it just use the initial value as a fixedValue. You should set something on the boundary, not just read data.
Another thing. Are you sure you check the same face? Because you go through all the faces, and are you sure that the supported boundary values are on the same face?

And you don't have to use std. In OF you can just use:
Info<< "Your text" << endl;

Benben · July 12, 2018, 09:55

Hello simrego, thank you for your help.

Indeed, i dont set anything in this example, it is only a test case to see if i am able to get the right values for the properties i am looking for.

I figured out that the problem was that i indexed properties read from inside the mesh with the same label (faceID in my code) as the properties from the boundary.
In this case, how can i get the adequate label object for the cell next to my boundary ?

What i need to find is how i can acess the "label" object of the cell from which "patchInternalField" is taking the value.

Benben · July 12, 2018, 10:08

I found that fvPatchField.patch().faceCells() returns the list of the labels of the cells associated to the patch faces. I will try to see if :

Code:

//Getting the list of faces that are part of the boundary patch
const fvPatch& boundaryPatch = patch();//boundary patch
const vectorField& Cf = boundaryPatch.Cf();//List of face centers
 
//Reading the "H2" fields of the cell next to each face of the patch
const scalarField& w_H2_inf = this->patchInternalField();
 
//Looping on the id of every faces of the patch
forAll(Cf, faceID)
{
         cellIDs = this->patch().faceCells()
         cellID = cellIDs[faceID]
         Info << "w_H2_inf = " << w_H2_inf[cellID] << endl;
}

Does the trick.

Benben · July 12, 2018, 10:28

I modified my boundary condition as follow :

Code:

fixedWalls
    {
        type            codedFixedValue;
        value           uniform 1.;
        redirectType    EvaporatingSurface;
        code
        #{
            //Getting the list of faces that are part of the boundary patch
            const fvPatch& boundaryPatch = patch();//boundary patch
            const vectorField& Cf = boundaryPatch.Cf();//List of face centers
            
            //read properties at the boundary patch
            scalarField& w_H2 = *this;
            const scalarField& w_N2 = this->patch().lookupPatchField<volScalarField, scalar>("N2");
            const scalarField& P = this->patch().lookupPatchField<volScalarField, scalar>("p");
            const scalarField& T = this->patch().lookupPatchField<volScalarField, scalar>("T");
            
            //Reading properties in the cells next to the boundary patch
            //const scalarField& w_H2_inf = this->patchInternalField();
            const scalarField& w_H2_inf = this->patch().lookupPatchField<volScalarField, scalar>("rho").patchInternalField();
            const scalarField& w_N2_inf = this->patch().lookupPatchField<volScalarField, scalar>("N2").patchInternalField();
            const scalarField& P_inf = this->patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();
            const scalarField& T_inf = this->patch().lookupPatchField<volScalarField, scalar>("T").patchInternalField();
            
            //Looping on the id of every faces of the patch
            //label faceID;
            forAll(Cf, faceID)
            {
              //Reading the ID of the cell next to the current face
              label faceCellID = this->patch().faceCells()[faceID];
              
              //printing the properties that were read, for every face of the patch
              std::cout << std::endl;
              std::cout << std::endl;
              std::cout << "w_H2 = " << w_H2[faceID] << std::endl;
              std::cout << "w_N2 = " << w_N2[faceID] << std::endl;
              std::cout << "T = " << T[faceID] << std::endl;
              std::cout << "P = " << P[faceID] << std::endl;
              std::cout << std::endl;
              std::cout << "w_H2_inf = " << w_H2_inf[faceCellID] << std::endl;
              std::cout << "w_N2_inf = " << w_N2_inf[faceCellID] << std::endl;
              std::cout << "T_inf = " << T_inf[faceCellID] << std::endl;
              std::cout << "P_inf = " << P_inf[faceCellID] << std::endl;
            }
        #};
    }

But the value returned from the cells next to faces are still wrong :/

simrego · July 12, 2018, 10:49

Hi!

Try this. I'm not so familiar with using internal cell value for BC, but i think it is working.
And as i see if you don't set anything on the boundary, the initial value will be used. (At least i guess that based on my trials.)

Code:

         code
        #{
            const scalarField& p =
                patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();            

            forAll(p, i)
            {
                Info<< p[i] << endl;
            }
        #};

Benben · July 12, 2018, 11:22

Thank you for your help simrego.

However the boundary condition i want to implement depends on properties at the boundary, and at the cell next to the boundary at the same time. (As if, for example, i wanted to set the pressure at the boundary so that pressure gradient be constant)

I need something that does that :

Code:

#{
    const scalarField& P_face = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();
    const scalarField& P_cell = patchInternalField();
    forAll(P_face, faceIndex)
    {
        CellIndex = patch().faceCells()[faceIndex];
        Info << "Pressure at the face of the boundary :" << P_face[faceIndex] << endl;
        Info << "Pressure in the cell next to the boundary :" << P_cell[CellIndex] << endl;
        Info << "Delta P at boundary :" << P_face[faceIndex] - P_cell[CellIndex] << endl;
    }
#};

I already tried something similar to this, but it didnt worked. The "faceIndex" i am getting doesnt appear to be the Index of the cell next to the face of ID faceID.

Benben · July 12, 2018, 11:29

I just tried it, and it seems to be working actually oO

The problem was that i was looping with for all on the wrong variable, aka :

Code:

forAll(Cf,i)
{
}

instead of

Code:

forAll(scalar,i)
{
}

I will try this with my bigger code, thank you very much for your help.

simrego · July 12, 2018, 11:34

Yes, I understand your problem. You can do like this:

Code:

 code
#{
        //Value on face:
        const scalarField& P_face = patch().lookupPatchField<volScalarField, scalar>("p");
        //Value on cell near patch
        const scalarField& P_cell = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField()
        forAll(P_face, faceIndex)
        {
             //Here do what you want with P_face and P_cell
         }


  #};

But if you have to define a given gradient, and value, you can use codedMixedFvPatch field where you can define a value, a gradient, and a fraction:
https://cpp.openfoam.org/v6/classFoa...d.html#details

You can loop on face centres or face values, i think you will get the same. If you check what is forAll:
https://cpp.openfoam.org/v6/UList_8H...0e9b03d336c6a9
It is just use the length of the 1st entry.

Benben · July 12, 2018, 11:44

Finally it didnt work

Here is my test code

Code:

fixedWalls
    {
      type            codedFixedValue;
      value           uniform 1.;
      redirectType    EvaporatingSurface;
      code
      #{
        const scalarField& P_face = patch().lookupPatchField<volScalarField, scalar>("p");
        const scalarField& P_cell = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();
        forAll(P_face, faceIndex)
        {
            //Index of the cell next to the face
            const label& CellIndex = patch().faceCells()[faceIndex];
            
            Info << "Pressure at the face of the boundary :" << P_face[faceIndex] << endl;
            Info << "Pressure in the cell next to the boundary :" << P_cell[CellIndex] << endl;
            
            const scalar Delta = P_face[faceIndex] - P_cell[CellIndex];
            Info << "Delta P at boundary :" << Delta << endl;
        }
      #};
    }

I was expecting to get a deltaP of 0 (as i initialized the P inside and at the boundary to the same values), but i get random values for the pressure in the cell next ot the boundary.

Your code seems to be working thought, here is my test code with your version :

Code:

fixedWalls
    {
      type            codedFixedValue;
      value           uniform 1.;
      redirectType    EvaporatingSurface;
      code
      #{
        const scalarField& P_face = patch().lookupPatchField<volScalarField, scalar>("p");
        const scalarField& P_cell = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();
        forAll(P_face, faceIndex)
        {           
            Info << "Pressure at the face of the boundary :" << P_face[faceIndex] << endl;
            Info << "Pressure in the cell next to the boundary :" << P_cell[faceIndex] << endl;
            
            const scalar Delta = P_face[faceIndex] - P_cell[faceIndex];
            Info << "Delta P at boundary :" << Delta << endl;
        }
      #};
    }

Will test it with some more variables, like H2 and N2

Benben · July 12, 2018, 12:05

It seems that this doesnt work with multiple variables, for example :

Code:

fixedWalls
    {
      type            codedFixedValue;
      value           uniform 1.;
      redirectType    EvaporatingSurface;
      code
      #{
        const scalarField& H2_face = *this;
        const scalarField& H2_cell = patchInternalField();
        const scalarField& P = patch().lookupPatchField<volScalarField, scalar>("p");
        
        const scalarField& N2_face = patch().lookupPatchField<volScalarField, scalar>("N2");
        const scalarField& N2_cell = patch().lookupPatchField<volScalarField, scalar>("N2").patchInternalField();
        const scalarField& P_inf = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();
        
        forAll(H2_face, faceIndex)
        {   
            Info << endl;
            Info << endl;
            Info << "N2 at the face of the boundary :" << N2_face[faceIndex] << endl;
            Info << "N2 in the cell next to the boundary :" << N2_cell[faceIndex] << endl;
            Info << endl;
            Info << "H2 at the face of the boundary :" << H2_face[faceIndex] << endl;
            Info << "H2 in the cell next to the boundary :" << H2_cell[faceIndex] << endl;
            Info << endl;
            Info << "P at the face of the boundary :" << P[faceIndex] << endl;
            Info << "P in the cell next to the boundary :" << P_inf[faceIndex] << endl;
        }
      #};
    }

sometimes returns abnormal things like :

Code:

N2 at the face of the boundary :0.5
N2 in the cell next to the boundary :101325
　
H2 at the face of the boundary :1
H2 in the cell next to the boundary :0
　
P at the face of the boundary :101325
P in the cell next to the boundary :101325

And it changes from one run of reactingFoam to another
Here is the test case :
reactingFoam_TestBug.zip

Benben · July 13, 2018, 03:59

I found out what was wrong !

I was actually doing some dirty type castings, and the compiler was not telling me anything about it.

I was storing the value returned by "patchInternalField()" into a "scalarField&". However, according to the doc it should return a "tmp<scalarField>". The compiler didnt tell me anuthing and when i tired to access the wrong object values, it was returning me some garbage.

Here is the working code :

Code:

fixedWalls
    {
        type            codedFixedValue;
        value           uniform 1.;
        redirectType    EvaporatingSurface;
        code
        #{
            //read properties at the boundary patch
            scalarField& w_H2 = *this;
            const scalarField& w_N2 = patch().lookupPatchField<volScalarField, scalar>("N2");
            const scalarField& P = patch().lookupPatchField<volScalarField, scalar>("p");
            const scalarField& T = patch().lookupPatchField<volScalarField, scalar>("T");
            
            //Reading properties in the cells next to the boundary patch
            const tmp<scalarField>& w_H2_inf = patchInternalField();
            const tmp<scalarField>& w_N2_inf = patch().lookupPatchField<volScalarField, scalar>("N2").patchInternalField();
            const tmp<scalarField>& P_inf = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField();
            const tmp<scalarField>& T_inf = patch().lookupPatchField<volScalarField, scalar>("T").patchInternalField();
            
            //Looping on the id of every faces of the patch
            forAll(w_H2, faceID)
            {
              //printing the properties that were read, for every face of the patch
              Info << endl;
              Info << endl;
              Info << "w_H2 = " << w_H2[faceID] << endl;
              Info << "w_N2 = " << w_N2[faceID] << endl;
              Info << "T = " << T[faceID] << endl;
              Info << "P = " << P[faceID] << endl;
              Info << endl;
              Info << "w_H2_inf = " << w_H2_inf.ref()[faceID] << endl;
              Info << "w_N2_inf = " << w_N2_inf.ref()[faceID] << endl;
              Info << "T_inf = " << T_inf.ref()[faceID] << endl;
              Info << "P_inf = " << P_inf.ref()[faceID] << endl;
            }
        #};

And the whole test case for people who have to deal with this issues :
reactingFoam_CodedFixedValue_example.zip

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Foam::error::PrintStack	almir	OpenFOAM Running, Solving & CFD	91	December 21, 2022 04:50
problems after decomposing for running	alessio.nz	OpenFOAM	7	March 5, 2021 04:49
creating an internal field?	maybee	OpenFOAM Programming & Development	2	February 4, 2021 17:15
[Other] dynamicTopoFVMesh and pointDisplacement	RandomUser	OpenFOAM Meshing & Mesh Conversion	6	April 26, 2018 07:30
Initial turbulent field	Tony Chen	Main CFD Forum	10	July 19, 1999 00:52

July 12, 2018, 08:28		#2
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 8	Here is a test case to reproduce : reactingFoam_TestBug.zip

July 12, 2018, 09:08		#3
simrego Senior Member anonymous Join Date: Jan 2016 Posts: 416 Rep Power: 14	Hi! I think this boundary condition is not a correct one. You set nothing, just read data. Honestly i don't know what's happening in the solver in this case, maybe it just use the initial value as a fixedValue. You should set something on the boundary, not just read data. Another thing. Are you sure you check the same face? Because you go through all the faces, and are you sure that the supported boundary values are on the same face? And you don't have to use std. In OF you can just use: Info<< "Your text" << endl;

July 12, 2018, 09:55		#4
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 8	Hello simrego, thank you for your help. Indeed, i dont set anything in this example, it is only a test case to see if i am able to get the right values for the properties i am looking for. I figured out that the problem was that i indexed properties read from inside the mesh with the same label (faceID in my code) as the properties from the boundary. In this case, how can i get the adequate label object for the cell next to my boundary ? What i need to find is how i can acess the "label" object of the cell from which "patchInternalField" is taking the value.

July 12, 2018, 10:49		#7
simrego Senior Member anonymous Join Date: Jan 2016 Posts: 416 Rep Power: 14	Hi! Try this. I'm not so familiar with using internal cell value for BC, but i think it is working. And as i see if you don't set anything on the boundary, the initial value will be used. (At least i guess that based on my trials.) Code: code #{ const scalarField& p = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField(); forAll(p, i) { Info<< p[i] << endl; } #};

July 12, 2018, 11:22		#8
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 8	Thank you for your help simrego. However the boundary condition i want to implement depends on properties at the boundary, and at the cell next to the boundary at the same time. (As if, for example, i wanted to set the pressure at the boundary so that pressure gradient be constant) I need something that does that : Code: #{ const scalarField& P_face = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField(); const scalarField& P_cell = patchInternalField(); forAll(P_face, faceIndex) { CellIndex = patch().faceCells()[faceIndex]; Info << "Pressure at the face of the boundary :" << P_face[faceIndex] << endl; Info << "Pressure in the cell next to the boundary :" << P_cell[CellIndex] << endl; Info << "Delta P at boundary :" << P_face[faceIndex] - P_cell[CellIndex] << endl; } #}; I already tried something similar to this, but it didnt worked. The "faceIndex" i am getting doesnt appear to be the Index of the cell next to the face of ID faceID.

July 12, 2018, 11:29		#9
Benben Member benoit favier Join Date: Jun 2017 Posts: 64 Rep Power: 8	I just tried it, and it seems to be working actually oO The problem was that i was looping with for all on the wrong variable, aka : Code: forAll(Cf,i) { } instead of Code: forAll(scalar,i) { } I will try this with my bigger code, thank you very much for your help.

July 12, 2018, 11:34		#10
simrego Senior Member anonymous Join Date: Jan 2016 Posts: 416 Rep Power: 14	Yes, I understand your problem. You can do like this: Code: code #{ //Value on face: const scalarField& P_face = patch().lookupPatchField<volScalarField, scalar>("p"); //Value on cell near patch const scalarField& P_cell = patch().lookupPatchField<volScalarField, scalar>("p").patchInternalField() forAll(P_face, faceIndex) { //Here do what you want with P_face and P_cell } #}; But if you have to define a given gradient, and value, you can use codedMixedFvPatch field where you can define a value, a gradient, and a fraction: https://cpp.openfoam.org/v6/classFoa...d.html#details You can loop on face centres or face values, i think you will get the same. If you check what is forAll: https://cpp.openfoam.org/v6/UList_8H...0e9b03d336c6a9 It is just use the length of the 1st entry.