[enoch]

I've in trouble with resolving the following problem.
I added a following lift force to liftDragCoeffs.H. There is no compilation error, but when I run a case, I've got Floating point exception. I found this problem takes place due to the lift force defined as follows. What's wrong with the below coding?

volVectorField shearL
(
Cl*scalar(1)/(beta*alpha + scalar(0.001))
*alpha*6.0*6.46/(4.0*Pi_kim)
*pow(rhob, 0.5)*pow(nub*rhob, 0.5)/da
*(symmGradUb & Ur)/pow(magsymmGradUb, 0.5)
);


==========
#0 Foam::error:rintStack(Foam::Ostream&) in "/home/OpenFOAM/openfoam210/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/home/OpenFOAM/openfoam210/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3
in "/home/OpenFOAM/openfoam210/platforms/linux64GccDPOpt/bin/mixtureBloodYeles"
#4
in "/home/OpenFOAM/openfoam210/platforms/linux64GccDPOpt/bin/mixtureBloodYeles"
#5 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#6
in "/home/OpenFOAM/openfoam210/platforms/linux64GccDPOpt/bin/mixtureBloodYeles"
Floating point exception

[kmooney]

A got way to approach these types of problems is to recompile the code in debug mode and re-run the case with a debugging program. In linux the process looks a little like this:

From your solver source code:

wclean
export WM_COMPILE_OPTION=Debug
wmake

Then from your case directory:

gdb nameOfYourSolver
r

Then once your get the floating point exception:

bt

That should return a lot more detail about where the floating point exception is occurring.

Cheers!
Kyle

[enoch]

Thanks for your tip. I found from gdb debugging that there might be a problem with "divided by zero".

Below is a screen from the gdb debugging.

Program received signal SIGFPE, Arithmetic exception.
0x000000000047d3a6 in Foam::divideOp3<double, double, double>:perator() (this=0x7fffffff47f0,
x=@0x7359f20, y=@0x7fffffff4800) at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/ops.H:133
133 Op(divide, x / y)
(gdb)

---
y could be locally zero, which is the volSclarField having unit dimensions.
If y is not-dimensional, I can easily deal it with
x/(y+0.0001)

But, when y(=pow(magsymmGradUb, 0.5)) has a dimension, how can I handle it to avoid "divided by zero"?


Thanks in advance for your help.

[kmooney]

Originally Posted by enoch

Thanks for your tip. I found from gdb debugging that there might be a problem with "divided by zero".

Below is a screen from the gdb debugging.

Program received signal SIGFPE, Arithmetic exception.
0x000000000047d3a6 in Foam::divideOp3<double, double, double>:perator() (this=0x7fffffff47f0,
x=@0x7359f20, y=@0x7fffffff4800) at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/ops.H:133
133 Op(divide, x / y)
(gdb)

---
y could be locally zero, which is the volSclarField having unit dimensions.
If y is not-dimensional, I can easily deal it with
x/(y+0.0001)

But, when y(=pow(magsymmGradUb, 0.5)) has a dimension, how can I handle it to avoid "divided by zero"?


Thanks in advance for your help.

I'm surprised that it didn't give you more information on the floating point exception. If your solver/library was compiled with the debug flag, it should give a full listing of the functions that were called at the time. Are you sure you entered 'bt' after gdb output that last line?

[enoch]

Oops. I forgot to type "bt". Below is the screen that I've got. Do I have a problem with "divided by zero"???

(gdb) bt
#0 0x000000000047d3a6 in Foam::divideOp3<double, double, double>:perator() (this=0x7fffffff47f0,
x=@0x7359f20, y=@0x7fffffff4800) at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/ops.H:133
#1 0x000000000046a6f9 in VectorSpaceOps<3, 0>:pVS<Foam::Vector<double>, Foam::VectorSpace<Foam::Vector<double>, double, 3>, double, Foam::divideOp3<double, double, double> > (vs=..., vs1=..., s=@0x7fffffff4800,
o=...) at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/VectorSpaceM.H:34
#2 0x0000000000454fbb in Foam:perator/<Foam::Vector<double>, double, 3> (vs=..., s=0)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/VectorSpaceI.H:574
#3 0x0000000000484d43 in Foam::divide<Foam::Vector<double> > (res=..., f1=..., f2=...)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/FieldFunctions.C:556
#4 0x0000000000474e3c in Foam::divide<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> (res=...,
gf1=..., gf2=...)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:564
#5 0x000000000045eec0 in Foam:perator/<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> (
tgf1=..., tgf2=...)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:564
#6 0x0000000000446ae9 in main (argc=1, argv=0x7fffffffd3b8) at liftDragCoeffs.H:41
(gdb)

[kmooney]

Originally Posted by enoch

Oops. I forgot to type "bt". Below is the screen that I've got. Do I have a problem with "divided by zero"???

