CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

lookupObject for <volScalarField>("p") interFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree4Likes
  • 2 Post By alexeym
  • 1 Post By hasansh1986
  • 1 Post By einstein_zee

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 24, 2017, 21:37
Default lookupObject for <volScalarField>("p") interFoam
  #1
New Member
 
Nuttita Pophet
Join Date: Oct 2016
Location: Mississippi, USA
Posts: 10
Rep Power: 9
nuttita is on a distinguished road
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
nuttita is offline   Reply With Quote

Old   January 25, 2017, 04:39
Default
  #2
Senior Member
 
floquation's Avatar
 
Kevin van As
Join Date: Sep 2014
Location: TU Delft, The Netherlands
Posts: 252
Rep Power: 20
floquation will become famous soon enough
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:
Code:
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:
    Code:
    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.
floquation is offline   Reply With Quote

Old   January 25, 2017, 12:20
Default
  #3
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to 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.
alexeym is offline   Reply With Quote

Old   January 25, 2017, 20:59
Default
  #4
New Member
 
Nuttita Pophet
Join Date: Oct 2016
Location: Mississippi, USA
Posts: 10
Rep Power: 9
nuttita is on a distinguished road
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:

Quote:
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:

Quote:
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
nuttita is offline   Reply With Quote

Old   January 26, 2017, 03:43
Default
  #5
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Hi,

@nuttita

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

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

Code:
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.
atulkjoy and PanPeter like this.
alexeym is offline   Reply With Quote

Old   January 27, 2017, 13:13
Default
  #6
New Member
 
Nuttita Pophet
Join Date: Oct 2016
Location: Mississippi, USA
Posts: 10
Rep Power: 9
nuttita is on a distinguished road
Now I see your point. Thanks for your help, alexeym.
nuttita is offline   Reply With Quote

Old   January 29, 2017, 12:06
Default
  #7
New Member
 
Join Date: May 2016
Posts: 2
Rep Power: 0
hasansh1986 is on a distinguished road
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.


Code:
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:
Code:
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,
CorbinMG likes this.
hasansh1986 is offline   Reply With Quote

Old   January 29, 2017, 12:28
Default
  #8
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to 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?
alexeym is offline   Reply With Quote

Old   January 29, 2017, 13:35
Default
  #9
New Member
 
Join Date: May 2016
Posts: 2
Rep Power: 0
hasansh1986 is on a distinguished road
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))).
hasansh1986 is offline   Reply With Quote

Old   November 19, 2019, 04:29
Default
  #10
Member
 
David Andersson
Join Date: Oct 2019
Posts: 46
Rep Power: 6
sippanspojk is on a distinguished road
Quote:
Originally Posted by floquation View Post
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:
    Code:
    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

Last edited by sippanspojk; November 21, 2019 at 03:21.
sippanspojk is offline   Reply With Quote

Old   November 19, 2019, 10:59
Default
  #11
Member
 
Hosein
Join Date: Nov 2011
Location: Germany
Posts: 93
Rep Power: 14
einstein_zee is on a distinguished road
Quote:
Originally Posted by sippanspojk View Post
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:
Code:
const polyMesh& mesh_ = this->mesh_;
const Time& time_ = this->time_;
Zhiheng Wang likes this.
einstein_zee is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
InterFoam stops after deltaT goes to 1e14 francesco_b OpenFOAM Running, Solving & CFD 9 July 25, 2020 07:36
interFoam vs. simpleFoam channel flow comparison DanM OpenFOAM Running, Solving & CFD 12 January 31, 2020 16:26
interFoam (HELYX-OS) pressure boundary conditions SFr OpenFOAM Running, Solving & CFD 8 June 23, 2016 17:36
k-e & GAMG interFoam Schemitisation Stability Issue JFM OpenFOAM Running, Solving & CFD 3 December 1, 2015 06:58
interFoam in parallel gooya_kabir OpenFOAM Running, Solving & CFD 0 December 9, 2013 06:09


All times are GMT -4. The time now is 00:33.