[S.Colucci]

Dear Philip,
thanks for your reply! I tried by simplifying the code, but I get the same error

lnInclude/Henry.C:132:53: error: passing ‘const Foam::fvPatchField<double>’ as ‘this’ argument of ‘void Foam::fvPatchField<Type>:perator=(const Foam::UList<T>&) [with Type = double]’ discards qualifiers [-fpermissive]
IndexTest[0][0].boundaryField()[patchi] = newPatchValues;

This is the code:

// -----------------------------------------------
PtrList<PtrList<volScalarField> > IndexTest(1);

IndexTest.set(0, new PtrList<volScalarField>(2));

IndexTest[0].set
(
0,
new volScalarField
(
IOobject
(
"IndexTest" ,
this->pair_.phase1().time().timeName(),
this->pair_.phase1().mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
this->pair_.phase1().mesh(),
dimensionedScalar("IndexTest", dimless, 1)
)
);

Info<<"IndexTest[0][0] "<<IndexTest[0][0]<<endl;

// correct BC
forAll(IndexTest[0][0].boundaryField(), patchi)
{
scalarField newPatchValues(IndexTest[0][0].boundaryField()[patchi].size(),scalar(0));

forAll(IndexTest[0][0].boundaryField()[patchi],i)
{
newPatchValues[i] = floor(IndexTest[0][0].boundaryField()[patchi][i]);
}

IndexTest[0][0].boundaryField()[patchi] = newPatchValues;
}
Info<<"IndexTest[0][0] "<<IndexTest[0][0]<<endl;
// -----------------------------------------------

I'm using OpenFOAM 4, the solver is reactingEulerFoam, the code above is in applications/solvers/multiphase/reactingEulerFoam/InterfacialCompositionModels/interfaceCompositionModels/Henry/Henry.C in the member function

template<class Thermo, class OtherThermo>
Foam::tmp<Foam::volScalarField>
Foam::interfaceCompositionModels::Henry<Thermo, OtherThermo>::Yf
(
const word& speciesName,
const volScalarField& Tf
) const
{

The only const I see is in the function, but I do not know if it related or not...

Simone

[S.Colucci]

I tried also to copy the code above outside the function, directly in a *.H in the folder reactingTwoPhaseEulerFoam in order to substitute

IOobject
(
"IndexTest" ,
this->pair_.phase1().time().timeName(),
this->pair_.phase1().mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),

with

IOobject
(
"IndexTest",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),

but I get the same error

error: passing ‘const Foam::fvPatchField<double>’ as ‘this’ argument of ‘void Foam::fvPatchField<Type>:perator=(const Foam::UList<T>&) [with Type = double]’ discards qualifiers [-fpermissive]
IndexTest[0][0].boundaryField()[patchi] = newPatchValues;

Any idea?

Thanks

Simone

[albet]

Dear Philip,
I have followed your comments on how to update the BC contents with coding in solver.
I want to write a Dirichlet BC with nonuniform values for omega at the wall. The omega should be calculated from friction velocity for each cell of the wall (I calculated it from U field). My solver is simpleFoam. I use the method you mentioned below for calculating omega in the simpleFoam solver.
I have written my code but the problem is I don't know how can I insert the calculated values in omega field?
(BTW I am not experienced in coding in c++)

any hints are really appreciated.
Many thanks in advance.
Amir

Originally Posted by bigphil

Hi thomas99,

As marupio said you could write your own boundary condition,

alternatively, you could just include a header file which updates the boundary contact inside your solver solution loop,
for example I sometimes use the "fixedGradient" boundary condition in my solid stress solver to apply a traction, this means I must update this fixedGradient every solution iteration.
The header file might look something like the following (where U is the volVectorField I solve for):

label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest_name");
if
(
  U.boundaryField()[patchID].type()
   == fixedGradientFvPatchVectorField::typeName
)
{
          fixedGradientFvPatchVectorField& Upatch =
        refCast<fixedGradientFvPatchVectorField>(U.boundaryField()[patchID]);
      
      Upatch.gradient() = (.......calculate something here.....);
}

Hope it helps,
Philip

[bigphil]

Hi Amir,

Can you give your code that you have so far? And where are you inserting this code?

Philip

[albet]

Originally Posted by bigphil

Hi Amir,

Can you give your code that you have so far? And where are you inserting this code?

Philip

Thank you for your reply,
I don't think my code structure is easily understandable.
First I calculated the velocity and the distance from the wall for the centre of the cells on the wall.
Then, I used the dimensionless van Driest near-wall velocity profile and used an iterative method (Newton-Raphson method) to calculate the friction velocity in each cell face of the wall, and calculated omega at the wall by the method that has been mentioned in Wilcox (2006) for both smooth and rough walls.
the problem is I don't know how can I insert the calculated omega into omega field. Because the omega has not been defined in the main solver and it is defined in turbulence models.
I simply tried to define a volScalarField for omega in the createFields.H and saving the calculated omega in it, but it seems that the solver does not use it.
I started my code exactly after the starting loop of the simpleFoam

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

         // I started my code from here.

Maybe I should write a new BC or write my code in a different place. I read about codedFixedValue BC as well, but I don't know that it can be used in this case!
Thank you very much for your help in advance,

Regards,
Amir

[bigphil]

Hi Amir,

The cleanest method would be to write a new boundary condition for omega where you can directly update the boundary values for omega.

However, it is possible to do it by adding code to your solver: one way to do it would be to "look up" the omega field from the turbulence model using the objectRegistry (this keeps track of all the significant objects) and then update it using your proposed method (however, this method requires const_cast which is evil: const cast is evil).
For example, you could do the following:

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // Lookup omega field from the objectRegistry and use const_cast to
        // remove the "read only" protection
        volScalarField& omega =
            const_cast<volScalarField&>
            (
                mesh.lookupObject<volScalarField>("omega")
            );

        // Find the patch of interest
        const label patchID = mesh.boundaryMesh().findPatchID("myInterestingPatch");

        // Check the patch was found
        if (patchID == -1)
        {
            FatalError
                << "patch was not found!" << abort(FatalError);
        }

        // Cast the patch field to let the solver know what class the boundary condition is
        // I am guessing it is a wall function patch for omega?
        fixedInternalValueFvPatchScalarField& omegaPatch =
            refCast<fixedInternalValueFvPatchScalarField>(omega.boundaryField()[patchID]);

        // Set values on/near the patch
        forAll(omegaPatch, faceI)
        {
            // Set values here
            // omegaPatch.refValue()[faceI] = ...;
        }


        // code for solving U and p equations etc.
        // ...
    }

Also, you would need to add the following header file at the top of the solver e.g. just under #include "fvCFD.H":

#include "fixedInternalValueFvPatchFields.H"

Note: I have not tried to compile the code above so there may be some errors but the overall method should work.

But if you are trying to create an omega wall function boundary condition, I would suggest making a copy of omegaWallFunctionFvPatchScalarField and renaming to myOmegaWallFunctionFvPatchScalarField and then you can edit the code in the updateCoeffs() function.

Philip

[albet]

Originally Posted by bigphil

Hi Amir,

The cleanest method would be to write a new boundary condition for omega where you can directly update the boundary values for omega.

However, it is possible to do it by adding code to your solver: one way to do it would be to "look up" the omega field from the turbulence model using the objectRegistry (this keeps track of all the significant objects) and then update it using your proposed method (however, this method requires const_cast which is evil: const cast is evil).
For example, you could do the following:

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // Lookup omega field from the objectRegistry and use const_cast to
        // remove the "read only" protection
        volScalarField& omega =
            const_cast<volScalarField&>
            (
                mesh.lookupObject<volScalarField>("omega")
            );

        // Find the patch of interest
        const label patchID = mesh.boundaryMesh().findPatchID("myInterestingPatch");

        // Check the patch was found
        if (patchID == -1)
        {
            FatalError
                << "patch was not found!" << abort(FatalError);
        }

        // Cast the patch field to let the solver know what class the boundary condition is
        // I am guessing it is a wall function patch for omega?
        fixedInternalValueFvPatchScalarField& omegaPatch =
            refCast<fixedInternalValueFvPatchScalarField>(omega.boundaryField()[patchID]);

        // Set values on/near the patch
        forAll(omegaPatch, faceI)
        {
            // Set values here
            // omegaPatch.refValue()[faceI] = ...;
        }


        // code for solving U and p equations etc.
        // ...
    }

