CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   nuEff()()[patchI] and nuEff() (https://www.cfd-online.com/Forums/openfoam-programming-development/167767-nueff-patchi-nueff.html)

Tobi March 8, 2016 07:16

nuEff()()[patchI] and nuEff()
 
Dear all,

I have a simple question about the following. In the fixedShearStress BC we find this:
Code:

    const label patchI = patch().index();

    const surfaceScalarField& phi =
        db().lookupObject<surfaceScalarField>(phiName_);

    scalarField nuEff;
    if (phi.dimensions() == dimVelocity*dimArea)
    {
        const incompressible::turbulenceModel& turbModel =
            db().lookupObject<incompressible::turbulenceModel>
            (
                "turbulenceModel"
            );

        nuEff = turbModel.nuEff()()[patchI];
    }

As far as I understand it, turbModel.nuEff()()[patchI] should return the values of each face at the patch [patchI]. Therefore the scalarField nuEff should have a size of the nFaces.

At least nuEff is empty:

Code:

nuEff
0()

If I output the nuEff() field like that:
Code:

Info<< turbModel.nuEff() << endl;
I get the list of all boundary's. So far so good. But the access to the patchI is somehow destroyed.

If I use this one, that is exactly the same, I get some complete different values:
Code:

        scalarField nuEff2(turbModel.nuEff()()[patch().index()]);

        Info<< "nuEff2: " << nuEff2 << endl;


Result:

nuEff2
124  (boundary has 500 faces)
(
6.91743e-310
2.3774e-316
4.24399e-314
2.34004e-316
8.3003e-321
4.74303e-322
.
.
.

And if I output it directly I get:
Code:

Info<< "nueff: " << turbModel.nuEff()()[patchI] << endl;


Result:

nueff: 124.601

In a programming point of view this makes not real sense to me (especially the first two):

Code:

        scalarField nuEff2(turbModel.nuEff()()[patch().index()]);
        nuEff = turbModel.nuEff()()[patchI];

Okay there are differences but the field "turbModel.nuEff()()[patchI]" has to be the same as "turbModel.nuEff()()[patch().index()]".

At least the BC will fail because we use nuEff to divide:
Code:

    operator==(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc)));
Any advice?


Finally, if we output the turbModel.nuEff() we get the correct values:
Code:

