CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

InterFoam based solver running into floating point error on restarting simulation

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 21, 2021, 17:02
Default InterFoam based solver running into floating point error on restarting simulation
  #1
Member
 
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6
Venky_94 is on a distinguished road
I'm running a bubble growth simulation on a custom solver based on interFoam, and I'm facing a weird issue. My code is running well the first time, but if I try to restart the simulation later from an earlier point, it doesn't start up and runs into a floating point exception.

In this particular case, the job ran till around 2.5s and the job exited because I had set a 24 hr timer on the cluster. And later I tried to restart the simulation, it kept throwing a floating point exception. And the weird thing is it kept happening even when I restarted from earlier save points of 0.1s.

I'm not sure how to debug this. Any support would be appreciated.

Code:
Reading field p_rgh

Reading field U

Using dynamicCode for codedFixedValue parabolicinlet at line 112370 in "/scratch/ganeshvt/VOF_Nodes1/0.1/U/boundaryField/airinlet"
Reading/calculating face flux field phi

Reading transportProperties

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
Selecting surfaceTensionModel constant
No phase change: phaseChangeProperties not found
Selecting phaseChange model none
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib64/libc.so.6"
#3  ? in "/lib64/libm.so.6"
#4  Foam::pow(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#5  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::pow<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#6  ? at ??:?
#7  __libc_start_main in "/lib64/libc.so.6"
#8  ? at ??:?
Floating point exception
Venky_94 is offline   Reply With Quote

Old   November 22, 2021, 10:15
Default
  #2
Senior Member
 
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17
matejfor is on a distinguished road
What is the first field it is reading? Apparently it is failing reading this field, is it written at the end correctly? Have you tried to restart from a different time?
matejfor is offline   Reply With Quote

Old   November 22, 2021, 11:01
Default
  #3
Member
 
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6
Venky_94 is on a distinguished road
Quote:
Originally Posted by matejfor View Post
What is the first field it is reading? Apparently it is failing reading this field, is it written at the end correctly? Have you tried to restart from a different time?
Thanks for you response.
  1. How could I check which field is being read in first?
  2. Yes, the fields are written out correctly and I'm even able to post process the case upto the time it has run without any issues.
  3. And no, I'm unable to restart the simulation from any other time step.
Venky_94 is offline   Reply With Quote

Old   November 22, 2021, 12:04
Default
  #4
Senior Member
 
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17
matejfor is on a distinguished road
well maybe not a field, as U and P_rgh are already read. What does in your transportProperties file contain? Last info it has read correctly is: "Selecting phaseChange model none".
How does the file looks like below it? Is the code trying to read the entry using correct read function? Meaning reading what it is supposed to read?

Show us the source code and the transportProperties file.
matejfor is offline   Reply With Quote

Old   November 22, 2021, 12:57
Default
  #5
Member
 
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6
Venky_94 is on a distinguished road
Quote:
Originally Posted by matejfor View Post
well maybe not a field, as U and P_rgh are already read. What does in your transportProperties file contain? Last info it has read correctly is: "Selecting phaseChange model none".
How does the file looks like below it? Is the code trying to read the entry using correct read function? Meaning reading what it is supposed to read?

Show us the source code and the transportProperties file.
I'm attaching the transportProperties file and the link to the solver source. I'm using OpenFoam 9 btw.

Solver: https://github.com/Venky-94/interFoamMod
Attached Files
File Type: zip transportProperties.zip (512 Bytes, 3 views)
Venky_94 is offline   Reply With Quote

Old   November 22, 2021, 13:48
Default
  #6
Senior Member
 
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17
matejfor is on a distinguished road
OK, the code in createFields.H is interesting. After reading transportProperties, you are trying to read in the density.

Code:
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
const volScalarField& nuModel1 = mixture.nuM1();
const volScalarField& nuModel2 = mixture.nuM2();


// Need to store rho for ddt(rho, U)
volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
		IOobject::AUTO_WRITE
    ),
	rho2 + (rho1 - rho2)*pow(alpha1, 1.5*(nuModel1/nuModel2) + 0.75)
	//alpha1*rho1 + alpha2*rho2
);
On the initial run, do you have a density file, or is it calculated?
Density should be stored as there is an AUTO_WRITE so it will be written as all the others. If not, density is calculated.

The equation is: scalar2 + (scalar1 - scalar2)*pow(volScalarField, (1.5*volScalarField/volScalarField) + 0.75)



I would bet, looking at the error message, someone does not like the arguments of the pow function. Is the density actually written down?

It makes sense to have some Info statements among constructors to easily identify what goes wrong.

Hope this helps.
matejfor is offline   Reply With Quote

Old   November 22, 2021, 20:00
Default
  #7
Member
 
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6
Venky_94 is on a distinguished road
Quote:
Originally Posted by matejfor View Post
OK, the code in createFields.H is interesting. After reading transportProperties, you are trying to read in the density.

Code:
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
const volScalarField& nuModel1 = mixture.nuM1();
const volScalarField& nuModel2 = mixture.nuM2();


