CFD Online Discussion Forums - Creating a temporary volScalarField in the calculation

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Creating a temporary volScalarField in the calculation (https://www.cfd-online.com/Forums/openfoam-programming-development/183265-creating-temporary-volscalarfield-calculation.html)

Creating a temporary volScalarField in the calculation

I'm trying to create a new nonNewtonian model and in my .C file, I have a private member function to calculate nu:

Code:

Foam::tmp<Foam::volScalarField>

Foam::viscosityModels::ShearThinning::calcNu() const

{



forAll(strainRate()(),cellI)

{

if ((strainRate()()[cellI] > 0.17) && (strainRate()()[cellI] < 700))

{

    return

    k_*((-0.00795*log(strainRate()/m_) + 0.344)/(sqr(log(strainRate()/m_)) + 9.02*log(strainRate()/m_) + 25.0) + (-0.0128*log(strainRate()/m_) + 0.398)/(sqr(log(strainRate()/m_)) + 6.11*log(strainRate()/m_) + 14.5))/rho_;

}

else

{

    return nu0;

}

}



}

where nu0 is a COSTANT value. I know to make the code work, I need to define nu0 as a volScalarField as the type of calcNu() is volScalarField. I want to create nu0 as a temperary volScalarField only for this calculation. Can anyone suggest how can I do it?

Thanks,

Dear Shuang,

I would not do it as you suggested because if you make a temporary volScalarField lead to memory allocation and deallocation in each time-step which is not a good way. It slows down your simulation and is waste of power. I think nowadays you should not have problems with memory space, so there is no need doing that. Back to your question. I think you know that your code cannot work because of not defined nu0 variable. However, please use code tags to visualize the c++ stuff. To implement a temporary volScalarField nu0 here you have to do it like that:

Code:



//- Copy-Constructor (but the boundaries conditions are similar to k_'s one)

volScalarField nu0 = k_;



//- And then re-initialize nu0 or do what you want. Changing BC, values etc.



//- Using normal constructor

volScalarField nu0

(

    IOobject

    (

        "nu0",

        .

        .

        .

        .

     ),

     k_.mesh()

     //- Maybe boundary conditions

     .

     .

     .

);

Since you’re trying to create a new Non-Newtonian model your class should be inherited from viscosityModel class. So, have access to U_ field which can be used for creating your volScalarField in return of your method.

Code:

// Your If statement

else

{

       return  tmp<volScalarField>

               (

                   new volScalarField

                   (

                       IOobject

                       (

                           "nu0",

                           U_.time().timeName(),

                           U_.db(),

                           IOobject::NO_READ,

                           IOobject::NO_WRITE,

                           false

                       ),

                       U_.mesh(),

                       dimensionedScalar("nu0", dimViscosity*dimDensity, nu0)

                   )

                );

}

Where nu0 is constant.

Cheers,

Thank you Hasan. This is helpful

Hi Hasan,

I have another practical question, I hope you can help out.

I create a temporary nu0 field using your approach and initialise the elements value to be 0. What I want to do right now is to update each element value based on the filter I applied. My code is below:

Quote:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::ShearThinning::calcNu() const
{
return
tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"nu0",
U_.time().timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
U_.mesh(),
dimensionedScalar("nu0",dimensionSet(0,2,-1,0,0,0,0),0)
)
);

forAll(nu0,cellI)
{
if (strainRate()()[cellI] > 0.17)
{
nu0[cellI] = strainRate()()[cellI]*0.00001;
}
else
{
nu0[cellI] = strainRate()()[cellI]*0.000004;
}
}
}

However, when I compile it, I got the following error.

Quote:

ShearThinning.C: In member function ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::viscosityModels::ShearThinning::calcNu() const’:
ShearThinning.C:87:1: error: ‘nu0’ was not declared in this scope

My question is
1. why do I need to declare nu0 again given I just created?
2. how can I declare then?

I hope you can help out.

Thanks,

Since you haven't defined any variable named nu0_, you're getting compile error. If you want to perform some operation on nu0_ you need to define it using its constructor.

Code:

volScalarField nu0_

                   (

                       IOobject

                       (

                           "nu0_",

                           U_.time().timeName(),

                           U_.db(),

                           IOobject::NO_READ,

                           IOobject::NO_WRITE,

                           false

                       ),

                       U_.mesh(),

                       dimensionedScalar("nu0_", dimViscosity, nu0)

                   )



                );

Then perform your operations on nu0_ and finally return it. Keep in mind shouldn't write any code after return of your method.

Hi Hasan,