(gdb) bt
#0 0x000000000047d3a6 in Foam::divideOp3<double, double, double>:perator() (this=0x7fffffff47f0,
x=@0x7359f20, y=@0x7fffffff4800) at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/ops.H:133
#1 0x000000000046a6f9 in VectorSpaceOps<3, 0>:pVS<Foam::Vector<double>, Foam::VectorSpace<Foam::Vector<double>, double, 3>, double, Foam::divideOp3<double, double, double> > (vs=..., vs1=..., s=@0x7fffffff4800,
o=...) at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/VectorSpaceM.H:34
#2 0x0000000000454fbb in Foam:perator/<Foam::Vector<double>, double, 3> (vs=..., s=0)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/VectorSpaceI.H:574
#3 0x0000000000484d43 in Foam::divide<Foam::Vector<double> > (res=..., f1=..., f2=...)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/FieldFunctions.C:556
#4 0x0000000000474e3c in Foam::divide<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> (res=...,
gf1=..., gf2=...)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:564
#5 0x000000000045eec0 in Foam:perator/<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> (
tgf1=..., tgf2=...)
at /home/OpenFOAM/openfoam210/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:564
#6 0x0000000000446ae9 in main (argc=1, argv=0x7fffffffd3b8) at liftDragCoeffs.H:41
(gdb)

It looks like the culprit is at line #41 of liftDragCoeffs.H:

#6  0x0000000000446ae9 in main (argc=1, argv=0x7fffffffd3b8) at liftDragCoeffs.H:41

[enoch]

The line 41 is );, which is not the culprit. I think pow(magsymmGradUb, 0.5) has local zero values in some cells. How can deal with this situation?

line35 volVectorField shearLift
line36 (
line37 Cl*scalar(1)/(beta*alpha + scalar(0.001))
line38 *alpha*6.0*6.46/(4.0*Pi_kim)
line39 *pow(rhob, 0.5)*pow(nub*rhob, 0.5)/da
line40 *(symmGradUb & Ur)/pow(magsymmGradUb, 0.5)
line41 );

[kmooney]

Originally Posted by enoch

The line 41 is );, which is not the culprit. I think pow(magsymmGradUb, 0.5) has local zero values in some cells. How can deal with this situation?

line35 volVectorField shearLift
line36 (
line37 Cl*scalar(1)/(beta*alpha + scalar(0.001))
line38 *alpha*6.0*6.46/(4.0*Pi_kim)
line39 *pow(rhob, 0.5)*pow(nub*rhob, 0.5)/da
line40 *(symmGradUb & Ur)/pow(magsymmGradUb, 0.5)
line41 );

Try this:
pow(magsymmGradUb+VSMALL, 0.5)

If that's what you're trying to do in the 'beta*alpha + scalar(0.001)' expression you might want to use VSMALL there as well.

[enoch]

When I added VSMALL into magsymmGradUb, i.e, magsymmGradUb+VSMALL, I have a dimensional problem shown below.
How can I define dimesnions for VSMALL which is same as that of magsymmGradUb???


--> FOAM FATAL ERROR:
LHS and RHS of + have different dimensions
dimensions : [0 0 -1 0 0 0 0] + [0 0 0 0 0 0 0]


From function operator+(const dimensionSet&, const dimensionSet&)
in file dimensionSet/dimensionSet.C at line 514.

FOAM aborting

[kmooney]

I'm thinking something like this might work (might not be correct)...

dimensionedScalar mySmallScalar
(
"mySmallScalar",
dimensionSet(0,0,-1,0,0,0,0),
VSMALL
);

[enoch]

Dear Kyle Mooney,

Thanks for your speedy replies. It works with the new lift model in the code.
Believe it or not, I've spent about three months with this problem to add the additional force to the code. Thanks to you, I've learnt how to debug a code and a trick on how to deal with singularity. Thanks again. Thanks, OF users for sharing your experience with those who just start learn OF like me.

[aamo]

hi all,

I am simulating a backward step using piso and i have used periodic boundary condition in the span wise direction of the flow. I have modified the std ke model a bit, which is working when i used wall condition instead of periodic bc. here is the error which i ma getting..

"Time = 0.005

Courant Number mean: 1.39002e+47 max: 1.53594e+52
DILUPBiCG: Solving for Ux, Initial residual = 0.999997, Final residual = 0.999997, No Iterations 1001
DILUPBiCG: Solving for Uy, Initial residual = 0.999997, Final residual = 0.999997, No Iterations 1001
DILUPBiCG: Solving for Uz, Initial residual = 0.999996, Final residual = 0.999996, No Iterations 1001
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.0358884, No Iterations 3
time step continuity errors : sum local = 6.85443e+84, global = -4.19787e+68, cumulative = -4.19787e+68
DICPCG: Solving for p, Initial residual = 1, Final residual = 5.28645e-07, No Iterations 18
time step continuity errors : sum local = 2.01227e+99, global = 3.06191e+88, cumulative = 3.06191e+88
#0 Foam::error:rintStack(Foam::Ostream&) in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
#5 Foam::fvMatrix<double>::solve() in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/bin/pisoFoam"
#6 Foam::incompressible::RASModels::kEpsilon::correct () in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libincompressibleRASModels.so"
#7
in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/bin/pisoFoam"
#8 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#9
in "/home/yudhast/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/bin/pisoFoam"
Floating point exception"


I am very new to openfoam so please help