// Need to store rho for ddt(rho, U)
volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
		IOobject::AUTO_WRITE
    ),
	rho2 + (rho1 - rho2)*pow(alpha1, 1.5*(nuModel1/nuModel2) + 0.75)
	//alpha1*rho1 + alpha2*rho2
);
On the initial run, do you have a density file, or is it calculated?
Density should be stored as there is an AUTO_WRITE so it will be written as all the others. If not, density is calculated.

The equation is: scalar2 + (scalar1 - scalar2)*pow(volScalarField, (1.5*volScalarField/volScalarField) + 0.75)



I would bet, looking at the error message, someone does not like the arguments of the pow function. Is the density actually written down?

It makes sense to have some Info statements among constructors to easily identify what goes wrong.

Hope this helps.
That actually makes a lot of sense. I remember that I added the AUTO_WRITE option to it because I needed the density field for post-processing.

My density field is being written out at every time step now, and is available in the reconstructed time step folder as well. Any inputs on how to rectify this issue?
Venky_94 is offline   Reply With Quote

Old   November 23, 2021, 14:54
Default
  #8
Member
 
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6
Venky_94 is on a distinguished road
So I changed the expression back to original (as in interFoam) in the createFields dictionary, and the simulation restarts without any error.

However I'm new to OpenFoam and I just want to understand something here. The expression here is only important when the simulation is starting up right? i.e. Density would be calculated as per the expression contained in alphaEqnSubCycle dictionary from the second iteration onwards and the expression within createFields dictionary wouldn't influence the simulation.

Code:
// createFields.H

Need to store rho for ddt(rho, U)
volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
		IOobject::AUTO_WRITE
    ),
	//rho2 + (rho1 - rho2)*pow(alpha1, 1.5*(nuModel1/nuModel2) + 0.75)
	alpha1*rho1 + alpha2*rho2
);
Code:
alphaEqnSubCycle.H

rho == (rho2 + (rho1 - rho2)*pow(limitedAlpha1, 1.5*(nuModel1/nuModel2) + 0.75));
Venky_94 is offline   Reply With Quote

Old   November 23, 2021, 17:21
Default
  #9
Senior Member
 
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17
matejfor is on a distinguished road
I'll start from a broader look. The
Code:
volScalarField rho (argument1, agument2);
is a constructor of the density field.
One of the predefined constructors for volField is having the first argument as IOobject.
the second argument is either mesh object - then the field is read from the file located in runTime.timeName() time directory. If it is an expression it is evaluated based on the expression. In your case the READ_IF_PRESENT allows to read the values if the rho files is present in the time dir you are starting from, otherwise the expression is used. As you have said previously.

This is, however, only happening at the start, as it is part of the createFields.H file, which is located outside of the solution loop and is used only at the initialisation to set the initial density field.

The values are then overwritten by the alphaEqnSubCycle.H as the alpha equation is solved the first thing at every time step, including the first one. So every iteration the density rho - which is used only for the time derivative in the momentum equation is updated correctly and there is no need to initialise it with your strange power equation.

As for the error itself, it is the SIGFPE error - floating point exception, meaning in most of the cases some fancy bold operation like division by zero. So I§d make sure your nuModel2 is initially far from zero in every cell and then you can return to your power expression.
matejfor is offline   Reply With Quote

Old   November 23, 2021, 17:53
Default
  #10
Member
 
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6
Venky_94 is on a distinguished road
Thanks a lot for your wonderful explanation.

While I was reading your last paragraph about making sure nuModel2 is far from zero before returning to the expression, I realized what was causing the error. I was so focused on the viscosity values that I forgot to check on the other argument for the pow expression, which is the alpha.

I realized that the reason the expression was working fine within alphaEqnSubCycle dictionary and not within the constructor was that I had this below expression within the dictionary to avoid negative values of alpha, whereas I was using alpha values as it is in the createFields dictionary.

Code:
const volScalarField limitedAlpha1
    (
        "limitedAlpha1",
        min(max(alpha1, scalar(0)), scalar(1))
    );
So I moved that expression to the createFields dictionary now and my solver is starting up and running fine as well. Silly me

Thanks a ton for your support throughout.
Venky_94 is offline   Reply With Quote

Reply

Tags
floating point exception, interfoam, print stack, sigfpe

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
ERORR floating point exception fluent udf cudau.95 Fluent UDF and Scheme Programming 0 August 4, 2021 23:11
[blockMesh] Floating Point Exception while generating wedge based mesh jns-v OpenFOAM Meshing & Mesh Conversion 9 July 8, 2021 06:36
Problem with an old Simulation FrankW CFX 3 February 8, 2016 05:28
Hypersonic Sim (running at Mach 7): floating point exception has occurred bungusbeefcake STAR-CCM+ 10 March 31, 2015 07:59
[snappyHexMesh] snappyHexMesh and cyclic boundaries Ruli OpenFOAM Meshing & Mesh Conversion 2 December 9, 2013 07:51


All times are GMT -4. The time now is 02:20.