
[Sponsors] 
October 2, 2013, 07:42 
Foam::pow function error

#1 
New Member
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5 
Hello again forum,
I've been having trouble with a power function I've implemented in a modified compressible kepsilon 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 kepsilon 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 = 1e05 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.14264e07, global = 2.41546e08, cumulative = 2.41546e08 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_64linuxgnu/libc.so.6" #3 in "/lib/x86_64linuxgnu/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/moffat2.2.0/platforms/linux64GccDPOpt/lib/libuserRAScomp.so" #7 in "/opt/openfoam220/platforms/linux64GccDPOpt/bin/rhoSimpleFoam" #8 __libc_start_main in "/lib/x86_64linuxgnu/libc.so.6" #9 in "/opt/openfoam220/platforms/linux64GccDPOpt/bin/rhoSimpleFoam" Floating point exception (core dumped) 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: 12 
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: 5 
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: 12 
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: 5 
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_) ) 

October 2, 2013, 09:37 

#6 
Senior Member
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12 
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: 5 
Fixed value (very turbulent flow) at twin air inlets and the compressible epsilon wall function at the boundaries with a value of 1e12 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: 12 
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: 715
Rep Power: 12 
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: 5 
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)); 

October 2, 2013, 10:01 

#11 
New Member
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5 
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: 12 
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: 5 
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 { marginbottom: 0.21cm; } Starting time loop Time = 1e05 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.14264e07, global = 2.41546e08, cumulative = 2.41546e08 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_64linuxgnu/libm.so.6 (gdb) bt #0 0x00007ffff314f9e0 in ?? () from /lib/x86_64linuxgnu/libm.so.6 #1 0x00007ffff7652912 in Foam::pow (s=4.8587039789792688e10, e=0.1333) at /home/moffat/FOAMDebug/OpenFOAM2.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/OpenFOAM2.2.1/src/OpenFOAM/lnInclude/scalarFieldField.C:94 #4 0x00007ffff6c4b8c8 in Foam::pow<Foam::fvPatchField, Foam::volMesh> (tPow=..., gsf=..., ds=...) at /home/moffat/FOAMDebug/OpenFOAM2.2.1/src/OpenFOAM/lnInclude/GeometricScalarField.C:274 #5 0x00007ffff6c66ffe in Foam::pow<Foam::fvPatchField, Foam::volMesh> (tgsf=..., ds=...) at /home/moffat/FOAMDebug/OpenFOAM2.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/OpenFOAM2.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) 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: 715
Rep Power: 12 
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: 5 
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 { marginbottom: 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(); // Recalculate thermal diffusivity alphat_ = mut_/Prt_; alphat_.correctBoundaryConditions(); tauC_ = 4.35*(k_/epsilon_); tauC_.correctBoundaryConditions(); den_ = 1/(((1yliq_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(); 

October 4, 2013, 10:44 

#16 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 715
Rep Power: 12 
Why not use either
Code:
forAll(yliq_, cellI) { if(yliq_[cellI] < 0) yliq_[cellI] = 0; } or Code:
yliqMin = dimensionedScalar("yliqMin", ...dimensions..., 0.0); yliq_ = max(yliq_, yliqMin); Code:
yliqMin = dimensionedScalar("yliqMin", ...dimensions..., 0.0); yliq_ = bound(yliq_, yliqMin); 

October 4, 2013, 11:31 

#17 
New Member
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 5 
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: 12 
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.


Thread Tools  
Display Modes  


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 OpenFOAM1.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 