[FaroukH]

Hello Foamers,

I am relatively new to OpenFOAM programming and C++ but have already had some experience with simple solver modifications and implementing coded function objects. It could be that I am stating the wrong issue here since I am still new. Thank you very much in advance, and I appreciate any hints that set me in the correct direction.

I have been trying to implement the Sigma (Singular Values based) sub-grid scale (SGS) LES model from Nicoud et al. by modifying the WALE model and referencing the Smagorinsky one that are already implemented in OpenFOAM. The work is largely based on someone else who did this exact same task in OpenFOAM 2.1 that does not compile with new versions of OpenFOAM. What I have done so far is try and figure out how the implementation would work in OpenFOAM 7 where a lot of syntax has changed but still get a lot of errors.

I am not sure where my errors are coming from, and I have tried learning more C++ and about OpenFOAM programming. The main error I am getting for all my declared functions is that the prototype does not match any class and the redefinition gives an error. This leads me to believe it’s coming from my incorrect implementations of the BasicTurbulenceModel template. An example of such an error can be seen below.

In file included from SigmaSV.H:161:0,
                 from SigmaSV.C:22:
SigmaSV.C:35:21: error: prototype for ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::LESModels::SigmaSV<BasicTurbulenceModel>::k(const Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> >&) const’ does not match any in class ‘Foam::LESModels::SigmaSV<BasicTurbulenceModel>’
 tmp<volScalarField> SigmaSV<BasicTurbulenceModel>::k
                     ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from SigmaSV.C:22:0:
SigmaSV.H:134:37: error: candidates are: virtual Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::LESModels::SigmaSV<BasicTurbulenceModel>::k() const
         virtual tmp<volScalarField> k() const

The section that I believe is causing trouble in the .C file. It is the section dedicated for protected member functions:

template<class BasicTurbulenceModel>
tmp<volScalarField> SigmaSV<BasicTurbulenceModel>::k
(
    const tmp<volTensorField>& gradU
) const
{
    volSymmTensorField D(symm(gradU));

    volScalarField a(this->Ce_/this->delta());
    volScalarField b((2.0/3.0)*tr(D));
    volScalarField c(2*Ck_*this->delta()*(dev(D) && D));

    return volScalarField::New
    (
        IOobject::groupName("k", this->alphaRhoPhi_.group()),
        sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a))
    );
}

