# Foam::pow function error

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

 October 2, 2013, 07:42 Foam::pow function error #1 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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&, Foam::UList const&, double const&) in "/opt/openfoam220/platforms/linux64GccDPOpt/lib/libOpenFOAM.so" #5 Foam::tmp > Foam::pow(Foam::tmp > const&, Foam::dimensioned 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.

 October 2, 2013, 08:45 #2 Senior Member   Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 20 And what about denliq, as you are dividing by that? Also, I would expect your solver to complain about epslion, instead of epsilon.

 October 2, 2013, 09:02 #3 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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.

 October 2, 2013, 09:19 #4 Senior Member   Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 20 Can you show you complete correct function? Especially your yvap equation, as it is not solving anything, according to your logfile.

 October 2, 2013, 09:30 #5 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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 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.

 October 2, 2013, 09:37 #6 Senior Member   Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 20 How do you initialize epsilon?

 October 2, 2013, 09:41 #7 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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.

 October 2, 2013, 09:47 #8 Senior Member   Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 20 I mean the internalField. If that is zero, it would explain the division by zero you are probably encountering.

 October 2, 2013, 09:53 #9 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 37 What is the sign of Code: `(den_*yliq_)/denliq` ? As you receive error from libm I guess you're out of the function domain.

 October 2, 2013, 09:54 #10 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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.

October 2, 2013, 10:01
#11
New Member

Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 14
Quote:
 Originally Posted by alexeym 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.

 October 2, 2013, 10:07 #12 Senior Member   Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 20 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.

 October 4, 2013, 06:53 #13 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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 (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 (tPow=..., gsf=..., ds=...) at /home/moffat/FOAMDebug/OpenFOAM-2.2.1/src/OpenFOAM/lnInclude/GeometricScalarField.C:274 #5 0x00007ffff6c66ffe in Foam::pow (tgsf=..., ds=...) at /home/moffat/FOAMDebug/OpenFOAM-2.2.1/src/OpenFOAM/lnInclude/GeometricScalarField.C:326 #6 0x00007ffff6c6635b in Foam::pow (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.

 October 4, 2013, 08:56 #14 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 37 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?

 October 4, 2013, 10:23 #15 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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 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.

 October 4, 2013, 10:44 #16 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 37 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.

 October 4, 2013, 11:31 #17 New Member   Dominic Join Date: Jan 2011 Location: Leeds, UK Posts: 25 Rep Power: 14 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.

 October 4, 2013, 11:35 #18 Senior Member   Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 20 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.