[faab]

Dear FOAMers,

I am trying to solve the temperature-induced pressure gradient in liquid helium (confined in a channel) occurring from a discharge of high energy. For this, I need to solve for the coupled momentum and energy equations and compute the thermodynamic properties using an external fortran module. For this purpose, I am building a compressible solver using the chtMultiRegionFoam solver as a reference, and am now able to run the solver but it keeps crashing at the time step "0.107" (see error from the log hereunder) no matter what parameters I set in the fvSolution and fvSchemes dictionaries:

 
Region: hell_channel Courant Number mean: 3.6443497e-06 max: 6.7134329e-05
Time = 0.117

Solving for superfluid helium region hell_channel
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.012916636, Final residual = 2.1946216e-14, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0080626801, Final residual = 9.8190616e-15, No Iterations 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.00011588099, Final residual = 9.1147225e-17, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 8.5073936e-05, Final residual = 1.2611083e-16, No Iterations 1
DILUPBiCG:  Solving for hi, Initial residual = 2.5778752e-08, Final residual = 2.5778752e-08, No Iterations 0
Min/max T:3.0016428 3.0445642
GAMG:  Solving for p, Initial residual = 0.0002906377, Final residual = 3.8736628e-09, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (hell_channel): sum local = 8.5116225e-07, global = 4.9213372e-13, cumulative = -0.0010918572
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib64/libc.so.6"
#3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PCG::solve(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>::solve(Foam::dictionary const&) at ~/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:82
#10  
 at ~/OpenFOAM/OpenFOAM-2.3.0/applications/solvers/heatTransfer/CompressibleHeSolver/./fluid/pEqn.H:78 (discriminator 1)
#11  __libc_start_main in "/lib64/libc.so.6"
#12  
 at ??:?
Floating point exception

The error comes from the following equation in the code:

 
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
        {
            fvScalarMatrix pEqn
            (
                pDDtEqn
              - fvm::laplacian(rhorAUf, p)
            );
            pEqn.solve
            (
                mesh.solver
                (
                    p.select
                    (
                        (
                           oCorr == nOuterCorr-1
                        && corr == nCorr-1
                        && nonOrth == nNonOrthCorr
                        )
                    )
                )
            );

And the pressure is computed using the following sets:

 
p
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        smoother         GaussSeidel;
        cacheAgglomeration  true;
        nCellsInCoarsestLevel 60;
        agglomerator     faceAreaPair;
        mergeLevels      1;
    }
   
    pFinal
    {
        $p;
        tolerance        1e-7;
        relTol           0;
    }

The boundary conditions to the external boundaries of my domain are zeroGradient because the liquid is confined in a volume, the system is then closed. And the log for the previous time step before the simulation crashes does not seem to indicate that the solution and/or residuals are diverging:

 
Region: hell_channel Courant Number mean: 3.5834383e-06 max: 6.7071608e-05
Time = 0.116

Solving for superfluid helium region hell_channel
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.013165873, Final residual = 2.1176503e-14, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0081476332, Final residual = 9.3894979e-15, No Iterations 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.00011809902, Final residual = 8.8268892e-17, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 8.5991334e-05, Final residual = 1.2143285e-16, No Iterations 1
DILUPBiCG:  Solving for hi, Initial residual = 2.5511043e-08, Final residual = 2.5511043e-08, No Iterations 0
Min/max T:3.0016456 3.0448337
GAMG:  Solving for p, Initial residual = 0.00027577967, Final residual = 2.4234806e-10, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (hell_channel): sum local = 8.1632713e-07, global = -3.0635965e-12, cumulative = -0.0010918571
GAMG:  Solving for p, Initial residual = 1.9836363e-06, Final residual = 2.2204109e-11, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (hell_channel): sum local = 8.1632298e-07, global = -3.5110486e-11, cumulative = -0.0010918572
Solving for solid region solid1
DICPCG:  Solving for T, Initial residual = 0.0027234834, Final residual = 4.5801085e-08, No Iterations 6
Min/max T:min(T) [0 0 0 1 0 0 0] 3.0018298 max(T) [0 0 0 1 0 0 0] 29.987348
Solving for solid region solid2
DICPCG:  Solving for T, Initial residual = 0.0026282575, Final residual = 4.0420011e-08, No Iterations 8
Min/max T:min(T) [0 0 0 1 0 0 0] 3.0019453 max(T) [0 0 0 1 0 0 0] 29.985524
ExecutionTime = 16.67 s  ClockTime = 17 s

I have tried everything I could think of and really need some insight, if anyone could help me that would be highly appreciated.
Thank you all !