CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   OpenFOAM floating point Error (https://www.cfd-online.com/Forums/openfoam-programming-development/173072-openfoam-floating-point-error.html)

upuli June 13, 2016 02:21

OpenFOAM floating point Error
 
Dear members

I run a customized solver and when I run it in the designed case it gave the error shown below. I could not figure out what is the error in the code.


upuli@upuli-To-be-filled-by-O-E-M:~/tutorial/particle$ myreactingFoam
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.3.0-f5222ca19ce6
Exec : myreactingFoam
Date : Jun 13 2016
Time : 11:46:32
Host : "upuli-To-be-filled-by-O-E-M"
PID : 6856
Case : /home/upuli/tutorial/particle
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


Reading g
Reading physicalProperties

Reading field rho

Reading field p

Reading field U

Reading field k

Reading field epsilon

Reading field alphat

Reading/calculating face flux field phi

Creating field dpdt

Creating field kinetic energy K

Reading field G

Reading field Tg

Reading field Ts

Reading field mMoisture

Reading field mcellulose

Reading field mhemicellulose

Reading field mlignin

Reading field mWood

Reading field mchar

Reading field mash

Reading field YO2

Reading field YH2

Reading field YCO2

Reading field YCO

Reading field YH2O

Reading field YCH4

Reading field YN2

Courant Number mean: 1.66357e-05 max: 0.00231433

PIMPLE: Operating solver in PISO mode


Starting time loop

Courant Number mean: 1.66357e-05 max: 0.00231433
deltaT = 6e-05
Time = 6e-05

diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for Ts, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mWood, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mchar, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mash, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mMoisture, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 6.38914e-07, No Iterations 24
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 6.56505e-07, No Iterations 24
DILUPBiCG: Solving for Tg, Initial residual = 1, Final residual = 5.12893e-17, No Iterations 1
DILUPBiCG: Solving for YN2, Initial residual = 1, Final residual = 2.99007e-17, No Iterations 1
DILUPBiCG: Solving for YH2O, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for YCO2, Initial residual = 1, Final residual = 7.22578e-17, No Iterations 1
DILUPBiCG: Solving for YCO, Initial residual = 1, Final residual = 3.65492e-17, No Iterations 1
DILUPBiCG: Solving for YH2, Initial residual = 1, Final residual = 3.23385e-17, No Iterations 1
DILUPBiCG: Solving for YCH4, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for YO2, Initial residual = 1, Final residual = 3.48572e-17, No Iterations 1
DICPCG: Solving for p, Initial residual = 3.02161e-08, Final residual = 3.02161e-08, No Iterations 0
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0.701402, global = 0.701402, cumulative = 0.701402
DICPCG: Solving for p, Initial residual = 3.02161e-08, Final residual = 3.02161e-08, No Iterations 0
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0.701402, global = 0.701402, cumulative = 1.4028
ExecutionTime = 0.74 s ClockTime = 1 s

Courant Number mean: 3.72766e-30 max: 4.2478e-28
deltaT = 7.2e-05
Time = 0.000132

diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for Ts, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mWood, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mchar, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mash, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for mMoisture, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 9.99078e-07, No Iterations 74
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 9.99329e-07, No Iterations 74
DILUPBiCG: Solving for Tg, Initial residual = 1, Final residual = 2.83412e-17, No Iterations 1
DILUPBiCG: Solving for YN2, Initial residual = 2.68917e-08, Final residual = 2.68917e-08, No Iterations 0
DILUPBiCG: Solving for YH2O, Initial residual = 1, Final residual = 1.14458e-09, No Iterations 1
DILUPBiCG: Solving for YCO2, Initial residual = 1, Final residual = 8.29709e-10, No Iterations 1
DILUPBiCG: Solving for YCO, Initial residual = 1, Final residual = 7.88245e-10, No Iterations 1
DILUPBiCG: Solving for YH2, Initial residual = 1, Final residual = 7.8634e-10, No Iterations 1
DILUPBiCG: Solving for YCH4, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for YO2, Initial residual = 1, Final residual = 1.47638e-10, No Iterations 1
DICPCG: Solving for p, Initial residual = 1, Final residual = 1.8251e-16, No Iterations 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.41561e+16, global = 1, cumulative = 2.4028
DICPCG: Solving for p, Initial residual = 9.92675e-17, Final residual = 9.92675e-17, No Iterations 0
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 4.10304e+15, global = 1, cumulative = 3.4028
ExecutionTime = 1.48 s ClockTime = 1 s

Courant Number mean: 5.87408e-25 max: 8.96153e-21
deltaT = 8.64e-05
Time = 0.0002184

#0 Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OSspecific/POSIX/printStack.C:221
#1 Foam::sigFpe::sigHandler(int) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OSspecific/POSIX/signals/sigFpe.C:108
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 in "/lib/x86_64-linux-gnu/libm.so.6"
#4 pow in "/lib/x86_64-linux-gnu/libm.so.6"
#5 Foam::pow(double, double) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/doubleFloat.H:78
#6 Foam::pow(Foam::Field<double>&, Foam::UList<double> const&, double const&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/fields/Fields/scalarField/scalarField.C:119 (discriminator 2)
#7 void Foam::pow<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensioned<double> const&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricScalarField.C:274
#8
at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricScalarField.C:328
#9
at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricScalarField.C:350 (discriminator 1)
#10
at ~/mysolver/myreactingFoam/variables.H:374 (discriminator 1)
#11 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12
at ??:?
Floating point exception (core dumped)


when I look at the variable written in variables.H:374 it is given as below

volScalarField Nu
(
"Nu",
0.3+(0.62*pow(Re,0.5)*pow(Pr,0.33)*pow(1+pow(Re/282000,0.625),0.8)/pow(1+pow(0.4/Pr,0.66),0.25))
);


I couldnot find what is the error in the above written lines.Can please someone help me to fix it.

Thanks

A_Pete June 13, 2016 06:57

Look for any division by "0" in your whole code. What is the value of your Prandtl number Pr, for example?

Code:

Pout<< "Pr = " << Pr << endl;
Put "Info<<" or "Pout<<" statements at the positions of your code, where it aborts. Really look for any possible division by zero. Hard to give an answer without knowing what's in the code.

hk318i June 14, 2016 10:26

Generally speaking, in any division operation add a VSMALL to the denominator. It is a good practice.

Code:

volScalarField Nu
    (
      "Nu",
    0.3+(0.62*pow(Re,0.5)*pow(Pr,0.33)*pow(1+pow(Re/282000,0.625),0.8)
    /pow(1+pow(0.4/(Pr+VSMALL),0.66),0.25))
    );

Try this to check if Pr is causing this problem or not.

Bw,
Hassan

upuli June 15, 2016 22:45

Dear all
Thank you for the answers. I added small value to the denominator . But it still give the error in the same way

volScalarField Nu
(
"Nu",
0.3+(0.62*pow(Re,0.5)*pow(Pr,0.33)*pow(1+pow(Re/282000,0.625),0.8)/pow(1+pow(0.4/(Pr+Prsmall),0.66),0.25))
);

hk318i June 16, 2016 07:52

Hello,

I created a small program to test your formula. It works fine with VSMALL and doesn't work without it as expected. Also, I noticed that, this formula is really sensitive to Pr value and independent of Re value. Here is the code you can use it to try different values for Pr and Re. Regarding the error, it could be due to something else.

Bet wishes,
Hassan


Code:

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
    label size = 10;
    scalarField Pr(size, 0.0);
    scalarField Re(size, 0.0);
    scalarField Nu(size, 0.0);

    Nu = 0.3+(0.62*Foam::pow(Re,0.5)*Foam::pow(Pr,0.33)
        * Foam::pow(1.0+Foam::pow(Re/282000,0.625),0.8)
        / Foam::pow(1+Foam::pow(0.4/(Pr+VSMALL),0.66),0.25));

    Info << Pr << endl;
    Info << Nu << endl;

    Info<< "End\n" << endl;

    return 0;
}

