[A_Pete]

Hi everyone!

I have implemented the correction of turbulent viscosity by Reboud in interface regions for cavitating flows. Physically it does everything I want pretty well, but from time to time my solutions will abort due to floating point exception errors.
For me this is an indication that I didn't prevent the solver from dividing by zero somewhere. So I found a spot where this might happen during the solution, altough it shouldn't and implemented a small prevention.
But that still doesn't get the job done, because I will still get the same FPEs and have no idea where to look for. This is often accompanied by omega getting very large first. But since this doesn't happen with other turbulence models in my simulations something has to be wrong with my code.

What I have done is take the standard kOmegaSST turbulence model of OpenFOAM and modify it the way I need it.

The error says this:

[0] #0  Foam::error::printStack(Foam::Ostream&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #1  Foam::sigFpe::sigHandler(int) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #2   in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #3  double Foam::sumProd<double>(Foam::UList<double> const&,
Foam::UList<double> const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #4  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&,
unsigned char) const in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #5  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
[0] #6  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #7  Foam::fvMatrix<double>::solve() in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #8  Foam::SolverPerformance<double>
Foam::solve<double>(Foam::tmp<Foam::fvMatrix<double> > const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #9  Foam::incompressible::RASModels::modKOmegaSST::correct() in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libincompressibleRASModels.so"
[0] #10  
[0]  in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #12  
[0]  in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[prop-04:28417] *** Process received signal ***
[prop-04:28417] Signal: Floating point exception (8)
...

Something seems to go wrong at #9 in the correct() method of modKOmegaSST, but I can't find it. Below are what I suppose to be the interesting parts of my code regarding this issue.


In the header file I added this:

// under member functions
//- Calculate Reboud's correction and return the normalized density function

virtual volScalarField normDensFunc();

And in the source file the interesting parts are the following:

// definition of the calculated function under member functions

volScalarField modKOmegaSST::normDensFunc()
{
    //- Reading fields for Reboud's correction
    const volScalarField& alpha1 = mesh_.lookupObject<volScalarField>("alpha1");
    const volScalarField& alpha2 = mesh_.lookupObject<volScalarField>("alpha2");

    //- Lookup transportProperties dictionary to recieve fluid properties
    const dictionary& transportProperties =
        mesh_.lookupObject<IOdictionary>("transportProperties");
    const dictionary& phase1 = transportProperties.subDict("phase1");
    const dictionary& phase2 = transportProperties.subDict("phase2");
    const dimensionedScalar& rho1 = phase1.lookup("rho");
    const dimensionedScalar& rho2 = phase2.lookup("rho");

    //- Calculation of mixture density
    const volScalarField rho_m = alpha1*rho1 + alpha2 * rho2;

    //- Calculation of density function
    const volScalarField f_rho =
        pow(( rho_m - rho2 ), nReb_) / pow((rho1-rho2), (nReb_-1)) + rho2;
    
    //- Calculation of ratio of density function and mixture density
    const volScalarField fRhoNorm = 
        f_rho / max(rho_m, rho1);

    //const volScalarField& test = mesh_.lookupObject<volScalarField>("alpha1");
    return fRhoNorm; 
}

...

void modKOmegaSST::correct()
{
...

// function is called and used for calculation of nut_ below

    //- Get normalized density function for Reboud's correction
    const volScalarField& fRhoNorm_(normDensFunc());

    // Re-calculate viscosity; standard calculation is the same, just w/o fRhoNorm
    nut_ = a1_*k_*fRhoNorm_/max(a1_*omega_, b1_*F23()*sqrt(S2));
    nut_.correctBoundaryConditions();

}

So if anybody has any idea what I could try out I would be very happy and thankful.


Thanks in advance.

Best regards,
Andreas

[A_Pete]

Any hints?

[GPesch]

I dont actually have a suggestion of how to avoid the FPE but I just want to give you a hint on how to find at which part of the code there is something going wrong:

a. Either get a debug version of OpenFOAM in which you can directly find out which line is causing the error.

b. Use Info<<"..."<<endl; to write some lines into the terminal at relevant parts of the code. From that you can see up to which point your code is still running. From that you can narrow done the line of code which produces the error..
That was helping me always a lot (since I didnt had a debug version of OF)

[Kappa]

Hey A_Pete,

I got the same problem. Can you tell me how you fixed it?

Thanks