[nuttita]

Hi everyone,

I'm trying to implement Coulomb type viscosity model and I use interFoam solver. In calcNu() in my Coulomb.c, I have

const volScalarField& pp=U_.db().objectRegistry::lookupObject<volScalarF ield>("p");

volScalarField nu__=mu0_/(max(sr(),dimensionedScalar("VSMALL",
dimless/dimTime,VSMALL)))*pp/rho_;

When I compiled the library, there is no error. However, when I included the library to run with interFoam I got an error:

--> FOAM FATAL ERROR:

request for volScalarField p from objectRegistry region0 failed
available objects of type volScalarField are

4
(
alpha.water
(1.41421*mag(symm(grad(U))))
p_rgh
alpha.air
)

I understand that interFoam solves for p_rgh (p-rho gh) instead of p, so there is no need for p to present in /0. However, p is written in the later time directories which means p has been calculated.

Any advice how to obtain field p?

Thank a lot!
Nuttita

[floquation]

I don't have time to trace the entire src code for you, but this is my hypothesis (check for yourself if it's true):

In createFields.H of interFoam the following is found:

Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
    IOobject
    (
        "p_rgh",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

(...)

volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    p_rgh + rho*gh
);

That is, both volScalarField's are constructed, but a different constructor is called.
p_rgh is constructed based on a mesh. p is constructed based on a tmp<volScalarField>.
Now, I do not know where exactly it happens, but I do know that Fields (e.g. volScalarField) should register themselves in the objectRegistry. This happens by the regIOobject class, which is a (distant) parent of volScalarField.
The "mesh" is the objectRegistry. That is, no objectRegistry is passed in the constructor of the p-field, hence it cannot register itself.

Therefore, you cannot obtain p.

Two different workarounds:

1. Change createFields.H of interFoam such that p is constructed in the same way as p_rgh. Then I reckon it may be accessed in exactly the same way as you'd access p_rgh
2. Construct p yourself when you need it:
   

   volScalarField p
   (
       IOobject
       (
           "p",
           runTime.timeName(),
           mesh,
           IOobject::NO_READ,
           IOobject::NO_WRITE,
           false
       ),
       p_rgh + rho*gh
   );

   Evidently, this piece of code will need to be adapted as to obtain runTime, mesh, p_rgh, rho and gh correctly. Give it a try.

[alexeym]

Hi all,

@nuttita

calcNu method is called during:
1. construction of incompressibleTwoPhaseMixture.
2. mixture.correct() call.

In the first case pressure field is not available yet, as it is created after mixture object in createFields.H. But it is available on the first mixture.correct() call (just before construction of velocity equation).

So to avoid error, you can check if pressure field is available using foundObject method (http://cpp.openfoam.org/v4/a01730.ht...85a5590fbebd84). If pressure is not there yet, return 0; otherwise use it.

[nuttita]

Hi Kevin and Alexey, thank you for your quick replies.

Kevin, you're right I should construct p the same way as p_rgh in createFields.H. However, the location to construct p is also important. As Alexey mentioned calcNu method is call during construction of incompressibleTwoPhaseMixture. So I need to construct p before that call. Here is my solution.

OLD createFields.H

code:

Info<< "Reading transportProperties\n" << endl;
immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);

volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());

const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();


// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2,
alpha1.boundaryField().types()
);
rho.oldTime();


// Mass flux
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);


// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, mixture)
);


#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"


volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);

my createFields.H

code:

volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

#include "createPhi.H"


Info<< "Reading transportProperties\n" << endl;
immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);

volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());

const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();


// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2,
alpha1.boundaryField().types()
);
rho.oldTime();


// Mass flux
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);


// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, mixture)
);


#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"



p = p_rgh + rho*gh;

Now, it works and I can access p. However, this way I need to have p in /0.

Any discussion or better solution are welcome!

Thanks!
Nuttita

[alexeym]

Hi,

@nuttita

Unfortunately you have missed my point. You keep original createFields.H BUT in your model change

const volScalarField& pp=U_.db().objectRegistry::lookupObject<volScalarField>("p");

to

