if function and turbulence->k()

Zack · December 21, 2016, 18:38

Hi All

I am trying to use if function for the non-zero value of turbulence->k():

if (turbulence->k() != 0.0)
{
do ....
}
else
{
do ...
}

What is the proper way of writing that? I tried different form, such as:

if (turbulence->k() Foam::operator!= 0.0)

which after compiling it gave me the below error:

myCreateFields.H:7:21: error: expected ) before Foam
if (turbulence->k() Foam::operator!= 0.0)
^
myCreateFields.H:7:41: error: could not convert Foam::turbulenceModel::k() from Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > to bool
if (turbulence->k() Foam::operator!= 0.0)
^

Thanks.
Zack

floquation · December 22, 2016, 03:27

Please always write code between [CODE]-tags. Without it, it is barely readable.

The function "turbulence->k()" is defined as:

Code:

        //- Return the turbulence kinetic energy
        virtual tmp<volScalarField> k() const = 0;

In other words, it returns a volScalarField, which is a field. It is not a single value, hence you cannot compare it with a single value (!=0.0) like you are trying.

With that being said, your question cannot be answered, unless you explain to us what exactly it is that you want:
Should this if-function be evaluated for each cell of the domain? Do you want to check if k=0 somewhere? Should k=0 everywhere?

Zack · December 22, 2016, 08:53

Thank you Kevin for your reply and help.

I need to check the k value at each single cell of the domain.

floquation · December 22, 2016, 09:02

Very well. I haven't done this before, but I guess the following should work:

Code:

            for (int celli = 0; celli < mesh.nCells(); celli++){
                if (turbulence->k()[celli] != 0)
                {
                    // Stuff.
                }
                else
                {
                    // Other stuff.
                }                
            }

Code:

            forAll(mesh.cells(), celli){
                if (turbulence->k()[celli] != 0)
                {
                    // Stuff.
                }
                else
                {
                    // Other stuff.
                }
            }

Zack · December 22, 2016, 10:02

Kevin

For both cases I am getting the same error:

Code:

In file included from FSCModifiedbEqn.H:203:0,
                 from FSCModifiedFoam.C:108:
newStModel.H:4:20: error: no match for operator[] (operand types are Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > and int)
 if (turbulence->k()[celli] != 0)
                    ^

floquation · December 22, 2016, 11:29

Hm. I think that is because the k()-function returns a tmp<volScalarField>, which is a "smart-pointer". We must first dereference it to obtain the volScalarField, of which we would like to take the celli'th element. (Again, I haven't done this before, I'm learning these things right now.)

Try the following (let me know whether they work):

Code:

turbulence->k()()[celli] // get const reference

Code:

turbulence->k().ref()[celli] // get non-const reference

See also:

Quote:

${FOAM_SRC}/OpenFOAM/memory/tmp/tmp.H

Zack · December 24, 2016, 17:45

Kevin

Using both of the suggestions, it seems they both can solve the issue with the if-statement. But It seems this causes another issue!
Imagine that in my bEqn.H file I need to solve:

Code:

 fvm::laplacian( rho*Dtt, b)

and I have Dtt which is defined as volScalarField in "myCreateField.H" file.
Before adding the if-statement, the bEqn can recognize the Dtt and I can compile it with no error! But now, after adding if-statement, using both suggestions, I am getting below error:

Code:

 Dtt was not declared in this scope

floquation · January 4, 2017, 16:49

Sorry for being slow -- holidays.

Based on your post, I'm not sure what could possibly go wrong. If a variable is defined outside the if-condition, it is automatically defined inside the if-condition.
Could you post the code that you tried?

Zack · January 6, 2017, 16:09

myCreateFields.H

Code:

Info<<"Reading myCreateFields.H File"<<endl;

 
scalar Cmu=0.09;
volScalarField up = uPrimeCoef*sqrt((2./3.)*turbulence->k());

forAll(mesh.cells(), cellI)
{
if (turbulence->k()()[cellI] != scalar(0.0)){
	#include "Fields1.H"
}
else{
	#include "Fields2.H"
}
}


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

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

