CFD Online URL
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Foam::pow function error

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree3Likes
  • 2 Post By alexeym
  • 1 Post By Bernhard

Reply
 
LinkBack Thread Tools Display Modes
Old   October 2, 2013, 07:42
Default Foam::pow function error
  #1
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
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

Code:
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.
351Cleveland is offline   Reply With Quote

Old   October 2, 2013, 08:45
Default
  #2
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
And what about denliq, as you are dividing by that?

Also, I would expect your solver to complain about epslion, instead of epsilon.
Bernhard is offline   Reply With Quote

Old   October 2, 2013, 09:02
Default
  #3
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
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.
351Cleveland is offline   Reply With Quote

Old   October 2, 2013, 09:19
Default
  #4
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
Can you show you complete correct function? Especially your yvap equation, as it is not solving anything, according to your logfile.
Bernhard is offline   Reply With Quote

Old   October 2, 2013, 09:30
Default
  #5
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
Certainly, I've included the yliq equation - yvap currently isn't fully implemented so won't be solving anything at this stage.

Code:
 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.
351Cleveland is offline   Reply With Quote

Old   October 2, 2013, 09:37
Default
  #6
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
How do you initialize epsilon?
Bernhard is offline   Reply With Quote

Old   October 2, 2013, 09:41
Default
  #7
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
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.
351Cleveland is offline   Reply With Quote

Old   October 2, 2013, 09:47
Default
  #8
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
I mean the internalField. If that is zero, it would explain the division by zero you are probably encountering.
Bernhard is offline   Reply With Quote

Old   October 2, 2013, 09:53
Default
  #9
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 715
Rep Power: 12
alexeym will become famous soon enough
What is the sign of

Code:
(den_*yliq_)/denliq
?

As you receive error from libm I guess you're out of the function domain.
alexeym is offline   Reply With Quote

Old   October 2, 2013, 09:54
Default
  #10
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
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:

Code:
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 is offline   Reply With Quote

Old   October 2, 2013, 10:01
Default
  #11
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
Quote:
Originally Posted by alexeym View Post
What is the sign of

Code:
(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.
351Cleveland is offline   Reply With Quote

Old   October 2, 2013, 10:07
Default
  #12
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
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.
Bernhard is offline   Reply With Quote

Old   October 4, 2013, 06:53
Default
  #13
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
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:

Code:
                     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.
351Cleveland is offline   Reply With Quote

Old   October 4, 2013, 08:56
Default
  #14
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 715
Rep Power: 12
alexeym will become famous soon enough
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?
alexeym is offline   Reply With Quote

Old   October 4, 2013, 10:23
Default
  #15
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
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.

Code:
                     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.
351Cleveland is offline   Reply With Quote

Old   October 4, 2013, 10:44
Default
  #16
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 715
Rep Power: 12
alexeym will become famous soon enough
Why not use either

Code:
forAll(yliq_, cellI)
{
    if(yliq_[cellI] < 0)
        yliq_[cellI] = 0;
}
if you like cycles

or

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

Code:
yliqMin = dimensionedScalar("yliqMin", ...dimensions..., 0.0);
yliq_ = bound(yliq_, yliqMin);
?
351Cleveland and Thamali like this.
alexeym is offline   Reply With Quote

Old   October 4, 2013, 11:31
Default
  #17
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5
351Cleveland is on a distinguished road
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.
351Cleveland is offline   Reply With Quote

Old   October 4, 2013, 11:35
Default
  #18
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
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.
Thamali likes this.
Bernhard is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Pressure outlet boundary condition rolando OpenFOAM Running, Solving & CFD 58 November 12, 2014 14:51
Ansys Fluent 13.0 UDF compilation problem in Window XP (32 bit) Yogini Fluent UDF and Scheme Programming 7 October 3, 2012 08:24
Building OpenFOAM1.7.0 from source ata OpenFOAM Installation 42 May 14, 2012 21:48
latest OpenFOAM-1.6.x from git failed to compile phsieh2005 OpenFOAM Bugs 25 February 9, 2010 05:37
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 18:51


All times are GMT -4. The time now is 19:35.