So I declared the nu0 in my .H file and my .C code is as below:

Quote:

// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::ShearThinning::calcNu()
{

return
tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"nu0",
U_.time().timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
U_.mesh(),
dimensionedScalar("nu0",dimensionSet(0,2,-1,0,0,0,0),0)
)
);

forAll(nu0(),cellI)
{
if ((strainRate()()[cellI] > 0.17)&&(strainRate()()[cellI] < 700))
{
nu0()[cellI] = 0.5*((-0.00795*log(strainRate()()[cellI]) + 0.344)/(sqr(log(strainRate()()[cellI])) + 9.02*log(strainRate()()[cellI]) + 25.0) + (-0.0128*log(strainRate()()[cellI]) + 0.398)/(sqr(log(strainRate()()[cellI])) + 6.11*log(strainRate()()[cellI]) + 14.5))/1060;
}
if (strainRate()()[cellI] < 0.17)
{
nu0()[cellI] = 0.046/1060;
}
if (strainRate()()[cellI] > 700)
{
nu0()[cellI] = 0.004/1060;
}
}
}

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

Foam::viscosityModels::ShearThinning::ShearThinnin g
(
const word& name,
const dictionary& viscosityProperties,
const volVectorField& U,
const surfaceScalarField& phi
)
:
viscosityModel(name, viscosityProperties, U, phi),
ShearThinningCoeffs_(viscosityProperties.subDict(t ypeName + "Coeffs")),