.
.
.

    movingWall
    {
        type            calculated;
        value          nonuniform List<scalar>
500
(
0.285685
0.285685
0.285685
0.285685
0.285685
0.285685
0.285685
0.285685
.
.
.


I also made a bug report now because it seems really strange. http://openfoam.org/mantisbt/view.php?id=2020

alexeym March 8, 2016 08:22

Hi,

Maybe I am missing something but, turbulenceModel::nuEff() returns tmp<volScalarField>, so nuEff()() (http://foam.sourceforge.net/docs/cpp...2d975c8876d3d6) is volScalarField and nuEff()()[index] should be scalar.

Yet there is turbulenceModel::nuEff(label patchI), which returns viscosity on a patch.

Could you post whole code snippet (so it could be compiled and tested)?

Tobi March 8, 2016 08:41

Hey Bruno,

thanks for the reference about ()() ... here is the link to the file.
If you are correct and ()()[patchI] returns a scalar, why we use a scalarField nuEff than and can you tell me what scalar is returned? Mean value of the patch or sth. like that?

http://www.openfoam.com/documentatio...36_source.html

alexeym March 8, 2016 09:04

Hi,

Not quite Bruno, but OK ;)

Do you talk about line 112 of fixedShearStressFvPatchVectorField.C? It uses turbulenceModel::nuEff(label patchI), which returns tmp<scalarField> (http://www.openfoam.com/documentatio...4f85f6644d2d6a).

Tobi March 8, 2016 09:14

1 Attachment(s)
Aaa its version 3.0 Here it the file I was talking about...

Tobi March 8, 2016 09:17

Quote:

Originally Posted by alexeym (Post 588659)
Hi,

Not quite Bruno, but OK ;)

Do you talk about line 112 of fixedShearStressFvPatchVectorField.C? It uses turbulenceModel::nuEff(label patchI), which returns tmp<scalarField> (http://www.openfoam.com/documentatio...4f85f6644d2d6a).

As I read
Code:

Return the effective viscosity on patch.
it is that i have a scalarField with the size of faces and each field stores one face value, am I right?

alexeym March 8, 2016 09:31

You have scalarField with the size of the patch and it is has scalar for each face value.

Though BC uses this - turbModel.nuEff()()[patchI] - strange construction, turbulenceModel class has nuEff(label patchI) method even in 2.3.1.

Tobi March 8, 2016 09:56

Hi,

Henry changed it: https://github.com/OpenFOAM/OpenFOAM...1fe57e8a140c7c

For 3.0.x it is working now but for 2.3.1 I get an error if I use the code you (or henry) mentioned:
Code:

fixedShearStress/fixedShearStressFvPatchVectorField.C: In member function 'virtual void Foam::fixedShearStressFvPatchVectorField::updateCoeffs()':
fixedShearStress/fixedShearStressFvPatchVectorField.C:136:55: error: no matching function for call to 'Foam::incompressible::turbulenceModel::nuEff(const label&) const'
      const scalarField& nuEff2(turbModel.nuEff(patchI));
                                                      ^
fixedShearStress/fixedShearStressFvPatchVectorField.C:136:55: note: candidate is:
In file included from fixedShearStress/fixedShearStressFvPatchVectorField.C:31:0:
/home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H:199:37: note: virtual Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::incompressible::turbulenceModel::nuEff() const
        virtual tmp<volScalarField> nuEff() const = 0;
                                    ^
/home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H:199:37: note:  candidate expects 0 arguments, 1 provided
fixedShearStress/fixedShearStressFvPatchVectorField.C:136:26: warning: unused variable 'nuEff2' [-Wunused-variable]
      const scalarField& nuEff2(turbModel.nuEff(patchI));
                          ^

So at last, there is no method for nuEff(label)

alexeym March 8, 2016 10:47

Hi,

Even in 2.3.1 there are two types of turbulence models ones in src/turbulenceModels and another in src/TurbulenceModels (surprise, surprise). First type does not have nuEff(const label), second type has it and for RASModel it is implemented as:

Code:

        //- Return the effective viscosity on patch
        virtual tmp<scalarField> nuEff(const label patchi) const
        {
            return this->nut(patchi) + this->nu(patchi);
        }

nu(patchi) is implemented as (example from Newtonian.H):

Code:

        //- Return the laminar viscosity for patch
        tmp<scalarField> nu(const label patchi) const
        {
            return nu_.boundaryField()[patchi];
        }

nut(patchi) is implemented as (example from eddyViscosity.H):

Code:

        //- Return the turbulence viscosity on patch
        virtual tmp<scalarField> nut(const label patchi) const
        {
            return nut_.boundaryField()[patchi];
        }

so your (well, OpenFOAM's) quite frightening

Code:

nuEff = turbModel.nuEff()()[patchI]
could be in fact

Code:

scalarField nuEff(turbModel.nuEff().boundaryField()[patchI]);
Or you can use transportModel's ability to return nu for a given patch using

Code:

transport().nu(patchI)
so it would be

Code:

scalarField nuEff(transport().nu(patchI) + nut().boundaryField()[patchI]);

Tobi March 8, 2016 11:58

Hi,

I will try it tomorrow but today I tried it with boundaryField() without success. Maybe I made a simple mistake. Thanks in advance.

chriss85 March 11, 2016 07:08

I think there might also be problems with memory management. You receive a tmp<volScalarField> without storing it, and only access the boundary. I had a similar issue before and for me it helped to initialize a copy of the volScalarField before accessing the patch values like this:

volScalarField copy(turb->nuEff());

Afterwards you can use the boundaries of the copied field without worrying about the scope of the classes contained in the temporary object. In my case I even regarded it as a bug or atleast a design error.

Tobi March 11, 2016 07:29

Hey,

thanks for the hint. At last the code has to be:
Code:

nuEff = turbModel.nuEff()().boundaryField()[patchI];
I think I will also make a copy (:

thomasArk47 March 11, 2016 18:16

strange things explained
 
Hello Tobi.
First thanks for your tuto and videos. Good job which could help some people and help OF to be more widely used than it is:)

For your problem described here, Alexeym gave you good answer (using turbModel.nuEf()().boundaryField()[patchI]).

To complete the post, I can give you some answers to the "strange" things you have mentioned in your first post. All things are related to the fact that turbModel.nuEff()()[patchI] (which is clearly a bug) is a scalar as mentioned by alexeym.

A) What happens is the following:
scalarField nuEff
--> create a scalarField of size 0
nuEff = turbModel.nuEff()()[patchI]
--> it is the UList<Type>::operator =(Type& v) (with Type=scalar) which assign all values of the List to v. But do nothing here since nuEff is of size 0! So nuEff is still empty list and Info << nuEff gives 0().

B) Other thing: your "scalarField nuEff2(turbmodel.nuEff()()[patchI])"...
In fact, you try to construct a scalarField with a constructor scalarField (scalar& v)
But this constructor doesn't exist...
You have scalarField(label i, scalar v) (build a scalarField of size i and all values equal v) but clearly it is not this one you use. In fact, you use the base constructor scalarField(label i) which construct a scalarField of size i and no prescribed values (so compiler dependant / something like VSMALL for most compilers). And what is your i? The integer part of the scalar turbModel.nuEff()()[patchI] (124 for your example).

These things explain what you have seen.

Tobi March 11, 2016 19:32

Dear Thomas,

thanks for completing the topic and the feedback to my webpage. I am familiar with c++ but the error reported here was from a collegue in the german openfoam forum, that I just investigated during other investigations (: So at last I wanted to report you the bug, learn new things (thanks to you and Alexeym) but I did not investigate into the code deeply and of course I am not a FOAM - code expert (what function will return what or what is doing that). Its good that there are people like you that will clear the questions.

So just one question. If I will create my scalarField nuEff with the size of the real field (lets say 200), and use
Code:

nuEff = turbModel.nuEff()()[patchI]
will it work? I could test it but its 01:32 a.m. and I am too tierd (:

Nice dreams,

thomasArk47 March 12, 2016 15:30

Well if you write:

"
scalarField nuEff(200) // create a scalarField with size 200 (not yet inialized here)
....
nuEff = turbModel.nuEff()()[patchI] // assign all values of the scalarField as equal to the scalar on the right; comes from UList::operator =(scalar &v); doesn't change the size of the container as said in th previous post
"

at end, all the 200 values of the nuEff list will be equal to the same scalar value turbModel.nuEff()()[patchI]. It takes sense from informatic point of view but it is clearly not what you want I think. You should prefer the scalarField nuEff being equal to the scalarField turbModel.nuEff()().boundaryField()[patchI]. For this, you must use alexeym correct expression:

"
scalarField nuEff // create zero sized List (size doesn't matter here -- see later)
nuEff = turbModel.nuEff()().boundaryField()[patchI] // it is the UList::operator =(UList& f) function. This function resizes the calling list (nuEff here) to the size of the called one (f in the prototype / turbModel..... in the example -- this is why the instanciation of nuEff as a zero sized list doesn't matter) and equals the values.
"

gaza October 26, 2017 08:10

Quote:

Originally Posted by thomasArk47 (Post 589349)
Well if you write:

"
scalarField nuEff(200) // create a scalarField with size 200 (not yet inialized here)
....
nuEff = turbModel.nuEff()()[patchI] // assign all values of the scalarField as equal to the scalar on the right; comes from UList::operator =(scalar &v); doesn't change the size of the container as said in th previous post
"

at end, all the 200 values of the nuEff list will be equal to the same scalar value turbModel.nuEff()()[patchI]. It takes sense from informatic point of view but it is clearly not what you want I think. You should prefer the scalarField nuEff being equal to the scalarField turbModel.nuEff()().boundaryField()[patchI]. For this, you must use alexeym correct expression:

"
scalarField nuEff // create zero sized List (size doesn't matter here -- see later)
nuEff = turbModel.nuEff()().boundaryField()[patchI] // it is the UList::operator =(UList& f) function. This function resizes the calling list (nuEff here) to the size of the called one (f in the prototype / turbModel..... in the example -- this is why the instanciation of nuEff as a zero sized list doesn't matter) and equals the values.
"

Hi
I have a problem with nuEff() because as you said it creates list with size zero.

How to force nuEff() to create a list of mesh size?

gaza October 26, 2017 08:45

Basically I overriden nu() to return my volScalarField.
And Now I get the error in nuEff() = this->nut() + this->nu()
that nut() returns zero size list and nu() returns the size equal to e.g. 1400
So the sizes are not equal during summation.

So I am wondering how to force nut() to return the same list size as nu()??

gaza October 27, 2017 03:40

Hi
In the RASModel.H we have a function nuEff()
Code:

        //- Return the effective viscosity

        virtual tmp<volScalarField> nuEff() const

        {

            return tmp<volScalarField>

            (

                new volScalarField

                (

                    IOobject::groupName("nuEff", this->U_.group()),

                     this->nut() + this->nu()

                )

            );

        }

I have overriden nu() as
Code:

Foam::tmp<Foam::volScalarField> nu()
{
    return tmp<volScalarField>(dynamicVisc/density);
}

where dynamicVisc is a volScalarField and density is a dimensionedScalar.

When in the UEqn.H I have a code
Code:

turbulence->divDevReff()
nuEff() is run and it sums nut() and nu(). However there is an error because nut() returns zero size list and nu() returns non zero size list (it is equal to the mesh size).

How can I change the code to nut() returns also mesh size list?

gaza November 4, 2017 05:06

Quote:

Originally Posted by gaza (Post 669373)
Hi
In the RASModel.H we have a function nuEff()
Code:

        //- Return the effective viscosity

        virtual tmp<volScalarField> nuEff() const

        {

            return tmp<volScalarField>

            (

                new volScalarField

                (

                    IOobject::groupName("nuEff", this->U_.group()),

                     this->nut() + this->nu()

                )

            );

        }

I have overriden nu() as
Code:

Foam::tmp<Foam::volScalarField> nu()
{
    return tmp<volScalarField>(dynamicVisc/density);
}

where dynamicVisc is a volScalarField and density is a dimensionedScalar.

When in the UEqn.H I have a code
Code:

turbulence->divDevReff()
nuEff() is run and it sums nut() and nu(). However there is an error because nut() returns zero size list and nu() returns non zero size list (it is equal to the mesh size).

How can I change the code to nut() returns also mesh size list?

Hi
I found how to solve the problem.
I had something like this
Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    volScalarField nu
    (
        IOobject
        (
            "nuLocal",
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        eta_/rho_
    );
    return nu;
}

It compiled however during calculation I got errors like:
Code:

Selecting incompressible transport model Helium
#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OSspecific/POSIX/printStack.C:218
#1  Foam::sigSegv::sigHandler(int) at ~/OpenFOAM/OpenFOAM-v1606+/src/OSspecific/POSIX/signals/sigSegv.C:56
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::List<double>::List(Foam::List<double>&, bool) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/List.C:166
#4  Foam::Field<double>::Field(Foam::Field<double>&, bool) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/Field.C:223
#5  Foam::DimensionedField<double, Foam::volMesh>::DimensionedField(Foam::IOobject const&, Foam::DimensionedField<double, Foam::volMesh>&, bool) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/DimensionedField.C:203 (discriminator 2)
#6  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricField(Foam::IOobject const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/GeometricField.C:498
#7  Foam::viscosityModels::Helium::Helium(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ~/OpenFOAM/przemek-v1606+/src/transportModels/incompressible/viscosityModels/Helium/Helium.C:276 (discriminator 86)
#8  Foam::viscosityModel::adddictionaryConstructorToTable<Foam::viscosityModels::Helium>::New(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ~/OpenFOAM/przemek-v1606+/src/transportModels/incompressible/lnInclude/viscosityModel.H:96 (discriminator 2)
#9  Foam::viscosityModel::New(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/transportModels/incompressible/viscosityModels/viscosityModel/viscosityModelNew.C:58
#10  Foam::singlePhaseTransportModel::singlePhaseTransportModel(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/transportModels/incompressible/singlePhaseTransportModel/singlePhaseTransportModel.C:50 (discriminator 12)
#11  ? at ~/OpenFOAM/OpenFOAM-v1606+/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/../buoyantBoussinesqSimpleFoam/readTransportProperties.H:7
#12  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#13  ? at ??:?
Segmentation fault (core dumped)

where Helium is my viscosityModel.

When I changed the nuCalc() function like this:

Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    volScalarField nu
    (
        IOobject
        (
            "nuLocal",
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        eta_/rho_
    );

    return max
    (
        nuMin_,
        min
        (
            nuMax_,
            nu
        )
    );
}

everything works well (at least solver runs). In my code eta_ is a volScalarField and rho_ is a dimensionedScalar.

Can anybody explain why this second solution works and the first doesn't?

Tobi November 5, 2017 14:43

Can you tell us how you access the field nuHelium() in your solver.

gaza November 5, 2017 16:09

Hi Tobi,
I have a class that derives from singlePhaseTransportModel class.
Inside this class is a pointer to viscosity model. This pointer was private and now I changed it into protected to be able to use it. So I have an access to viscosity model and nu is calculated inside helium viscosity model class as shown in my previous post.

Tobi November 5, 2017 16:50

Okay, and how do you access the function? I mean, you should have something like that somewhere:

Code:

volScalarField foo = pointerToModel->nuHelium();
Actually I am interested in that line you should have somewhere (I guess in another form).

gaza November 6, 2017 02:35

yes it is inside singlePhaseTransportModel
Code:

Foam::tmp<Foam::volScalarField>

 Foam::singlePhaseTransportModel::nu() const

 {

    return viscosityModelPtr_->nu();

 }

you can see it in doxygen. I changed only private to protected for viscosityModelPtr_.

Why are you asking?

Tobi November 6, 2017 03:38

Hi, actually I ment something different but it does not matter. I guess I know the answer to your question:

Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    volScalarField nu
    (
        IOobject
        (
            "nuLocal",
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        eta_/rho_
    );
    return nu;
}

This code should be not working based on the fact that the object nu will be destroyed after the return statement. The returned tmp<> pointer points to a non-allocated memory address and thus you should get a run time error. However, your second example is working because you are returning the object returned by the max() function. The max() function will build a tmp<> object that will not be destroyed and thus, the object survives after returning.


To get your first approach running, this should be done:

Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    tmp<volScalarField> nu
    (
        IOobject
        (
            "nuLocal",
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        eta_/rho_
    );
    return nu;
}

Here, you build the tmp object which knows how many pointers point on the object. If you return it, you have to store the returned object in a tmp<> object again, in order to keep the object alive. After the end of calcNu() is reached, the returned tmp<> object survives and you should not get an run time error.

I hope I formulated it in a good way :) ...

This will work:
Code:

tmp<volScalarField> tnu = viscosityPtr_->nuHelium();

//- Do other stuff

This will fail:
Code:

volScalarField nu = viscosityPtr_->nuHelium();

gaza November 6, 2017 07:14

Hi Tobi,
Thank you very much for explanation, now it is clear for me.

Unfortunately I still have an error with nuEff() function.
This function is run in the line
Code:

turbulence->divDevReff(U)
in the UEqn.H. I have this error now
Code:

Courant Number mean: 0 max: 0
deltaT = 9.999e-06
PIMPLE: iteration 1


--> FOAM FATAL ERROR:
 Field<scalar> f1(14000), Field<scalar> f2(0) and Field<scalar> f3(14000)
    for operation f1 = f2 + f3

    From function void Foam::checkFields(const Foam::UList<T>&, const Foam::UList<Key>&, const Foam::UList<Type3>&, const char*) [with Type1 = double; Type2 = double; Type3 = double]
    in file /home/przemek/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/FieldM.H at line 75.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OSspecific/POSIX/printStack.C:218
#1  Foam::error::abort() at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/error.C:246
#2  Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/errorManip.H:85 (discriminator 4)
#3  void Foam::checkFields<double, double, double>(Foam::UList<double> const&, Foam::UList<double> const&, Foam::UList<double> const&, char const*) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/FieldM.H:75
#4  void Foam::add<double, double>(Foam::Field<Foam::typeOfSum<double, double>::type>&, Foam::UList<double> const&, Foam::UList<double> const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/FieldFunctions.C:770
#5  void Foam::add<double, double, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::typeOfSum<double, double>::type, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:961
#6  Foam::tmp<Foam::GeometricField<Foam::typeOfSum<double, double>::type, Foam::fvPatchField, Foam::volMesh> > Foam::operator+<double, double, Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:961 (discriminator 9)
#7  Foam::RASModel<Foam::IncompressibleTurbulenceModel<Foam::transportModel> >::nuEff() const at ~/OpenFOAM/OpenFOAM-v1606+/src/TurbulenceModels/incompressible/../turbulenceModels/lnInclude/RASModel.H:226 (discriminator 8)
#8  Foam::linearViscousStress<Foam::RASModel<Foam::IncompressibleTurbulenceModel<Foam::transportModel> > >::divDevRhoReff(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const at ~/OpenFOAM/OpenFOAM-v1606+/src/TurbulenceModels/incompressible/../turbulenceModels/lnInclude/linearViscousStress.C:101
#9  Foam::IncompressibleTurbulenceModel<Foam::transportModel>::divDevReff(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const at ~/OpenFOAM/OpenFOAM-v1606+/src/TurbulenceModels/incompressible/lnInclude/IncompressibleTurbulenceModel.C:116
#10  ? at ~/OpenFOAM/przemek-v1606+/applications/solvers/heatTransfer/buoyantBoussinesqSuperFluidPimpleFoam/UEqn.H:35
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  ? at ??:?
Aborted (core dumped)

It results from that nuEff = nut() + nu() and nut() is size of zero and I do not know why.

Tobi November 6, 2017 07:30

There are too less information on what you changed. So it is not possible to give any advice, sorry.

gaza November 6, 2017 07:46

Quote:

Originally Posted by gaza (Post 670577)
Hi Tobi,
Thank you very much for explanation, now it is clear for me.

Unfortunately I still have an error with nuEff() function.
This function is run in the line
Code:

turbulence->divDevReff(U)
in the UEqn.H. I have this error now
Code:

Courant Number mean: 0 max: 0
deltaT = 9.999e-06
PIMPLE: iteration 1


--> FOAM FATAL ERROR:
 Field<scalar> f1(14000), Field<scalar> f2(0) and Field<scalar> f3(14000)
    for operation f1 = f2 + f3

    From function void Foam::checkFields(const Foam::UList<T>&, const Foam::UList<Key>&, const Foam::UList<Type3>&, const char*) [with Type1 = double; Type2 = double; Type3 = double]
    in file /home/przemek/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/FieldM.H at line 75.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OSspecific/POSIX/printStack.C:218
#1  Foam::error::abort() at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/error.C:246
#2  Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/errorManip.H:85 (discriminator 4)
#3  void Foam::checkFields<double, double, double>(Foam::UList<double> const&, Foam::UList<double> const&, Foam::UList<double> const&, char const*) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/FieldM.H:75
#4  void Foam::add<double, double>(Foam::Field<Foam::typeOfSum<double, double>::type>&, Foam::UList<double> const&, Foam::UList<double> const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/FieldFunctions.C:770
#5  void Foam::add<double, double, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::typeOfSum<double, double>::type, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:961
#6  Foam::tmp<Foam::GeometricField<Foam::typeOfSum<double, double>::type, Foam::fvPatchField, Foam::volMesh> > Foam::operator+<double, double, Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ~/OpenFOAM/OpenFOAM-v1606+/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:961 (discriminator 9)
#7  Foam::RASModel<Foam::IncompressibleTurbulenceModel<Foam::transportModel> >::nuEff() const at ~/OpenFOAM/OpenFOAM-v1606+/src/TurbulenceModels/incompressible/../turbulenceModels/lnInclude/RASModel.H:226 (discriminator 8)
#8  Foam::linearViscousStress<Foam::RASModel<Foam::IncompressibleTurbulenceModel<Foam::transportModel> > >::divDevRhoReff(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const at ~/OpenFOAM/OpenFOAM-v1606+/src/TurbulenceModels/incompressible/../turbulenceModels/lnInclude/linearViscousStress.C:101
#9  Foam::IncompressibleTurbulenceModel<Foam::transportModel>::divDevReff(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const at ~/OpenFOAM/OpenFOAM-v1606+/src/TurbulenceModels/incompressible/lnInclude/IncompressibleTurbulenceModel.C:116
#10  ? at ~/OpenFOAM/przemek-v1606+/applications/solvers/heatTransfer/buoyantBoussinesqSuperFluidPimpleFoam/UEqn.H:35
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  ? at ??:?
Aborted (core dumped)

It results from that nuEff = nut() + nu() and nut() is size of zero and I do not know why.

Ok I found the origin of the error. I had the code in UEqn.H
Code:

...
    nut = turbulence->nut();

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::ddt(U) + fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)   
    ==
        fvOptions(U)
    );

If I comment the line
Code:

nut = turbulence->nut()
everything works well, but I do not know why???

Tobi November 6, 2017 07:49

How did you define your object nut. Can you give the definition of that object? And where do you use nut afterwards?

gaza November 6, 2017 08:01

nut definition in createFields.H:

Code:

volScalarField nut
(
    IOobject
    (
        "nut",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    turbulence->nut() 
);

and I do not use it before line nut = turbulence->nut();

Tobi November 6, 2017 08:08

This cannot work because turbulence->nut() will return a tmp<> object which will be destoyed in your case :). So build your own field such as:

Code:

volScalarField nut
(
    IOobject
    (
        "nut",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("nut", dimensionSet(...), value(0))
);

then you should get access to the field such as:
Code:


// Keep tmp alive
tmp<volScalarField> tnut = turbulence->nut();

// Get reference to nut object
volScalarField& nutIF = nut;

// Set nut to be tnut from turbulence model
nutIF = tnut();

Or simply use:
Code:

tmp<volScalarField> tnut = turbulence->nut();
without creating a volScalarField in the createFields.H

Zeppo November 6, 2017 16:44

Quote:

Originally Posted by Tobi (Post 670553)
To get your first approach running, this should be done:

Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    tmp<volScalarField> nu
    (
        IOobject
        (
            "nuLocal",
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        eta_/rho_
    );
    return nu;
}


This should be corrected a little bit:
Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    tmp<volScalarField> nu
    (
        new volScalarField
        (
            IOobject
            (
                "nuLocal",
                U_.time().timeName(),
                U_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            eta_/rho_
        );
    )
    return nu;
}

To my understanding it is essential that the object of volScalarField type is to be created in dynamic memory (on the heap) with operator new. If so, tmp can store a pointer to some place in dynamic memory which will stay untouched while tmp is returned from calcNu. On the other hand, when you initialize tmp with a pointer/reference to an automatic object (i.e., in simple words, you don't use operator new to create an object) tmp stores a pointer to a place in automatic memory (on the stack) which will be "erased" as soon as you return from calcNu and it yields undefinite behaviour.

Tobi November 6, 2017 18:10

Quote:

On the other hand, when you initialize tmp with a pointer/reference to an automatic object (i.e., in simple words, you don't use operator new to create an object) tmp stores a pointer to a place in automatic memory (on the stack) which will be "erased" as soon as you return from calcNu and it yields undefinite behaviour.
Is is like that? I don't think so because the tmp<> class knows how much pointer points to the object. Therefore, if one returns the object such as:
Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    tmp<volScalarField> nu
    (
        IOobject
        (
            "nuLocal",
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        eta_/rho_
    );
    return nu;
}

while taking the returned pointer such as:
Code:

tmp<volScalarField> tnut = viscouseModel->calcNu();
The object should be still alive, even if we leave the calcNu() function because a pointer still points onto the object. But out of the box I cannot proof that. Bruno told me once the following:

Quote:

"tmp" is a pointer wrapper class that handles reference counting depending on the operator that is used.

You should see several examples in OpenFOAM, such as:
const tmp<volScalarField> tnuEff = turbulence->nuEff();
const volScalarField& nuEff = tnuEff();
"tnuEff" is a self-destructing capsule that stores a pointer (it's a something something "RAII" trick in C++). When the destructor of "tmp" is called, it will take the pointer inside with it.

If you do this:
turbulence->nuEff()()
It destroys its content after self-destruction, namely when the destructor is called for the object that was returned by "return tmp<volScalarField> nuEff", because no one is holding it after this line of code is processed.

That's why you must keep a live object of "tmp" inside the current method/function.


The "autoptr" class that sometimes is used in C++ for managing pointers does this as well.

Zeppo November 7, 2017 16:12

Any automatic object created within the body of a function can live untill the execution flow leaves the function. Then a destructor of the object is invoked rendering the object basically dead. A pointer to the automatic object can be incapsulated into a smart pointer object and returned from the function, but it will be useless anyway because the object it points to is dead. If you want an object to stay alive when the execution flow leaves the function you should create it with operator new as a dynamically allocated memory object.

Tobi November 7, 2017 17:27

Hmmm,.... actually without the new keyword it is not working at all based on the fact that the pointer has to point to an explicit allocated memory build by new. So you are right. Thanks for correcting me. For those who want to test it:
Code:

Foam::tmp<Foam::volScalarField> myCalc (const Foam::fvMesh& mesh)             
{                                                                             
    tmp<volScalarField> T                                                     
    (                                                                         
        new volScalarField                                                     
        (                                                                     
            IOobject                                                           
            (                                                                 
                "T",                                                           
                mesh.time().timeName(),                                       
                mesh,                                                         
                IOobject::MUST_READ,                                           
                IOobject::AUTO_WRITE                                           
            ),                                                                 
            mesh                                                               
        )                                                                     
    );                                                                         
                                                                               
    return T;                                                                 
}                                                                             
                                                                               
int main(int argc, char *argv[])                                               
{                                                                             
    #include "setRootCase.H"                                                   
    #include "createTime.H"                                                   
    #include "createMesh.H"                                                   
       
    //- Working - keeping object alive                                                                 
    tmp<volScalarField> re = myCalc(mesh);                                                                                                       
    Info<< re() << endl;               

    //- Run-Time Error - object destroyed, reference point to some non-allocated memory
    const volScalarField& foo = myCalc(mesh);                                 
    Info<< foo << endl; 
   
                                       
                                                                               
                                                                               
    Info<< "End\n" << endl;                                                   
                                                                               
    return 0;                                                                 
}


Zeppo November 8, 2017 16:56

Hi, Tobias
Code:

//- Run-Time Error - object destroyed, reference point to some non-allocated memory
    const volScalarField& foo = myCalc(mesh);                                 
    Info<< foo << endl;

Are you sure this piece of code issues an error? I didn't try it personally, but I believe it should be ok. myCalc returns a temporary object which is assigned to a const reference foo. And c++ rules make it possible: a constant reference can be "initialized" with a temporary object. The temporary object doesn't die right away, it's lifetime prolongates up to the constant reference lifetime. And they both will live untill the execution flow leaves the scope where the constant reference was defined. This code might not work for one reason though: myCalc returns tmp and foo expects volScalarField, I am not sure tmp can "converts" to volScalarField implicitly. If you do the conversion explicitly (with operator()) it should work like a breeze:
Code:

const volScalarField& foo =  myCalc(mesh)();

Tobi November 8, 2017 17:01

Hi,

yes it will give an runtime error. If one test the above code, the error appears. I was struggling about that in my openComfort library too because here I need access to the nut field which comes from the turbulence model. Right at the beginning I did this one:
Code:

const volScalarField& nut = turbulence->nut();
However, this code failed always. Therefore, I replaced it with a copy operator which worked (obviously):
Code:

const volScalarField nut = turbulence->nut()
I was not happy about copying the object and asked Bruno. He gave me the explanation, that we have to keep a object of the tmp<> class alive by doing that:
Code:

const tmp<volScalarField> tnut = turbulence->nut()
Otherwise the object returned by the nut() function will be destroyed. You might be right that it will work with smart pointers like the auto_ptr in c++ but with the FOAM tmp<> class actually not (tested it).

Zeppo November 8, 2017 17:09

Hopefully I could test it tomorrow and then get back to this thread.

Zeppo November 9, 2017 16:35

I played around with the code and have to confirm that
Code:

const volScalarField& nut = turbulence->nut()();
fails at runtime. You have to use a pair of round brackets at the end of this expression to obtain (a reference to) volScalarField out of tmp. Otherwise it won't even compile let alone run.

Tobi November 9, 2017 18:47

Hmmm.... :/
Sure, it is obvious now.
After the operator () we have access to the object and can take the reference. My way was wrong. Should think more before I reply.

As always Sergei. Good point. Thanks for sharing.

gaza November 11, 2017 08:45

Quote:

Originally Posted by Zeppo (Post 670639)
This should be corrected a little bit:
Code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Helium::calcNu()
{
    tmp<volScalarField> nu
    (
        new volScalarField
        (
            IOobject
            (
                "nuLocal",
                U_.time().timeName(),
                U_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            eta_/rho_
        );
    )
    return nu;
}

To my understanding it is essential that the object of volScalarField type is to be created in dynamic memory (on the heap) with operator new. If so, tmp can store a pointer to some place in dynamic memory which will stay untouched while tmp is returned from calcNu. On the other hand, when you initialize tmp with a pointer/reference to an automatic object (i.e., in simple words, you don't use operator new to create an object) tmp stores a pointer to a place in automatic memory (on the stack) which will be "erased" as soon as you return from calcNu and it yields undefinite behaviour.

Hi Tobi and Zeppo
Thank you guys for explanation all of these things and help. I solved my problem mentioned in #25.
It was caused by the bed use of tmp object.

I fixed the code according Zeppo solution as in the citation.
Also I added this line according to Tobi:
Code:

tmp<volScalarField> tnut = turbulence->nut();
alphaEff = turbulence->nu() + tnut();

Thank's for sharing :)


All times are GMT -4. The time now is 03:54.