[JulioPieri]

Dear Foamers,

I'm facing a very weird and enigmatic behavior in my transient twoPhaseEulerFoam simulation of a fluidized bed. I'm using OF v2012, and got the same result with v1812, both on WSL Ubuntu. I've been working on a series of similar cases for over 4 months now, and I have had pretty good results so far.

During a power shortage the other night, some cases stopped running on the server after 7.45s of simulation (it took 3 days to get to this point, hence why I want to solve this issue). This have happened before and I could easily restart them. However, this particular setup won't restart, always diverging on the second time iteration - no matter if I deleted one, two, ten, whatever, latest timesteps. Smaller timesteps also didn't help.

I've created a simple case that reproduces this. If I leave it be, fresh start, it runs smoothly with no problems. If I interrupt it anywhere, though, it doesn't restart. Even stranger: if I write every timestep calculated, I can advance with the solution:
it writes the first iteration -> diverges on second tstep -> I resume it -> it writes one more -> diverges again, etc. Looking at paraview, the solution evolves as I would expect, suggesting the setup is ok.

The error is the following:

PIMPLE: iteration 1
MULES: Solving for alpha.particles
MULES: Solving for alpha.particles
smoothSolver:  Solving for alpha.particles, Initial residual = 3.29388e-06, Final residual = 2.35206e-10, No Iterations 1
alpha.particles volume fraction = 0.31103  Min(alpha.particles) = 0  Max(alpha.particles) = 0.600047
Constructing momentum equations
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libpthread.so.0
#3  Foam::tmp<Foam::GeometricField<Foam::typeOfMag<Foam::Vector<double> >::type, Foam::fvPatchField, Foam::volMesh> > Foam::mag<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#4  Foam::phasePair::magUr() const at ??:?
#5  Foam::phasePair::Re() const at ??:?
#6  Foam::heatTransferModels::RanzMarshall::K() const at ??:?
#7  Foam::BlendedInterfacialModel<Foam::heatTransferModel>::K() const at ??:?
#8  Foam::twoPhaseSystem::Kh() const at ??:?
#9  ? at ??:?
#10  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#11  ? at ??:?
Floating point exception

I traced back the error to the function magUr() inside phasePair.C, which simply is:

Foam::tmp<Foam::volScalarField> Foam::phasePair::magUr() const
{
    return mag(phase1().U() - phase2().U());
}

Therefore, no obvious invalid math operation.

If I change the Heat Transfer Model from RanzMarshall to the other available, called spherical, the error changes and get closer to what I believe might be related to the cause:

PIMPLE: iteration 1
MULES: Solving for alpha.particles
MULES: Solving for alpha.particles
smoothSolver:  Solving for alpha.particles, Initial residual = 3.44837e-06, Final residual = 2.4905e-10, No Iterations 1
alpha.particles volume fraction = 0.31103  Min(alpha.particles) = 0  Max(alpha.particles) = 0.600048
Constructing momentum equations
smoothSolver:  Solving for e.particles, Initial residual = 8.98554e-08, Final residual = 2.86103e-08, No Iterations 1
smoothSolver:  Solving for e.air, Initial residual = 0.0395403, Final residual = 3.89446e-07, No Iterations 2
min T.particles 1122.95
min T.air 426
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libpthread.so.0
#3  Foam::scalarProduct<double, double>::type Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PCG::scalarSolve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<double>&, Foam::Field<double> const&) const at ??:?
#6  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#7  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#8  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#9  Foam::fvMatrix<double>::solveSegregatedOrCoupled(Foam::dictionary const&) at ??:?
#10  Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const at ??:?
#11  ? at ??:?
#12  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#13  ? at ??:?
Floating point exception

This pointed me to a possible problem with the pressure step of the PIMPLE algorithm. I then decided to run the case with Floating Point Exception trapping desabled:

export FOAM_SIGFPE=false

Which pointed again at possible pressure related problem:

[...] at second timestep:
Constructing momentum equations
smoothSolver:  Solving for e.particles, Initial residual = 8.98554e-08, Final residual = 2.86103e-08, No Iterations 1
smoothSolver:  Solving for e.air, Initial residual = 0.0395403, Final residual = 3.89446e-07, No Iterations 2
min T.particles 1122.95
min T.air 426
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = nan, No Iterations 1000
GAMG:  Solving for p_rgh, Initial residual = nan, Final residual = nan, No Iterations 1000

For some reason, the solver seems to not be able to carry the information over the second timestep, but only when the simulation is restarted. It seems to be initializing p_rgh as inf or NaN (therefore the Initial residual = 1). The first timestep might work because it reads from the HD instead of memory (my guess here).

So far, what I've tried:

• Changing every single parameter in my PIMPLE subDict - didn't help.
• Changing timeStep size (fixed or adjustable) to very small ones -didn't help.
• Changing from ascii to binary and increasing (alot) writePrecision - didn't help
• Coarsen/refine the mesh - refining the mesh got it worse (diverges even from a fresh start), while coarsening it made me able to restart once the simulation. If I interrupt it again, it doesn't restart.
• Increase/descrese velocity - didn't help.
• Change multiple constant-folder properties, models and parameters - didn't help.
• Not running setFields and starting the domain without particle phase - did run, and was able to restart. Maybe because of the absence of one phase.
• Visually inspect the field files saved - it all seems ok. No odd number or incipit-divergence stored
• Running with FOAM_SETNAN=true - didn't help.
• Change some parameters in GAMG solver, or changing the precoditioner - couldn't even start simulation, actually.

I've tried eveything I can think of, and I can't manage to overcome the problem. I believe I might be close to a possible OF bug: this might be a very very specific set of parameters - mesh, timescale, velocity, phase fraction, etc - that makes the solver go crazy. Remember, I've had several other successful cases running with very similar parameters.

And if I fresh start it, it goes well!! The difference from this case to the last one that was ok is that this one have a slightly smaller domain (same element size).

Attached you can find the compressed case (zip or tar). To reproduce the error, hit Allrun, stop after it writes the first timestep (0.001s), then hit "twoPhaseEulerFoam" to restart it. You can change the heat transfer model inside constant/phaseProperties.

Any help or insight is welcome.
Thank you very much!