const objectRegistry& db = U_.db();
if (db.foundObject<volScalarField>("p")) {
  const volScalarField& pp = db.lookupObject<volScalarField>("p");
  ...update nu__, etc...
} else {
  ...update nu__ assuming pp == 0, etc...
}

Or you can remove else part and do not update viscosity if there is no pressure.

[nuttita]

Now I see your point. Thanks for your help, alexeym.

[hasansh1986]

Dear everyone,

By taking a look at objectRegistry::lookupObject(const word& name) method implemented in ObjectRegistryTemplates.C you can see this method iterates through registered objects and returns the object with specified name and Type. If there is no match with the specified name and Type with those registered a Fatal Error happens and you will have your run stopped. So it’s essentially important to make sure that your object is registered in db before calling the lookupObject. For example if you call lookupObject() method before creation and registration of the object you will have Error. And on the other hand by calling the lookupObject() method after registration of the object you won’t have any problem.

const volScalarField& pp = U.db().lookupObject<volScalarField>("p"); //Incorrect
volScalarField p
(
   IOobject
   (
       "p",
       runTime.timeName(),
       mesh,
       IOobject::NO_READ,
       IOobject::AUTO_WRITE,       
   ),
   p_rgh + rho*gh
);
const volScalarField& pp = U.db().lookupObject<volScalarField>("p"); //Correct

volScalarField and other types such as volVectorField are GeometricFields which are inherited from regIOObject. In the constructor of these fields you pass an IOobject which is responsible for registration of the object. The registration process (checkIn method) is implemented in the regIOobject class (NOT in the volScalarField). In cases that you don’t want to register your field you can use false for registerObject() property of the IOobject (the default value for this property is true ) as following:

volScalarField p
(
   IOobject
   (
       "p",
       runTime.timeName(),
       mesh,
       IOobject::NO_READ,
       IOobject::AUTO_WRITE,
       false
   ),
   p_rgh + rho*gh
);
const volScalarField& pp = U.db().lookupObject<volScalarField>("p");//Yields error because the object is not registered

In this case you can see that p field will not be registered.

Cheers,

[alexeym]

Dear hasansh1986,

And what if we are not in control of field creation? For example, we are writing, library, which could be used with arbitrary solver?

[hasansh1986]

If existence of the object is crucial in the algorithm you're implementing you need to make sure that the object is already created (for instance you can create the object in constructor of you class). In these cases FATAL ERROR helps you to realise that something is going wrong in your implementation. On the other hand, if the object is optional in order to prevent FATAL ERROR it's a good practice to check if the object is registered or not (if(db.foundObject<volScalarField>(someName))).

[sippanspojk]

Originally Posted by floquation

Two different workarounds:

1. Change createFields.H of interFoam such that p is constructed in the same way as p_rgh. Then I reckon it may be accessed in exactly the same way as you'd access p_rgh
2. Construct p yourself when you need it:
   

   volScalarField p
   (
       IOobject
       (
           "p",
           runTime.timeName(),
           mesh,
           IOobject::NO_READ,
           IOobject::NO_WRITE,
           false
       ),
       p_rgh + rho*gh
   );

   Evidently, this piece of code will need to be adapted as to obtain runTime, mesh, p_rgh, rho and gh correctly. Give it a try.

Kevin, I like your second suggestion of constructing p when I need it. You also mention that I have to adapt my code so I can access mesh, runTime, etc... I am writing a new functionObject where I started off with foamNewFunctionObject, i.e. minimal amount of libraries are included.

How do I make mesh, runTime, etc. accessible within my code? And also, should I define my volScalarField in my .C or .H file?

Thankful for help here!
David

[einstein_zee]

Originally Posted by sippanspojk

Kevin, I like your second suggestion of constructing p when I need it. You also mention that I have to adapt my code so I can access mesh, runTime, etc... I am writing a new functionObject where I started off with foamNewFunctionObject, i.e. minimal amount of libraries are included.

How do I make mesh, runTime, etc. accessible within my code? And also, should I define my volScalarField in my .C or .H file?

Thankful for help here!
David

Hii there,

no idea what you are trying to do . Looking into this https://cpp.openfoam.org/v7/function...8H_source.html lines 131-134 access to Time and polyMesh is guaranteed. You may try this:

const polyMesh& mesh_ = this->mesh_;
const Time& time_ = this->time_;