CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   If-else condition over the whole mesh (http://www.cfd-online.com/Forums/openfoam-programming-development/79677-if-else-condition-over-whole-mesh.html)

vkrastev August 31, 2010 06:33

If-else condition over the whole mesh
 
Hi everybody, I have some troubles with finding the correct way to modify a source .C file in order to seek for positive/negative scalar values over the whole range of mesh cells. Namely, I'm trying to implement an improved version of the k-omega turbulence model, in which there is a correction of some of the closure coefficients that implies an if-else "check" over a field of scalar values. The following is the part of the source code containing the if-else condition (in the form I trivially suppose it could be correct):

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

tmp<volScalarField> kOmegaSI::kiKappa() const
{
return (1/(pow(omega_,3.0)))*((fvc::grad(k_))&(fvc::grad(ome ga_)));
}

tmp<volScalarField> kOmegaSI::fBetaStar() const
{
{
forAll(mesh.cells(),celli)
if (kiKappa()[celli] > scalar(0.0)) {
return ((scalar(680.0)+sqr(kiKappa()[celli]))/(scalar(400.0)+sqr(kiKappa()[celli])));
}
else {
return scalar(1.0);
}
}
}

And here it is the error message which appears after running the wmake utility:

kOmegaSI_1/kOmegaSI.C:56: error: ‘mesh’ was not declared in this scope
kOmegaSI_1/kOmegaSI.C:57: error: no match for ‘operator[]’ in ‘Foam::incompressible::RASModels::kOmegaSI::kiKapp a() const()[celli]’
kOmegaSI_1/kOmegaSI.C:58: error: no match for ‘operator[]’ in ‘Foam::incompressible::RASModels::kOmegaSI::kiKapp a() const()[celli]’
kOmegaSI_1/kOmegaSI.C:58: error: no match for ‘operator[]’ in ‘Foam::incompressible::RASModels::kOmegaSI::kiKapp a() const()[celli]’
kOmegaSI_1/kOmegaSI.C:61: error: conversion from ‘double’ to non-scalar type ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ requested

So, can anybody give me some hints about my problem? Of course, if necessary, I could add additional information about the .C (and the corresponding .H header) file.

Thank you in advance

ngj August 31, 2010 07:14

Hi

There are several issues in your code.

1. You have a function kiKappa which is called for every single element in the mesh. The function, however, computes the sought for values over the entire mesh. I would move kiKappa into fBetaStar unless kiKappa is needed elsewhere.

2. mesh is termed mesh_ in the turbulence models. See error on line 56.

3. Your return statements are wrong. You return a scalar, however the function fBetaStar is declared as returning a tmp<volScalarField>. So you most loop over a given field, populate it and return it as the correct type.

Best regards,

Niels

vkrastev August 31, 2010 09:02

Quote:

Originally Posted by ngj (Post 273462)
Hi

There are several issues in your code.

1. You have a function kiKappa which is called for every single element in the mesh. The function, however, computes the sought for values over the entire mesh. I would move kiKappa into fBetaStar unless kiKappa is needed elsewhere.

2. mesh is termed mesh_ in the turbulence models. See error on line 56.

3. Your return statements are wrong. You return a scalar, however the function fBetaStar is declared as returning a tmp<volScalarField>. So you most loop over a given field, populate it and return it as the correct type.

Best regards,

Niels

Hi Niels, first of all thanks for the reply. Next, I've tried to modify the code following the first two points that you have remarked and this is what I have obatined:

tmp<volScalarField> kOmegaSI::fBetaStar() const
{
volScalarField kiKappa=(1/(pow(omega_,3.0)))*((fvc::grad(k_))&(fvc::grad(ome ga_)));
{
forAll(mesh_.cells(),celli)
if (kiKappa[celli] > scalar(0.0)) {
return ((scalar(680.0)+sqr(kiKappa[celli]))/(scalar(400.0)+sqr(kiKappa[celli])));
}
else {
return scalar(1.0);
}
}
}

with the corresponding error message:

kOmegaSI_1/kOmegaSI.C: In member function ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::incompressible::RASModels::kOmegaSI::fBetaSt ar() const’:
kOmegaSI_1/kOmegaSI.C:59: error: conversion from ‘double’ to non-scalar type ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ requested
kOmegaSI_1/kOmegaSI.C:62: error: conversion from ‘double’ to non-scalar type ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ requested

So it seems like the first two points have been fixed, but honestly I do not understand what do you mean about "loop over a given field, populate it and return it as the correct type". Namely, I got the concept but I don't know how to practically realize it...if you could give me some practical advice about it I will be really grateful....anyway, thank you once again

Regards

Vesselin

ngj August 31, 2010 11:09

Perhaps you could do something like this:

Code:

tmp<volScalarField> kOmegaSI::fBetaStar() const
{
    volScalarField kiKappa=(1/(pow(omega_,3.0)))*((fvc::grad(k_))&(fvc::grad(ome  ga_)));
   
    forAll(mesh_.cells(),celli)
    {
        if (kiKappa[celli] > scalar(0.0)) {
            kiKappa[celli] =  ((scalar(680.0)+sqr(kiKappa[celli]))/(scalar(400.0)+sqr(kiKappa[celli])));
          }
        else {
            kiKappa[celli] = scalar(1.0);
        }
    }
    return kiKappa;
}

However, I have not worked extensively with the tmp<>, hence not utterly sure it will work.

/ Niels

vkrastev August 31, 2010 11:36

Quote:

Originally Posted by ngj (Post 273501)
Perhaps you could do something like this:

Code:

tmp<volScalarField> kOmegaSI::fBetaStar() const
{
    volScalarField kiKappa=(1/(pow(omega_,3.0)))*((fvc::grad(k_))&(fvc::grad(ome  ga_)));
   
    forAll(mesh_.cells(),celli)
    {
        if (kiKappa[celli] > scalar(0.0)) {
            kiKappa[celli] =  ((scalar(680.0)+sqr(kiKappa[celli]))/(scalar(400.0)+sqr(kiKappa[celli])));
          }
        else {
            kiKappa[celli] = scalar(1.0);
        }
    }
    return kiKappa;
}

However, I have not worked extensively with the tmp<>, hence not utterly sure it will work.

/ Niels

No error message from the compiler! Now I can test the model in order to see If it actually works in conformity to my expectations. Thank you very much!

Regards

Vesselin

vkrastev August 31, 2010 12:51

Bad news...the source code is compiled without any error report, but when I try to start a test case (flow around a 2D airfoil at a Reynolds number of about 10^5-10^6) using the simpleFoam incompressible solver, I get the following error message (which for me is quite obscure...):

Selecting incompressible transport model Newtonian
Selecting RAS turbulence model kOmegaSI
kOmegaSICoeffs
{
betaZeroStar 0.09;
betaZero 0.072;
alpha 0.52;
alphaK 0.5;
alphaOmega 0.5;
}


Starting time loop

Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.0802683, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.077118, No Iterations 1
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.00991627, No Iterations 302
DICPCG: Solving for p, Initial residual = 0.0759295, Final residual = 0.000746425, No Iterations 25
DICPCG: Solving for p, Initial residual = 0.00187137, Final residual = 1.84519e-05, No Iterations 128
DICPCG: Solving for p, Initial residual = 0.000254188, Final residual = 2.47855e-06, No Iterations 151
time step continuity errors : sum local = 1.13767e-05, global = 2.50864e-07, cumulative = 2.50864e-07
#0 Foam::error::printStack(Foam::Ostream&) in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigSegv::sigSegvHandler(int) in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 memcpy in "/lib/tls/i686/cmov/libc.so.6"
#4 Uninterpreted:
#5 std::basic_string<char, std::char_traits<char>, std::allocator<char> >::basic_string(std::string const&) in "/home/vesselin/OpenFOAM/ThirdParty-1.6/gcc-4.3.3/platforms/linux/lib/libstdc++.so.6"
#6 Foam::IOobject::IOobject(Foam::word const&, Foam::fileName const&, Foam::objectRegistry const&, Foam::IOobject::readOption, Foam::IOobject::writeOption, bool) in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#7 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator*<double, Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/applications/bin/linuxGccDPOpt/simpleFoam"
#8 Foam::incompressible::RASModels::kOmegaSI::betaSta r() const in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#9 Foam::incompressible::RASModels::kOmegaSI::kiOmega () const in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#10 Foam::incompressible::RASModels::kOmegaSI::fBeta() const in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#11 Foam::incompressible::RASModels::kOmegaSI::correct () in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#12 main in "/home/vesselin/OpenFOAM/OpenFOAM-1.6/applications/bin/linuxGccDPOpt/simpleFoam"
#13 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#14 _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/i386/elf/start.S:122
Segmentation fault

Here It is the portion of the code containing all the functions correlated with fBetaStar:

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

tmp<volScalarField> kOmegaSI::fBetaStar() const
{
volScalarField kiKappa=(1/(pow(omega_,3.0)))*((fvc::grad(k_))&(fvc::grad(ome ga_)));
forAll(mesh_.cells(),celli)
{
if (kiKappa[celli] > scalar(0.0)) {
kiKappa[celli]=((scalar(680.0)+sqr(kiKappa[celli]))/(scalar(400.0)+sqr(kiKappa[celli])));
}
else {
kiKappa[celli]=scalar(1.0);
}
}
return kiKappa;
}

tmp<volScalarField> kOmegaSI::betaStar() const
{
return (Cmu_*fBetaStar());
}

tmp<volScalarField> kOmegaSI::kiOmega() const
{
return (mag((1.0/pow((betaStar()*omega_),3.0))*(((skew(fvc::grad(U_ )))& (skew(fvc::grad(U_))))&&(symm(fvc::grad(U_))))));
}

tmp<volScalarField> kOmegaSI::fBeta() const
{
return ((scalar(1.0)+70.0*kiOmega())/(scalar(1.0)+80.0*kiOmega()));
}

I'll go on trying to understand what all of this actually means....

V.

salehi144 March 6, 2012 08:51

Thanks
 
I just wanted to thank every body here.
I had the same problem. Then i solved it by reading this thread.
Thanks.

u396852 November 8, 2012 07:36

error "mesh_ was not declared in this scope"
 
I am trying to loop over the full mesh. I tried the posted suggestions. I am using OpenFOAM 1.5 version and 'myinterFoam' solver. I am trying to modify my viscosity through library files. I want to loop over the full mesh to impliment a yield stress condition. when I try to use 'forALL(mesh_.cells(), celli)' it gives compilation errors saying "mesh_ was not declared in this scope". could anybody comment?


All times are GMT -4. The time now is 02:50.