// ************************************************************************* //


upuli June 20, 2016 03:19

Dear All
I was able to solve the error in Nu.But the following error appears


#0 Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OSspecific/POSIX/printStack.C:221
#1 Foam::sigFpe::sigHandler(int) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OSspecific/POSIX/signals/sigFpe.C:108
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/fields/Fields/scalarField/scalarField.C:98 (discriminator 2)
#4 double Foam::gSumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&, int) at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/FieldFunctions.C:560
#5 Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ~/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C:138 (discriminator 1)
#6 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C:169 (discriminator 3)
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:82
#8
at ~/mysolver/myreactingFoam/pEqn.H:76 (discriminator 1)
#9 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10
at ??:?
Floating point exception (core dumped)
My pEqn is

volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhog*rAU));

volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();

if (pimple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
(fvc::interpolate(rhog*HbyA) & mesh.Sf())
+ rhorAUf*fvc::ddtCorr(rhog, U, phi)
)/fvc::interpolate(rhog)
);
surfaceScalarField phie
(
"phie",
fvc::interpolate(solidporosity)*phid
);



while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(solidporosity*psi, p)
+ fvm::div(phie, p)
- fvm::laplacian(solidporosity*rhog*rAU, p)
==
gasgeneration
);



pEqn.solve(mesh.solver(p.select(pimple.finalInnerI ter())));

if (pimple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(rhog*HbyA) & mesh.Sf())
+ rhorAUf*fvc::ddtCorr(rhog, U, phi)
)
);

surfaceScalarField phiHbyB
(
"phiHbyB",
fvc::interpolate(solidporosity)*phiHbyA
);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(solidporosity*psi, p)
+ fvc::div(phiHbyB)
- fvm::laplacian(solidporosity*rhog*rAU, p)
==
gasgeneration
);

pEqn.solve(mesh.solver(p.select(pimple.finalInnerI ter())));

if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyB + pEqn.flux();
}
}
}

#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"

U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();

K = 0.5*magSqr(U);


All times are GMT -4. The time now is 10:08.