Also, you would need to add the following header file at the top of the solver e.g. just under #include "fvCFD.H":

#include "fixedInternalValueFvPatchFields.H"

Note: I have not tried to compile the code above so there may be some errors but the overall method should work.

But if you are trying to create an omega wall function boundary condition, I would suggest making a copy of omegaWallFunctionFvPatchScalarField and renaming to myOmegaWallFunctionFvPatchScalarField and then you can edit the code in the updateCoeffs() function.

Philip

Thank you Philip,
I have learned a lot from this post and this thread.
Actually, I want to use omega on the wall (y=0), and not as a wall function, that is why I have not tried changing the omegaWallFunction.
In this case, I should use codedFixedValue instead of fixedInternalValueFvPatchScalarField in your code?

Regards,
Amir

[bigphil]

Hi Amir,

OK, so if you intend to set the value of omega at the patch faces then I suggest you use the fixedValue patch type (for your patch in 0/omega); then the code becomes more simple:

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // Lookup omega field from the objectRegistry and use const_cast to
        // remove the "read only" protection
        volScalarField& omega =
            const_cast<volScalarField&>
            (
                mesh.lookupObject<volScalarField>("omega")
            );

        // Find the patch of interest
        const label patchID = mesh.boundaryMesh().findPatchID("myInterestingPatch");

        // Check the patch was found
        if (patchID == -1)
        {
            FatalError
                << "patch was not found!" << abort(FatalError);
        }

        // Take a reference to the patch field (values at the patch faces)
        scalarField& omegaPatch = omega.boundaryField()[patchID];

        // You might find some of these patch functions useful when setting the omega
        // values
        // Patch face centres
        const vectorField& patchCf = mesh.boundary()[patchID].Cf();
        // Patch face area vectors
        const vectorField& patchSf = mesh.boundary()[patchID].Sf();
        // Patch delta vectors (from face centre to closest cell centre)
        // Note: this returns a value not a reference (no '&')
        const vectorField patchDelta = mesh.boundary()[patchID].delta();

        // Set face values the patch
        forAll(omegaPatch, faceI)
        {
            // Set values here
            // It could be based on patchCf[faceI], patchSf[faceI], patchDelta[faceI], etc.
            // omegaPatch[faceI] = ...;
        }


        // code for solving U and p equations etc.
        // ...
    }