//calculating T_Tild
Info<< "Calculating T_Tild\n" << endl;
volScalarField T_Tild
(
	IOobject
	(
		"T_Tild",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	Tu*(b+((thermo.rhou()/rho)*(1-b)))
);

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

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

Fields1.H

Code:

volScalarField L = Foam::pow(Cmu,scalar(0.75))*pow(turbulence->k(),1.5)/turbulence->epsilon();

volScalarField tauTurb = L/up;
volScalarField tauChem = thermo.alpha()/(rho*pow(unstrainedLaminarFlameSpeed()(),2));
volScalarField Da = tauTurb/tauChem ;
volScalarField tO=0.1*L/up;

//Required to calculate Dt and Dtt
scalar PrandtlTurb=0.7;


//calculating St
Info<< "Calculating S_turb\n" << endl;
volScalarField S_turb
(
	IOobject
	(
		"S_turb",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	XiShapeCoef*up*pow(Da,0.25)
);

//calculating Dt
Info<< "Calculating  Dt\n" << endl;
volScalarField Dt
(
	IOobject
	(
		"Dt",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	(Cmu/PrandtlTurb)*(pow(turbulence->k(),2.)/turbulence->epsilon())
);

volScalarField tauTurbPrime = Dt/(Foam::pow(up,2.0));


//calculating S_turbTimeDependent

Info<< "Calculating S_turbTimeDependent \n" << endl;
volScalarField S_turbTimeDependent
(
	IOobject
	(
		"S_turbTimeDependent",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	
	S_turb*Foam::sqrt(mag(1.+tauTurbPrime/(runTime-tO)*(Foam::exp(-1.*(runTime-tO)/tauTurbPrime)-1.)))
	
);


//calculating Dtt
Info<< "Calculating Dtt\n" << endl;
volScalarField Dtt
(
	IOobject
	(
		"Dtt",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	
	Dt * (1-Foam::exp(-1.*(runTime-tO)/tauTurbPrime))
	//Dt * (1-exp((-(runTime+tO))/tauTurbPrime))
);

Fields2.H

Code:

volScalarField L = Foam::pow(Cmu,scalar(0.75))*pow(turbulence->k(),1.5)*turbulence->epsilon();
volScalarField tauTurb = L*up;
volScalarField tauChem = thermo.alpha()/(rho*pow(unstrainedLaminarFlameSpeed()(),2));
volScalarField Da = tauTurb/tauChem ;
volScalarField tO = 0.1*L*up;

//Required to calculate Dt and Dtt
scalar PrandtlTurb=0.7;

//calculating St
Info<< "Calculating  S_turb\n" << endl;
volScalarField S_turb
(
	IOobject
	(
		"S_turb",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	XiShapeCoef*up*pow(Da,0.25)
);

//calculating Dt
Info<< "Calculating Dt\n" << endl;
volScalarField Dt
(
	IOobject
	(
		"Dt",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	(Cmu/PrandtlTurb)*(pow(turbulence->k(),2.0)*turbulence->epsilon())
);

volScalarField tauTurbPrime = Dt*(pow(up,2.0));

//calculating S_turbTimeDependent
Info<< "Calculating  S_turbTimeDependent \n" << endl;
volScalarField S_turbTimeDependent
(
	IOobject
	(
		"S_turbTimeDependent",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	
	S_turb*up
	
);


//calculating Dtt
Info<< "Calculating  Dtt\n" << endl;
volScalarField Dtt
(
	IOobject
	(
		"Dtt",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	
	Dt *up
	
);

bEqn.H

Code:

....
// Create b equation
    // ~~~~~~~~~~~~~~~~~
    fvScalarMatrix bEqn
    (
        fvm::ddt(rho, b)
       + mvConvection->fvmDiv(phi, b)
       + fvm::div(phiSt, b, "div(phiSt,b)")
       - fvm::Sp(fvc::div(phiSt), b)
       - fvm::laplacian( rho*Dtt, b)
       + fvm::Sp(rho*Foam::exp(-activT/T_Tild)/(tr*(1+Dtt/alphab)), b)

	 //-------End--------//
    );
.....

floquation · January 9, 2017, 03:38

Your code may be reduced to the following statement:

Code:

if ( condition )
{
    int a = 0;
}
else{
    int a  = 1;
}
a = a + 1; // error, a not declared

Since you define your variable (here,

) inside the if-condition, it is not known outside the if-condition.
A correct statement would be:

Code:

int a;
if ( condition )
{
    a = 0;
}
else{
    a  = 1;
}
a = a + 1; // OK!

In your case, you'd want three include files:

Code:

#include "Fields_init.H"
if (turbulence->k()()[cellI] != scalar(0.0)){
    #include "Fields_set1.H"
}
else{
    #include "Fields_set2.H"
}

Fields_init.H declares the fields that must be known outside the if-condition.
Fields_set#.H initialise those fields by assigning values to them.

(#include "Fields_init.H" should, of course, be outside your for-loop as well, if those variables should be known outside of the scope of your for-loop.)

Zack · January 22, 2017, 19:16

Hi Kevin

Sorry for late reply and thanks for your help. It worked fine...

One ore question for you. If I need to set the value of any volume scalar field to zero and I consider the unit for that as well, how can I do that?

Lets say that I have VolScalarField B which has (m/s) as its unit and I need to define its value at special condition as zero:

Code:

volScalarField B1 = scalar (0.0);  //unit m/s

floquation · January 23, 2017, 04:51

Code:

    // Just to get myself a volScalarField, you won't need this line:
    volScalarField alpha = this->lookupObject<volScalarField>("alpha.water");

    // Overwrite the entire volScalarField with a constant value of specified dimension:
    alpha=dimensioned<scalar>("name",dimensionSet(0,1,-1,0,0,0,0),scalar(0.0));
    // This will result in a run-time error, because the dimensions of alpha are not m/s, meaning it works.

I urge you to install an IDE. An IDE will allow you to fairly quickly figure out how to do such things by simply looking at the constructors and "operator=" of the entity you wish to create (here: volScalarField). Besides, it has auto-completion and find-declaration functionality to quickly find the relevant pieces of source code.
Personally, I use Eclipse.

December 21, 2016, 18:38	if function and turbulence->k()	#1
Zack New Member Join Date: Sep 2012 Posts: 23 Rep Power: 13	Hi All I am trying to use if function for the non-zero value of turbulence->k(): if (turbulence->k() != 0.0) { do .... } else { do ... } What is the proper way of writing that? I tried different form, such as: if (turbulence->k() Foam::operator!= 0.0) which after compiling it gave me the below error: myCreateFields.H:7:21: error: expected ) before Foam if (turbulence->k() Foam::operator!= 0.0) ^ myCreateFields.H:7:41: error: could not convert Foam::turbulenceModel::k() from Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > to bool if (turbulence->k() Foam::operator!= 0.0) ^ Thanks. Zack

December 22, 2016, 03:27		#2
floquation Senior Member Kevin van As Join Date: Sep 2014 Location: TU Delft, The Netherlands Posts: 252 Rep Power: 20	Please always write code between [CODE]-tags. Without it, it is barely readable. The function "turbulence->k()" is defined as: Code: //- Return the turbulence kinetic energy virtual tmp<volScalarField> k() const = 0; In other words, it returns a volScalarField, which is a field. It is not a single value, hence you cannot compare it with a single value (!=0.0) like you are trying. With that being said, your question cannot be answered, unless you explain to us what exactly it is that you want: Should this if-function be evaluated for each cell of the domain? Do you want to check if k=0 somewhere? Should k=0 everywhere? Zack likes this.

December 22, 2016, 09:02		#4
floquation Senior Member Kevin van As Join Date: Sep 2014 Location: TU Delft, The Netherlands Posts: 252 Rep Power: 20	Very well. I haven't done this before, but I guess the following should work: Code: for (int celli = 0; celli < mesh.nCells(); celli++){ if (turbulence->k()[celli] != 0) { // Stuff. } else { // Other stuff. } } Code: forAll(mesh.cells(), celli){ if (turbulence->k()[celli] != 0) { // Stuff. } else { // Other stuff. } }

December 24, 2016, 17:45		#7
Zack New Member Join Date: Sep 2012 Posts: 23 Rep Power: 13	Kevin Using both of the suggestions, it seems they both can solve the issue with the if-statement. But It seems this causes another issue! Imagine that in my bEqn.H file I need to solve: Code: fvm::laplacian( rho*Dtt, b) and I have Dtt which is defined as volScalarField in "myCreateField.H" file. Before adding the if-statement, the bEqn can recognize the Dtt and I can compile it with no error! But now, after adding if-statement, using both suggestions, I am getting below error: Code: Dtt was not declared in this scope

January 9, 2017, 03:38		#10
floquation Senior Member Kevin van As Join Date: Sep 2014 Location: TU Delft, The Netherlands Posts: 252 Rep Power: 20	Your code may be reduced to the following statement: Code: if ( condition ) { int a = 0; } else{ int a = 1; } a = a + 1; // error, a not declared Since you define your variable (here, ) inside the if-condition, it is not known outside the if-condition. A correct statement would be: Code: int a; if ( condition ) { a = 0; } else{ a = 1; } a = a + 1; // OK! In your case, you'd want three include files: Code: #include "Fields_init.H" if (turbulence->k()()[cellI] != scalar(0.0)){ #include "Fields_set1.H" } else{ #include "Fields_set2.H" } Fields_init.H declares the fields that must be known outside the if-condition. Fields_set#.H initialise those fields by assigning values to them. (#include "Fields_init.H" should, of course, be outside your for-loop as well, if those variables should be known outside of the scope of your for-loop.) Zack likes this.

December 22, 2016, 08:53		#3
Zack New Member Join Date: Sep 2012 Posts: 23 Rep Power: 13	Thank you Kevin for your reply and help. I need to check the k value at each single cell of the domain.

January 4, 2017, 16:49		#8
floquation Senior Member Kevin van As Join Date: Sep 2014 Location: TU Delft, The Netherlands Posts: 252 Rep Power: 20	Sorry for being slow -- holidays. Based on your post, I'm not sure what could possibly go wrong. If a variable is defined outside the if-condition, it is automatically defined inside the if-condition. Could you post the code that you tried?

January 22, 2017, 19:16		#11
Zack New Member Join Date: Sep 2012 Posts: 23 Rep Power: 13	Hi Kevin Sorry for late reply and thanks for your help. It worked fine... One ore question for you. If I need to set the value of any volume scalar field to zero and I consider the unit for that as well, how can I do that? Lets say that I have VolScalarField B which has (m/s) as its unit and I need to define its value at special condition as zero: Code: volScalarField B1 = scalar (0.0); //unit m/s

January 23, 2017, 04:51		#12
floquation Senior Member Kevin van As Join Date: Sep 2014 Location: TU Delft, The Netherlands Posts: 252 Rep Power: 20	Code: // Just to get myself a volScalarField, you won't need this line: volScalarField alpha = this->lookupObject<volScalarField>("alpha.water"); // Overwrite the entire volScalarField with a constant value of specified dimension: alpha=dimensioned<scalar>("name",dimensionSet(0,1,-1,0,0,0,0),scalar(0.0)); // This will result in a run-time error, because the dimensions of alpha are not m/s, meaning it works. I urge you to install an IDE. An IDE will allow you to fairly quickly figure out how to do such things by simply looking at the constructors and "operator=" of the entity you wish to create (here: volScalarField). Besides, it has auto-completion and find-declaration functionality to quickly find the relevant pieces of source code. Personally, I use Eclipse. Zack likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
in k-epsilon wall function approach high Re turbulence models: question of velocity	romant	OpenFOAM Programming & Development	6	May 26, 2016 09:14
issue compiling new turbulence model	perplexed user	OpenFOAM Programming & Development	1	January 13, 2012 03:40
LiencubiclowRemodel	nzy102	OpenFOAM Bugs	14	January 10, 2012 08:53
channelFoam for a 3D pipe	AlmostSurelyRob	OpenFOAM	3	June 24, 2011 13:06
3-D turbulence model with wall function	Henry	Main CFD Forum	7	May 19, 2006 10:40