template<class BasicTurbulenceModel>
tmp<volScalarField> SigmaSV<BasicTurbulenceModel>::DSigma
(
    const tmp<volTensorField>& gradU
)
{
    dimensionedScalar smallev("smallev", dimensionSet(0,0,-2,0,0,0,0), SMALL);
    dimensionedScalar dimzero("dimzero", dimensionSet(0,0,-2,0,0,0,0), 0.0);

    volTensorField G(dev(T(gradU)) & dev(gradU)); //implemented correct implementation with dev2 and no dev at transposed (T() operator) gradU since the bug mentioned by schneider has been fixed

    volVectorField eigenvalues
    (
        IOobject
        (
            "eigenvalues",
            "eigenvalues",
            this->mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        this->mesh_,
        dimensionedVector("eigenvalues", dimensionSet(0,0,-2,0,0,0,0), vector::zero)
    );

    forAll(G, cellI)
    {
        eigenvalues[cellI] = calcEigenvalues(G[cellI]);
    }

    volScalarField sigma1(sqrt(max(eigenvalues.component(0),dimzero)));
    volScalarField sigma2(sqrt(max(eigenvalues.component(1),dimzero)));
    volScalarField sigma3(sqrt(max(eigenvalues.component(2),dimzero)));
    // max(.,0) is necessary because the eigenvalues may get slightly negative due to roundoff errors
    volScalarField dSigmaCalc(sigma3*(sigma1-sigma2)*(sigma2-sigma3) / (sqr(sigma1)+smallev));
    return  dSigmaCalc;// a mistake was detected, wherer sigma1 in the denominator has to be quared, however the square root was taken sqr
}


template<class BasicTurbulenceModel>
void SigmaSV<BasicTurbulenceModel>::correctNut()
{

    this->nut_ = sqr(Csigma*this->delta()) * this->DSigma(fvc::grad(this->U_));// no need for multiplying by rho since we are taken nut and this is done later automatically for compressible flows. Also a mistake with sqrt not a square for csigma*delta.
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

Furthermore, I am using the smagorinsky implementation for calculating the subgrid scale kinetic energy and diffusivity.

My Make, C, and H files can be found here (SigmaSV), as well as the original implementations (SingularValOrg):
https://www.dropbox.com/sh/7394svgvi...5z57IHkQa?dl=0

Thank you

[ztnuaa]

Dear Farouk,
I'm more familiar with OF-2.3.1 and OF-5.x , but I think in OF-7 it may be similar as OF-5.x when we want to add a new SGS model.

After check your files, I think the error did occur in your Make "files".

Here is some advice to add a new SGS model in the original OF code. (My test is under OF-5.x):

1. Copy your SigmaSV model files to:

src/TurbulenceModels/turbulenceModels/LES/SigmaSV

2. Modify those two files:
(1) Compressible turbulence model

src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C

(2) Incompressible turbulence model

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C

Add those two lines to the end of them:

#include "SigmaSV.H"
makeLESModel(SigmaSV);

[FaroukH]

Dear Timo,

thank you very much for the advice! I have been doing something similar by creating two new libraries for the Sigma sgs model by following this setup: http://hassankassem.me/posts/newturbulencemodel6/ and adapting it to Sigma.

I am trying to apply your explanation but I am not sure what to do next. Should I do this in a copy of openFOAM source file locally in my home folder and then compiling tubulenceModels?

Thank You,
Farouk

[ztnuaa]

Dear Farouk,

The tutorial you followed is okay too, it will compile a lib named *libmyincompressibleTurbulenceModels*, which you can include in the controlDict file.
However, you did something wrong apparently. I recommend you to check the steps again.

Back to my method, yes, you should do it on the OpenFOAM source code in your home folder, say ~/OpenFOAM/OpenFOAM-7/src. Next:
1. clean turbulenceModels lib

cd $FOAM_SRC/TurbulenceModels/turbulenceModels
wclean

2. compile the src folder and application folder

cd $FOAM_SRC
./Allwmake
cd  $FOAM_APP
./Allwmake

By the way, I don't know the problem you are solving. If you are using a compressible solver , then the tutorial post you followed is not okay since it only compiles for the incompressible turbulence models. You should take care of the releated Make files/options .

Timo

[FaroukH]

Dear Timo,

I will attempt these steps that you suggested. I beleive it would be a better approach. Thank you!

Farouk

[ztnuaa]

Dear Farouk,
I tried to debug your original code. Here are some corrections may be helpful to you.
1. SigmaSV.H line 82:83, the following correction will fix the compiling error you post.

    tmp<volScalarField> k(const volTensorField& gradU) const;
    tmp<volScalarField> DSigma(const volTensorField& gradU);

should be

    tmp<volScalarField> k(const tmp<volTensorField>& gradU) const;
    tmp<volScalarField> DSigma(const tmp<volTensorField>& gradU);

2. SigmaSV.C
When return volScalarField try use:

    
return tmp<volScalarField>
    (
        new volScalarField
        (
        IOobject::groupName("k", this->alphaRhoPhi_.group()),
        sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a))
        )
    );

NOT:

    return volScalarField::New
    (
        IOobject::groupName("k", this->alphaRhoPhi_.group()),
        sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a))
    );

And at line 267, I quoted it here:

Info <<<BasicTurbulenceModel> "\nIn SingularValues::calcEigenvalues: acosarg is out of range, value is " << acosarg << endl;

I think you are trying to output the model type name? It's not correct to use it like this. especially you have "<<<" which will make compiler complains because "<<" is Output "<<<" is nothing. Just don't output BasicTurbulenceModel will fix it.


Following the tutorial you mentioned in the last post, I had my test files attached here. Again, it's tested under OF-5.x and compile passed.

Timo

[FaroukH]

Dear Timo,

thank you very much for helping in finding the issues. I also addressed some other ones in the past week after learning a bit more about how the code is to be structured. I will run a few simulations to validate.

The only thing with following this tutorial that is linked, should I create two libraries? One for the compressible case, and the other for the incompressible case?

Farouk

[ztnuaa]

Dear Farouk,

I'm not sure whether this is a good way or not. And I don't know how to do it.

But it will work if you create two libs for compressible and incompressible turbulence respectively.

Timo

[FaroukH]

Hello Timo,


The code is compiling now. However I am running into one warning and and error while running the solver. I decided to use the method of implementing a new library since I do not have root access on the machine that will run the simulation and it takes too with ./Allmake. Furthermore, I tested the method by implementing the WALE model by just changing the name to CWALE and it works and runs.

