[Venky_94]

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.

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

[matejfor]

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?

[Venky_94]

Originally Posted by matejfor

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.

[matejfor]

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.

[Venky_94]

Originally Posted by matejfor

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

[matejfor]

OK, the code in createFields.H is interesting. After reading transportProperties, you are trying to read in the density.

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.

[Venky_94]

Originally Posted by matejfor

OK, the code in createFields.H is interesting. After reading transportProperties, you are trying to read in the density.

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]

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.

// 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
);

alphaEqnSubCycle.H

rho == (rho2 + (rho1 - rho2)*pow(limitedAlpha1, 1.5*(nuModel1/nuModel2) + 0.75));

[matejfor]

I'll start from a broader look. The

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.

[Venky_94]

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.

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.