There is no need to include any extra header files because we can directly manipulate the value on a patch without casting to the specific boundary condition type (though if boundary conditions other than fixedValue are used then they might overwrite this value).

Philip

[Kombinator]

Dear Philip,

Hope you have been monitoring this thread. I try to compile the following code (based on yours) for fixedGradient but I have the following error. I wonder if you could tell what the problem is.
I tryed "#include "fixedGradientFvPatchFields.H.H" but it didn't help.

Code:

**body of solver is here**

    label patchID = mesh.boundaryMesh().findPatchID("frontAndBack");
    if(patchID == -1)   
    {
      FatalError << "patch not found!" << exit(FatalError);      
    }
    
    if
      (
	PotE.boundaryField()[patchID].type()
	== fixedGradientFvPatchVectorField::typeName
      )
      
{
          fixedGradientFvPatchVectorField& PotEpatch =
	    refCast<fixedGradientFvPatchVectorField>(PotE.boundaryField()[patchID]);
	  
	  PotEpatch.gradient() = vector(0,0,0);
}

Error:

In file included from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/token.H:46:0,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/UListIO.C:28,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/UList.C:233,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/UList.H:484,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/List.H:43,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/labelList.H:48,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/UPstream.H:42,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/Pstream.H:42,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/parRun.H:35,
                 from /data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/finiteVolume/lnInclude/fvCFD.H:4,
                 from mhdEpot5Foam.C:37:
/data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/typeInfo.H: In instantiation of ‘To& Foam::refCast(From&) [with To = Foam::fixedGradientFvPatchField<Foam::Vector<double> >; From = const Foam::fvPatchField<double>]’:
mhdEpot5Foam.C:157:76:   required from here
/data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/typeInfo.H:110:35: error: cannot dynamic_cast ‘r’ (of type ‘const class Foam::fvPatchField<double>’) to type ‘class Foam::fixedGradientFvPatchField<Foam::Vector<double> >&’ (conversion casts away constness)
         return dynamic_cast<To&>(r);
                                   ^
/data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/typeInfo.H:119:35: error: cannot dynamic_cast ‘r’ (of type ‘const class Foam::fvPatchField<double>’) to type ‘class Foam::fixedGradientFvPatchField<Foam::Vector<double> >&’ (conversion casts away constness)
         return dynamic_cast<To&>(r);
                                   ^