I re-implmented the Sigma model into the WALE.C and WALE.H files. This is done since i need to calculate k() and epsilon() for some post processing that I do with a rhoPimpleFoam solver modified to output epsilon and k_sgs (the case I am applying this to showed better results under WALE than Smagorinsky). However, I do not need those values to directly calculate nut() with the Sigma model.


The warning that I get when running blockMesh:

--> FOAM Warning : 
    From function void* Foam::dlOpen(const Foam::fileName&, bool)
    in file POSIX.C at line 1251
    dlopen error : /home/itsnas/uuein/OpenFOAM/uuein-7/platforms/linux64GccDPInt32Opt/lib/libmycompressibleTurbulenceModel.so: undefined symbol: _ZTIN4Foam27compressibleTurbulenceModelE
--> FOAM Warning : 
    From function bool Foam::dlLibraryTable::open(const Foam::fileName&, bool)
    in file db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C at line 105
    could not load "libmycompressibleTurbulenceModel.so"

• Note that this warning does not necessarily effect the code, since when running rhoPimpleFoam, openFOAM recognizing the turbulence model Sigma and its relevant parameter CSigma
• It also occurs when running blockMesh for the custom WALE model and has no issue

The error that I get when running rhoPimpleFoam:

new cannot satisfy memory request.
This does not necessarily mean you have run out of virtual memory.
It could be due to a stack violation caused by e.g. bad use of pointers or an out of date shared library
Abort (core dumped)
Exit 134

• This does not occur when running the custom WALE variant

To test this, it is possible to use the pitzDaily compressible LES tutorial from openFOAM and modifying the following:

• changing LES model to Sigma in costant folder
• adding the the following above functions{...} in controlDict

libs                  
(
    "libmycompressibleTurbulenceModel.so"    
);

I have attached the file with the code I am using. The make and WALE base code work together (I tested it), however the implementation of the Sigma functions and parameters are causing the issues stated above.

[ztnuaa]

Dear FaroukH,

I downloaded your code and compiled under OF-7.

I think the error occurs at line 119:120 of Sigma.C

    volScalarField dSigmaCalc(sigma3*(sigma1-sigma2)*(sigma2-sigma3) / (sqr(sigma1)+smallev));
    return  dSigmaCalc;

If you want to return dSigmaCalc, try this:

    tmp<volScalarField> dSigmaCalc(sigma3*(sigma1-sigma2)*(sigma2-sigma3) / (sqr(sigma1)+smallev));
    return  dSigmaCalc;

Sorry, I'm not good at C++ too. I don't know why we must use tmp<volScalarField>.

Timo

[FaroukH]

Dear Time,

it works! I had some issues with negative eigen vectors and this fixes it at line 91:

  volTensorField G(dev(gradU.T()) & dev(gradU));

I will run it with the simulation I am working on and see how it goes! Thank you very much

Farouk

[ztnuaa]

Dear Farouk,

You are welcome.

I've been working on SGS models several months ago and SIGMA model is one of them.

According to Nicoud's paper: g is the local velocity graident tensor and G equals g^t g. And why you use the deviatoric part of g^t and g to build G? Are these two formula equivalent?

In my implementation I met some errors/bugs during the call of eigenValues (Maybe there are some bugs in OF-2.3.1 eigenValues). So I use method B(In Nicoud's paper) to calculate eigen values.

Timo

[FaroukH]

Dear Timo,

yes there are some error in eigenValues in openFOAM 2.3.1. However, I beleive those have been since fixed. The code I am referencing used the deviatoric parts of velocity gradients to eliminate issues with eigenValues. I will have to look into this more over the next two days. I will also take your advice try out Nicoud's Solver B settings.

Generally I have been referencing an implementation in openFOAM 2.x from 2015. There, a special eigenValues function was also used for their calculations.

Thank You,
Farouk

[Wang Shang]

Dear Zhang,
It's been a long time, but I'd like to ask if you implemented the Sigma SGS model and tested well. If possible, whether I can get the source code, I am very interested in testing this model.

[ztnuaa]

Hi Wang,
I didn't use SIGMA model in my simulations. However, there is a community implementation by :Zhi CHEN
Check his code if you like. Since he's from Peking University, I would like to trust its performance.

[Wang Shang]

Very useful. Thank you!