| A_Pete |
September 23, 2013 09:37 |
Floating point exceptions with own modified k-omega SST with Reboud's correction
Hi everyone!
I have implemented the correction of turbulent viscosity by Reboud in interface regions for cavitating flows. Physically it does everything I want pretty well, but from time to time my solutions will abort due to floating point exception errors.
For me this is an indication that I didn't prevent the solver from dividing by zero somewhere. So I found a spot where this might happen during the solution, altough it shouldn't and implemented a small prevention.
But that still doesn't get the job done, because I will still get the same FPEs and have no idea where to look for. This is often accompanied by omega getting very large first. But since this doesn't happen with other turbulence models in my simulations something has to be wrong with my code.
What I have done is take the standard kOmegaSST turbulence model of OpenFOAM and modify it the way I need it.
The error says this:
Code:
[0] #0 Foam::error::printStack(Foam::Ostream&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #1 Foam::sigFpe::sigHandler(int) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #2 in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #3 double Foam::sumProd<double>(Foam::UList<double> const&,
Foam::UList<double> const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #4 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&,
unsigned char) const in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[0] #5 Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
[0] #6 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #7 Foam::fvMatrix<double>::solve() in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #8 Foam::SolverPerformance<double>
Foam::solve<double>(Foam::tmp<Foam::fvMatrix<double> > const&) in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #9 Foam::incompressible::RASModels::modKOmegaSST::correct() in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libincompressibleRASModels.so"
[0] #10
[0] in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[0] #11 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #12
[0] in
"/home/peters/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/interPhaseChangeDyMFoam"
[prop-04:28417] *** Process received signal ***
[prop-04:28417] Signal: Floating point exception (8)
...
Something seems to go wrong at #9 in the correct() method of modKOmegaSST, but I can't find it. Below are what I suppose to be the interesting parts of my code regarding this issue.
In the header file I added this:
Code:
// under member functions
//- Calculate Reboud's correction and return the normalized density function
virtual volScalarField normDensFunc();
And in the source file the interesting parts are the following:
Code:
// definition of the calculated function under member functions
volScalarField modKOmegaSST::normDensFunc()
{
//- Reading fields for Reboud's correction
const volScalarField& alpha1 = mesh_.lookupObject<volScalarField>("alpha1");
const volScalarField& alpha2 = mesh_.lookupObject<volScalarField>("alpha2");
//- Lookup transportProperties dictionary to recieve fluid properties
const dictionary& transportProperties =
mesh_.lookupObject<IOdictionary>("transportProperties");
const dictionary& phase1 = transportProperties.subDict("phase1");
const dictionary& phase2 = transportProperties.subDict("phase2");
const dimensionedScalar& rho1 = phase1.lookup("rho");
const dimensionedScalar& rho2 = phase2.lookup("rho");
//- Calculation of mixture density
const volScalarField rho_m = alpha1*rho1 + alpha2 * rho2;
//- Calculation of density function
const volScalarField f_rho =
pow(( rho_m - rho2 ), nReb_) / pow((rho1-rho2), (nReb_-1)) + rho2;
//- Calculation of ratio of density function and mixture density
const volScalarField fRhoNorm =
f_rho / max(rho_m, rho1);
//const volScalarField& test = mesh_.lookupObject<volScalarField>("alpha1");
return fRhoNorm;
}
...
void modKOmegaSST::correct()
{
...
// function is called and used for calculation of nut_ below
//- Get normalized density function for Reboud's correction
const volScalarField& fRhoNorm_(normDensFunc());
// Re-calculate viscosity; standard calculation is the same, just w/o fRhoNorm
nut_ = a1_*k_*fRhoNorm_/max(a1_*omega_, b1_*F23()*sqrt(S2));
nut_.correctBoundaryConditions();
}
So if anybody has any idea what I could try out I would be very happy and thankful.
Thanks in advance.
Best regards,
Andreas |