/data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/typeInfo.H: In function ‘To& Foam::refCast(From&) [with To = Foam::fixedGradientFvPatchField<Foam::Vector<double> >; From = const Foam::fvPatchField<double>]’:
/data/OpenFoams/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/typeInfo.H:121:1: warning: control reaches end of non-void function [-Wreturn-type]
 }
 ^
make: *** [Make/linux64GccDPInt64Opt/mhdEpot5Foam.o] Error 1

I have found that the problem is in the following two lines (if I remove them everything is ok), but I don't know what to do. Thank you.

          fixedGradientFvPatchVectorField& PotEpatch =
	    refCast<fixedGradientFvPatchVectorField>(PotE.boundaryField()[patchID]);

UPDATE: I changed boundaryFielt to boundaryFieldRef, and now it complies. But did I do right thing (I don't understand why it helped). Thanks

Best regards,
Tar

[bigphil]

Hi Tar,

Yep, it seems using boundaryFieldRef is the correct thing to do, e.g. as discussed here link.

Philip

[Kombinator]

Dear Philip,

Thank you very much for your answer. Your link helped me a lot to understand the reason.

Best regards,
Tar

[Kombinator]

Dear Philip,

Sorry to bother you again. I faced one more problem. My code is following:

label patchID = mesh.boundaryMesh().findPatchID("frontAndBack");
  
if(patchID == -1)
  
  {
    FatalError << "patch not found!" << exit(FatalError);      
  }
    
if
  (
    PotE.boundaryFieldRef()[patchID].type()
    == fixedGradientFvPatchScalarField::typeName
  )
      
  {
    fixedGradientFvPatchScalarField& PotEpatch = refCast<fixedGradientFvPatchScalarField>(PotE.boundaryFieldRef()[patchID]);  
    forAll(PotEpatch, faceI)
      {
        PotEpatch.gradient()[faceI] = PotE;
      }
    
  }

But when I try to compile I get the following error:

dynamicBC.H: In function ‘int main(int, char**)’:
dynamicBC.H:23:35: error: cannot convert ‘Foam::volScalarField {aka Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>}’ to ‘double’ in assignment
       PotEpatch.gradient()[faceI] = PotE;
                                   ^
make: *** [Make/linux64GccDPInt64Opt/mhdEpot5Foam.o] Error 1

As I understood from different topics, this happens due to different dimensions of the variable (for instance if I substitute PotE (volScalarField) with 0 or 1 or any number - it complies, but only with numbers). But unfortunately I didn't find how to overcome this problem. I will be very thankful if you can advice something about this (may be another reason exists). Thanks in advance.

A also noticed that if I remove [faceI] from this line:

 PotEpatch.gradient()[faceI] = PotE;

then it complies. Maybe the problem is not connected only with dimensions.

Best regards,
Tar

[Kombinator]

One more update:

When I tried this

      PotEpatch.gradient()[faceI] = PotE.internalField()[faceI];

the code has complied. To be honest I don'y understand why and I suspect that it is not right, because I use internalField. Any help will be appreciated . Thanks in advance.

Best regards,
Tar

[bigphil]

Originally Posted by Kombinator

One more update:

When I tried this

      PotEpatch.gradient()[faceI] = PotE.internalField()[faceI];

the code has complied. To be honest I don'y understand why and I suspect that it is not right, because I use internalField. Any help will be appreciated . Thanks in advance.

Best regards,
Tar

Hi Tar,

What are you trying to set as the boundary gradient?
It is not clear to me what you are actually trying to do.

Philip

[Kombinator]

Originally Posted by bigphil

Hi Tar,

What are you trying to set as the boundary gradient?
It is not clear to me what you are actually trying to do.

Philip

Hi Philip,

Thank you for your answer. The general goal is implementation the following dynamic boundary condition:
= a
where is electric potential (Scalar Field, PotE in openFoam code, I acquire this value on the wall every iteration and then I update BC), a is constant, is tangential Laplacian at the wall.
But right now I try to implement something more simple in order to understand the algorithm. So of course, this equation

      PotEpatch.gradient()[faceI] = PotE.internalField()[faceI];

=
doesn't make sense. But I use this just to see how I can realize this connection in openFoam.

Best regards,
Tar

[bigphil]

OK, I understand.

Something like the following should work:

if
(
    PotE.boundaryFieldRef()[patchID].type()
 == fixedGradientFvPatchScalarField::typeName
)
{
    fixedGradientFvPatchScalarField& PotEpatch = 
        refCast<fixedGradientFvPatchScalarField>(PotE.boundaryFieldRef()[patchID]);  

    // PotEpatch is also a list of face values at the patch
    PotEpatch.gradient() = PotEpatch;
}

Obviously, this is an explicit implementation and outer iterations are required i.e. the boundary gradient here depends on the previous iteration of the boundary value as opposed to the current value. Also, for this style of approach, under-relaxation may be required to achieve convergence.

By the way, the following does not make sense:

PotE.internalField()[faceI]

because PotE.internalField() is a scalar field for the cells (not the faces); faceI is a boundary face index, which is not directly related to cell indices.

[Kombinator]

Originally Posted by bigphil

Something like the following should work:

Dear Philip,

Thank you very much for your advice. It works now and I understood my error and algorithm. I also would like to ask two more questions (hope last questions) if you don't mind:
1. In your suggestion I get uniform boundary condition, but I would like to get nonuniform. I wrote the following code to cope with this problem, it compiles and seems logical but I kindly ask you to take a look (maybe I still have some silly mistakes):

if
  (
    PotE.boundaryFieldRef()[patchID].type()
    == fixedGradientFvPatchScalarField::typeName
  )
      
  {
    fixedGradientFvPatchScalarField& PotEpatch = refCast<fixedGradientFvPatchScalarField>(PotE.boundaryFieldRef()[patchID]);  
    forAll(PotEpatch, faceI)
      
    {
      PotEpatch.gradient()[faceI] = PotEpatch[faceI];
    }
    
  }

2. The second question is about tangential and normal parts of the gradient on the wall. As I understood in openFoam we usually have an access and calculate only the normal part of the gradient. Is it possible to have an access and calculate the tangential part of the gradient on the wall (the gradient along the boundary)?
For instance I would like to realize something like this:
=

You wrote some information about the tangential part in this post Updating Boundary Conditions Each Iteration but to be honest I am not fully sure how to do this.

Thanks in advance again.

Best regards,
Tar

[bigphil]

Your proposed code for applying a spatially varying condition look fine to me.

As regards the tangential gradient at the boundary, OpenFOAM uses the tangential gradient at the closest cell-centre (so-called patch-internal-field) by default i.e. extrapolates tangential gradient from the internal field using a zero gradient approach. it is possible to directly calculate the tangential gradient at the boundary using, for example, the finite area method (you can google this for more information); however, I suggest you try the former method first.

An example approach could be:

// Calculate the gradient of PotE everywhere
const volVectorField gradPotE = fvc:grad(PotE);

fixedGradientFvPatchScalarField& PotEpatch =
    refCast<fixedGradientFvPatchScalarField>
    (
        PotE.boundaryFieldRef()[patchID]
    );

// Get patch unit normals
const vectorField nf = PotEpatch.patch().nf();

// Get tangential component of gradPotE at the boundary
const vectorField tangGradPotE =
    ((I - sqr(nf)) & gradPotE.boundaryFieldRef()[patchID].patchInternalField());

// Set the normal gradient equal to the magnitude of the tangential gradient for each face
forAll(PotEpatch, faceI)
{
    PotEpatch.gradient()[faceI] = mag(tangGradPotE[faceI]);
}

[Kombinator]

Dear Philip,

Thank you for your answer, it helps me a lot. But when I try to compile your code it shows the following error (openFoam 4.1):

dynamicBC.H: In function ‘int main(int, char**)’:
dynamicBC.H:27:48: error: passing ‘const volVectorField {aka const Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>}’ as ‘this’ argument of ‘Foam::GeometricField<Type, PatchField, GeoMesh>::Boundary& Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryFieldRef() [with Type = Foam::Vector<double>; PatchField = Foam::fvPatchField; GeoMesh = Foam::volMesh]’ discards qualifiers [-fpermissive]
     ((I - sqr(nf)) & gradPotE.boundaryFieldRef()[patchID].patchInternalField());
                                                ^
make: *** [Make/linux64GccDPInt64Opt/mhdEpot5Foam.o] Error 1

I have seen this problem before in this topic, so I removed "const" from all lines and now it complies (I hope it was the right solution).
Now I will try to do the same for the tangential laplacian. When I get something I will share my results in this topic. Thank you very much for your help again.

Best regards,
Tar

[bigphil]

Hi Tar,

Removing the "const" in this case will fix it; also, you could use boundaryField instead of boundaryFieldRef for the boundary of gradPotE. Both will act the same.

Philip