nu_
(
IOobject
(
"nu",
U_.time().timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
calcNu()
)
{}

So my idea is, create a temporary volFieldScalar field nu0 that is editable and accessible and then use for loop to edit the elements in nu0. nu0 will then return to calcNu(). Actual nu_ is constructed in the constructors using the value of calcNu().

There is no compile error any more but when I run the model. It turns out the program did not go through the for loop, so at any time step, the nu is the constant initialised in nu0.

Is there anything wrong in my logic? How could I modify it?

Thanks

Hi,

I did not check everything (reason below) but your function calcNu() will just return a volScalarField which is always zero. You should know that the return statement is the last execution in the function. What does that mean? It means that your forAll loop is senseless because you will never go to that point. Ergo, it will not work as you would like to have. By the way, I am sorry not to see that you want to build your own non-Newton model. Otherwise I would told you the same than Hasan (derive it from the non-Newtonian class) in order to have access to fields. However, I have the feeling that you are not aware with C++ and finally FOAM is nothing more than C++. To sum up:

Code:



//- Function that takes one argument (const scalar) and return one scalar

scalar myFunction(const scalar foo)

{

    //- What ever you want to do

    scalar temp = foo*sqrt(2);

    .

    .

    .



    //- The last operation in your function

    return temp;



    //- Everything that come after the return keyword will not be executed

    Info<< "Still in myFunction()??? \n" << endl;

}

At last one comment. Please ensure that your code has a nice format. There are no shiftings and therefore, it is not nice to read. Maybe I am the only one who do not like it but as far as I reply to a lot of threads over the years, it shows me - somehow - that the thread starter (or the one who ask) does not take time for the question and demonstrate the problem in a nice and clear way. On the contrary: But the people who give support (or wants to give) should spend time to your problem?

Quote:

Originally Posted by Tobi (Post 636196)

Please ensure that your code has a nice format. There are no shiftings and therefore, it is not nice to read. Maybe I am the only one who do not like it but (...)

I'm fairly sure that applies to everyone.
Unless I have a very good mood, I have a tendency to skip those topics altogether.

There is a post about using [CODE]-tags here, but I'm afraid that it being post #6 instead of #1, people simply do not see it. Nonetheless, it should be common forum sense... If you search before asking your question [which is mandatory], you'll also see others use it. Therefore, if someone does not use it, it is rather likely that that person did not search in the first place.

Help on temporary volScalarField

Dear Foamers, I found this conversation suitable for a question regarding a volScalarField private member function.
I'm updating an OpenFOAM k-epsilon code from the 2.3.0 version to the 4.0.
In doing so, when compiling my .C file, I get an error in correspondence of these lines:

Code:

tmp<volScalarField> kEpsilon1::gCmuSmooth()

{

   volScalarField Cmus=pow(ustar_,4)/sqr(kref_);  

    // Return with smoothed value

    return min(Cmus,CmuMax_);  // Here Cmu_ changed to Cmus_

}

where gCmuSmooth was declared in the .H file as:

Code:

{

tmp<volScalarField> gCmuSmooth();

}

which states: " 'template<class BasicTurbulenceModel> class Foam::RASModels::K=kEpsilon1' used without template parameters
tmp<volScalarField> gCmuSmooth() "

This part of the code, written in the same way, was working perfectly on the 2.3.0 version of OpenFOAM. Should I change something in the way I declare or initialise gCmuSmooth() ?
Any help or explanation would be extremely useful!

Quote:

Originally Posted by hasan_shetabivash (Post 635905)

Code:

// Your If statement

else

{

       return  tmp<volScalarField>

               (

                   new volScalarField

                   (

                       IOobject

                       (

                           "nu0",

                           U_.time().timeName(),

                           U_.db(),

                           IOobject::NO_READ,

                           IOobject::NO_WRITE,

                           false

                       ),

                       U_.mesh(),

                       dimensionedScalar("nu0", dimViscosity*dimDensity, nu0)

                   )

                );

}

Where nu0 is constant.

Cheers,

Hi, By Using this code i am trying to find min of volScalarField.
volScalarField MINprincipalstress
(
IOobject
(
"MINprincipalstress",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
min(principalstress, MINprincipalstress)

);
I have calculated the principal stress by using this code:
volScalarField principalstress
(
IOobject
(
"principalstress",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
(0.5 * (sigmaxx + sigmayy) + 0.5 * sqrt( sqr(sigmaxx- sigmayy) + 4 * sqr(sigmaxy)))

);

When i compile the code, it does not gives me any error. But when i run the simulation. after some iteration it gives me this error. I do not know, what should i do now. here is the error:

GAMG: Solving for Dx, Initial residual = 9.14894e-07, Final residual = 9.14894e-07, No Iterations 0
GAMG: Solving for Dy, Initial residual = 9.97011e-07, Final residual = 9.97011e-07, No Iterations 0
ExecutionTime = 6.99 s ClockTime = 7 s

Iteration: 9.994e-06

GAMG: Solving for Dx, Initial residual = 9.14894e-07, Final residual = 9.14894e-07, No Iterations 0
GAMG: Solving for Dy, Initial residual = 9.97011e-07, Final residual = 9.97011e-07, No Iterations 0
ExecutionTime = 6.99 s ClockTime = 7 s

Iteration: 9.996e-06

GAMG: Solving for Dx, Initial residual = 9.14894e-07, Final residual = 9.14894e-07, No Iterations 0
GAMG: Solving for Dy, Initial residual = 9.97011e-07, Final residual = 9.97011e-07, No Iterations 0
ExecutionTime = 6.99 s ClockTime = 7 s

Iteration: 9.998e-06

GAMG: Solving for Dx, Initial residual = 9.14894e-07, Final residual = 9.14894e-07, No Iterations 0
GAMG: Solving for Dy, Initial residual = 9.97011e-07, Final residual = 9.97011e-07, No Iterations 0
ExecutionTime = 6.99 s ClockTime = 7 s

Iteration: 1e-05

GAMG: Solving for Dx, Initial residual = 9.14894e-07, Final residual = 9.14894e-07, No Iterations 0
GAMG: Solving for Dy, Initial residual = 9.97011e-07, Final residual = 9.97011e-07, No Iterations 0
Max sigmaEq = 1.89073e+09
Max sigmazz = 2.35054e+08
Principal Stress Maximum = 9.42315e+08
Principal Stress Minimum = -4.72528e+08

--> FOAM FATAL ERROR: (openfoam-2312)
Different dimensions for 'min(a, b)'
dimensions : [1 -1 -2 0 0 0 0] != [0 0 0 0 0 0 0]

one more thing if i use std::min(principalstress, MINprincipalstress) it gives me an error this
GAMG: Solving for Dx, Initial residual = 9.14894e-07, Final residual = 9.14894e-07, No Iterations 0
GAMG: Solving for Dy, Initial residual = 9.97011e-07, Final residual = 9.97011e-07, No Iterations 0
ExecutionTime = 7.44 s ClockTime = 8 s

Iteration: 1e-05

GAMG: Solving for Dx, Initial residual = 9.14894e-07, Final residual = 9.14894e-07, No Iterations 0
GAMG: Solving for Dy, Initial residual = 9.97011e-07, Final residual = 9.97011e-07, No Iterations 0
Max sigmaEq = 1.89073e+09
Max sigmazz = 2.35054e+08
Principal Stress Maximum = 9.42315e+08
Principal Stress Minimum = -4.72528e+08

--> FOAM FATAL ERROR: (openfoam-2312)
Cannot dereference nullptr at index 0 in range [0,5)

can you let me know, what is issue here?