[351Cleveland]

Hello again forum,

I've been having trouble with a power function I've implemented in a modified compressible k-epsilon turbulence model.
The equation is as follows:

rEQ = Cr * Foam:ow((den_*yliq)/denliq_,0.133)*((st_/Foam:ow(epslion_,0.4)*63.1)

There are no errors during the compiling process and this identical equation has worked in a code based on the incompressible k-epsilon RAS model.

Cr is a constant and den, yliq , denliq and st are all volScalarFields; when attempting to solve, the following error is thrown up:

Starting time loop

Time = 1e-05

smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0637079, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 0.0881759, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 1, Final residual = 0.0805518, No Iterations 4
DILUPBiCG:  Solving for e, Initial residual = 0.985986, Final residual = 0.0452748, No Iterations 1
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0474283, No Iterations 16
time step continuity errors : sum local = 1.14264e-07, global = -2.41546e-08, cumulative = -2.41546e-08
rho max/min : 1.19018 1.1863
smoothSolver:  Solving for epsilon, Initial residual = 0.0494738, Final residual = 0.00120473, No Iterations 4
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 0.0654523, No Iterations 4
DILUPBiCG:  Solving for yliq, Initial residual = 1, Final residual = 0.0810717, No Iterations 171
DILUPBiCG:  Solving for bigSig, Initial residual = 1, Final residual = 0.0944473, No Iterations 253
DILUPBiCG:  Solving for yvap, Initial residual = 0, Final residual = 0, No Iterations 0
#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam220/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::sigFpe::sigHandler(int) in "/opt/openfoam220/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3   in "/lib/x86_64-linux-gnu/libm.so.6"
#4  Foam::pow(Foam::Field<double>&, Foam::UList<double> const&, double const&) in "/opt/openfoam220/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#5  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::pow<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::dimensioned<double> const&) in "/opt/openfoam220/platforms/linux64GccDPOpt/lib/libcompressibleRASModels.so"
#6  Foam::compressible::RASModels::yliqcomp::correct() in "/home/moffat/OpenFOAM/moffat-2.2.0/platforms/linux64GccDPOpt/lib/libuserRAScomp.so"
#7  
 in "/opt/openfoam220/platforms/linux64GccDPOpt/bin/rhoSimpleFoam"
#8  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#9  
 in "/opt/openfoam220/platforms/linux64GccDPOpt/bin/rhoSimpleFoam"
Floating point exception (core dumped)

I'm fairly sure this has to do with the pow term containing (den*yliq)/denliq
as when removing the yliq term, the equation solves normally with no errors, however yliq does not have any values at zero or negative as a filter function has been included to prevent this?

Can anybody shed some light on the problem?

Thanks.
Dom.

[Bernhard]

And what about denliq, as you are dividing by that?

Also, I would expect your solver to complain about epslion, instead of epsilon.

[351Cleveland]

denliq is set at a constant dimensioned value of 1000 - but is implemented as a field rather than a constant due to a number of other tests I was running, unfortunately substituting it with an integer or removing it altogether still results in the same problem.

I attribute epslion to my inconsistent typing It is epsilon in the equation itself.

[Bernhard]

Can you show you complete correct function? Especially your yvap equation, as it is not solving anything, according to your logfile.

[351Cleveland]

Certainly, I've included the yliq equation - yvap currently isn't fully implemented so won't be solving anything at this stage.

 rEQ_ = Cr_*(Foam::pow((den_*yliq_)/denliq_,0.1333))*(st_/(Foam::pow(epsilon_,0.4)*63.1));
  
 tmp<fvScalarMatrix> yliqEqn
      (
       fvm::ddt(rho_,yliq_)
       +fvm::div(phi_,yliq_)
       -fvm::Sp(fvc::div(phi_), yliq_)
       ==
       fvm::laplacian(DyliqEff(), yliq_)
      )

The transport equation was copied from the kinetic energy equation with the field changed - this equation solves with no issues which is what is puzzling.

[Bernhard]

How do you initialize epsilon?

[351Cleveland]

Fixed value (very turbulent flow) at twin air inlets and the compressible epsilon wall function at the boundaries with a value of 1e-12 and zeroGradient at the outlet - its my first compressible case so I'm adapting as best I can from the incompressible version.

[Bernhard]

I mean the internalField. If that is zero, it would explain the division by zero you are probably encountering.

[alexeym]

What is the sign of

(den_*yliq_)/denliq

?

As you receive error from libm I guess you're out of the function domain.

[351Cleveland]

Ah, my apologies for missing that out. That is also set at a fixed value of 0.3 - identical to one of the air inlet values. I've managed to get a modified form of the equation working as follows:

rEQ_ = Cr_*(Foam::pow((den_)/denliq_,0.1333))*(st_/(Foam::pow(epsilon_,0.4)*63.1));

However any attempt to add the yliq_ field to this results in the previous error - even when forcing all of the its values to be non-zero.

[351Cleveland]

Originally Posted by alexeym

What is the sign of

(den_*yliq_)/denliq

?

As you receive error from libm I guess you're out of the function domain.

From the working version of the equation in the above post, all values (without yliq included) return as positive, and all values for yliq calculated separately are filtered to either return zero or positive - currently all at a small value to try and get the equation working. However I'm not sure if this somehow produces an fpe error as it crashes on the first step before any results are written.

[Bernhard]

It is really difficult to tell what is going wrong, as you provide incomplete information. Clearly from the piece of code you are sharing, the problem arises after the solving (and construction) of yliqEqn. Please, share your whole code, or just try to old school debug with Info statements.

[351Cleveland]

After taking your advice, I built a debug version of openfoam 2.2.1 and ran my code using gdb as the debugger and found something very interesting from the error log:

                     P { margin-bottom: 0.21cm; }   Starting time loop  
 

 Time = 1e-05  
 

 smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0637079, No Iterations 4  
 smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 0.0881759, No Iterations 4  
 smoothSolver:  Solving for Uz, Initial residual = 1, Final residual = 0.0805518, No Iterations 4  
 DILUPBiCG:  Solving for e, Initial residual = 0.985937, Final residual = 0.045264, No Iterations 1  
 GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0474283, No Iterations 16  
 time step continuity errors : sum local = 1.14264e-07, global = -2.41546e-08, cumulative = -2.41546e-08  
 rho max/min : 1.19018 1.1863  
 smoothSolver:  Solving for epsilon, Initial residual = 0.0494738, Final residual = 0.00120473, No Iterations 4  
 smoothSolver:  Solving for k, Initial residual = 1, Final residual = 0.0654523, No Iterations 4  
 DILUPBiCG:  Solving for yliq, Initial residual = 1, Final residual = 0.0673715, No Iterations 254  
 DILUPBiCG:  Solving for bigSig, Initial residual = 1, Final residual = 0.094838, No Iterations 253  
 DILUPBiCG:  Solving for yvap, Initial residual = 0, Final residual = 0, No Iterations 0  
 

 Program received signal SIGFPE, Arithmetic exception.  
 0x00007ffff314f9e0 in ?? () from /lib/x86_64-linux-gnu/libm.so.6  
 (gdb) bt  
 #0  0x00007ffff314f9e0 in ?? () from /lib/x86_64-linux-gnu/libm.so.6  
 #1  0x00007ffff7652912 in Foam::pow (s=-4.8587039789792688e-10, e=0.1333)  
     at /home/moffat/FOAMDebug/OpenFOAM-2.2.1/src/OpenFOAM/lnInclude/doubleFloat.H:78 
 #2  0x00007ffff3f639c9 in Foam::pow (res=..., f1=..., s2=@0x7fffffff8cb0: 0.1333) at fields/Fields/scalarField/scalarField.C:118  
 #3  0x00007ffff6c4d434 in Foam::pow<Foam::fvPatchField> (f=..., f1=..., s=@0x7fffffff8cb0: 0.1333)  
     at /home/moffat/FOAMDebug/OpenFOAM-2.2.1/src/OpenFOAM/lnInclude/scalarFieldField.C:94 
 #4  0x00007ffff6c4b8c8 in Foam::pow<Foam::fvPatchField, Foam::volMesh> (tPow=..., gsf=..., ds=...)  
     at /home/moffat/FOAMDebug/OpenFOAM-2.2.1/src/OpenFOAM/lnInclude/GeometricScalarField.C:274 
 #5  0x00007ffff6c66ffe in Foam::pow<Foam::fvPatchField, Foam::volMesh> (tgsf=..., ds=...) 
     at /home/moffat/FOAMDebug/OpenFOAM-2.2.1/src/OpenFOAM/lnInclude/GeometricScalarField.C:326 
 #6  0x00007ffff6c6635b in Foam::pow<Foam::fvPatchField, Foam::volMesh> (tgsf=..., s=@0x7fffffffa098: 0.1333)  
     at /home/moffat/FOAMDebug/OpenFOAM-2.2.1/src/OpenFOAM/lnInclude/GeometricScalarField.C:350 
 #7  0x00007fffee839135 in Foam::compressible::RASModels::yliqcomp::correct (this=0x9f97e30) at yliqcomp/yliqcomp.C:826  
 #8  0x000000000044d252 in main (argc=1, argv=0x7fffffffd1e8) at rhoSimpleFoam.C:68  
 (gdb)

Although I'm not too familiar with the s and e synatx relating to the pow function in #1, I'm guessing that s is the sum of the terms and e is the exponent?
Since the s term term comes out negative despite the filter function attached to yliq - can I assume that it is not working correctly and producing -ve values for yliq? Or does the problem lie elsewhere in the term?

Thanks for the assistance.

[alexeym]

I'm not quite sure you've shown right piece of code. According to solver output yliq equation is solved. Can you show the code after solution of yvap equation?

[351Cleveland]

I've included the yliq equation and the correction routine along with the equations implemented immediately after in the .C file. Although the yliq field seems to solves without issue, it is only when being used to calculate rEQ_ that the problem seems to occur.

I'm beginning to wonder if the correction routine I put in (based on the number of cells in the mesh) only takes effect after everything is solved - as the output file for yliq contains no negative values at all.....but seems strange that it worked in exactly the same form for the incompressible version of the same model I made.

                     P { margin-bottom: 0.21cm; }    tmp<fvScalarMatrix> yliqEqn
       (
        fvm::ddt(rho_,yliq_)
        +fvm::div(phi_,yliq_)
        -fvm::Sp(fvc::div(phi_), yliq_)
        ==
        fvm::laplacian(DyliqEff(), yliq_)
        );
 

     yliqEqn().relax();
     solve(yliqEqn);
 

     //correction scheme
 

        int n;
     for (n=0; n<=1564000;)
       {
         if(yliq_[n]<0){
           yliq_[n]=0;
         }
         n++;
     }
 

 

     mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 

     // Re-calculate thermal diffusivity
     alphat_ = mut_/Prt_;
     alphat_.correctBoundaryConditions();
 

         tauC_ = 4.35*(k_/epsilon_);
     tauC_.correctBoundaryConditions();
 

        den_ = 1/(((1-yliq_-yvap_)/denair_)+(yliq_/denliq_)+(yvap_/rhovap_));
     den_.correctBoundaryConditions();
 

     //ph2_ = pow(yliq_,0.5);
 

     rEQ_ = Cr_*(Foam::pow((den_*yliq_)/denliq_,0.1333))*(st_/(Foam::pow(epsilon_,0.4)*63.1));
     rEQ_.correctBoundaryConditions();

Cheers.

[alexeym]

Why not use either

forAll(yliq_, cellI)
{
    if(yliq_[cellI] < 0)
        yliq_[cellI] = 0;
}

if you like cycles

or

yliqMin = dimensionedScalar("yliqMin", ...dimensions..., 0.0);
yliq_ = max(yliq_, yliqMin);

or

yliqMin = dimensionedScalar("yliqMin", ...dimensions..., 0.0);
yliq_ = bound(yliq_, yliqMin);

?

[351Cleveland]

It appears that's solved the problem nicely - looks like the routine I used wasn't working correctly and causing negative values to slip through into the following equations.

I would like to thank you all very much for your help, hopefully this is the last time the problem will show up.

Thanks again.
Dominic.

[Bernhard]

Dear Dominic, I want to emphasize that the problem was solved within 20 minutes after posting the relevant lines. Please keep in mind to carefully prepare you post when asking a question here. If you reread the thread, you will see that I already asked you to post the relevant piece of the code two days ago, which could have saved you a lot of time.