|
[Sponsors] |
InterFoam based solver running into floating point error on restarting simulation |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
November 21, 2021, 17:02 |
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 |
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 |
|
November 22, 2021, 10:15 |
|
#2 |
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17 |
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?
|
|
November 22, 2021, 11:01 |
|
#3 | |
Member
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6 |
Quote:
|
||
November 22, 2021, 12:04 |
|
#4 |
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17 |
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. |
|
November 22, 2021, 12:57 |
|
#5 | |
Member
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6 |
Quote:
Solver: https://github.com/Venky-94/interFoamMod |
||
November 22, 2021, 13:48 |
|
#6 |
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17 |
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 ); 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. |
|
November 22, 2021, 20:00 |
|
#7 | |
Member
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6 |
Quote:
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? |
||
November 23, 2021, 14:54 |
|
#8 |
Member
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6 |
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)); |
|
November 23, 2021, 17:21 |
|
#9 |
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17 |
I'll start from a broader look. The
Code:
volScalarField rho (argument1, agument2); 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. |
|
November 23, 2021, 17:53 |
|
#10 |
Member
Venkat Ganesh
Join Date: May 2020
Location: Cincinnati, Ohio
Posts: 49
Rep Power: 6 |
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)) ); Thanks a ton for your support throughout. |
|
Tags |
floating point exception, interfoam, print stack, sigfpe |
Thread Tools | Search this Thread |
Display Modes | |
|
|
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 |