OpenFOAM Floating Point Exception: Causes and Fix

A floating point exception (SIGFPE) in OpenFOAM means a field value has reached NaN or infinity. The solver crashes with a message like Floating point exception (core dumped) or #0 Foam::error::printStack. This is always a symptom — the underlying cause is one of a small set of known issues.

By Rayan Bhar, CFD Engineer · Updated July 2026

What the error actually means

OpenFOAM enables floating point exception trapping at runtime (controlled by the FOAM_SIGFPE environment variable, on by default). When any arithmetic produces NaN (division by zero, 0/0, sqrt of negative) or Inf (overflow), the OS sends a SIGFPE signal, OpenFOAM's sigFpe handler catches it, and the solver terminates immediately. The stack trace in the log tells you where the crash happened but not why.

A typical FPE crash log looks like this:

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  Foam::symGaussSeidelSmoother::smooth(...) at ??:?
Floating point exception (core dumped)

The crash inside symGaussSeidelSmoother means the smoother encountered NaN while solving the linear system — not that the smoother caused the problem. The NaN originated elsewhere and propagated into the solver. The diagnostic strategy is to find which field went NaN first, at which cell, and at which time step.

Reading the FPE block in your log: sigFpe, trapFpe and the stack trace

Three strings show up around floating point exceptions in OpenFOAM logs, and one of them fools almost everyone at least once:

trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE). is NOT an error. This line is printed in the header of every OpenFOAM log, on every healthy run, right after the case banner. It just announces that FPE trapping is switched on. If you are grepping a log for the word "exception" to find your crash, this banner is the false positive you will hit first — a run whose only "exception" line is the trapFpe banner did not crash at all.

Foam::sigFpe::sigHandler(int) is the real thing. When an actual FPE fires, the signal handler appears in the stack trace, typically like this (real log from a k-epsilon run blowing up at startup):

bounding k, min: -6.34e+74 max: 8.46e+78 average: 3.01e+75
[stack trace]
=============
#1  Foam::sigFpe::sigHandler(int) in libOpenFOAM.so
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  Foam::PCG::scalarSolve(...) in libOpenFOAM.so
#4  Foam::GAMGSolver::solveCoarsestLevel(...) in libOpenFOAM.so
=============
Floating point exception (core dumped)

Note the bounding k line just above the trace: k had already reached 10⁷⁴ before the crash. The bounding messages that precede a sigFpe are usually a better clue than the stack trace itself — they name the field that ran away and when.

Two debugging switches worth knowing: unsetting FOAM_SIGFPE lets the solver keep running through NaN instead of crashing (useful only to write out the poisoned fields for inspection — the results are garbage). Setting FOAM_SETNAN=true fills uninitialised memory with NaN so that any use of an uninitialised field crashes instantly instead of silently reading zeros — the sharpest tool for catching a missing field initialisation.

Startup crash vs mid-run FPE: two different diseases

Before hunting causes, check one thing in the log: did the solver complete any time steps before the FPE? The answer splits the diagnosis in two:

Zero time steps completed = a setup failure, not a numerical one. If the FPE (or a FOAM FATAL ERROR) appears during startup — while the log still shows "Create mesh", "Reading field p", "Reading transportProperties" — nothing diverged. The solver never got to solve anything. The cause is configuration: a missing or misnamed dict, a solver that doesn't match the case structure (for example a case built for one OpenFOAM fork run with a solver from the other), an invalid mesh, or a zero initialisation that poisons the very first evaluation. Courant number advice is useless here; read the startup error message itself, it names the missing piece.

N time steps completed, then FPE = a numerical blow-up. The run was viable and something accumulated: Courant number creeping up, unbounded scheme overshoots, turbulence quantities going negative. This is where sections 1, 2 and the bounding-message trail apply. The iteration count at the crash is itself a clue: a crash at iteration 3 points at initial conditions; a crash at iteration 400 after hours of stable running points at a slowly growing instability or an event (adaptive time step jump, a source term activating).

1. deltaT too large — Courant number exceeds 1

The most common cause. If the Courant number Co = U·Δt/Δx >> 1, the explicit time-stepping scheme becomes unstable. A single cell gets a velocity value that is unphysically large, which propagates and causes NaN within a few iterations.

Check the log: look for Courant Number max: just before the crash. If it's above 1 (or above 5–10 for PIMPLE), reduce deltaT in controlDict.

// controlDict — use adjustable time stepping
adjustTimeStep  yes;
maxCo           0.8;
maxDeltaT       1e-4;

Identifying the Courant number in the log

pimpleFoam and icoFoam print the Courant number at every time step before solving. Look for the output just before the crash:

Time = 0.0023

Courant Number mean: 0.342  max: 3.87

PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.231, ...
Floating point exception (core dumped)

A max Courant of 3.87 is the smoking gun. Reduce deltaT by a factor of 5 (here, from whatever it is to approximately 1e-5) and rerun. With adjustTimeStep yes and maxCo 0.8, the solver will maintain stability automatically.

2. Turbulence fields initialised at zero or negative

Setting k 0 or epsilon 0 in 0/k or 0/epsilon causes immediate division by zero in the turbulence model source terms. This is the second most common cause of FPE at t=0.

• Set k to a realistic estimate: 1.5 * (U * I)² with I = 0.05 (5%)
• Set epsilon from the mixing length: C_mu^0.75 * k^1.5 / L
• For omega: epsilon / (C_mu * k) with C_mu = 0.09
• Never use internalField uniform 0 for turbulence fields

The correct initialisation for a typical pipe flow at U=5 m/s with 5% turbulence intensity and D=0.05m hydraulic diameter:

// k = 1.5 * (5 * 0.05)^2 = 0.09375 m²/s²
// L = 0.07 * 0.05 = 0.0035 m
// epsilon = 0.09^0.75 * 0.09375^1.5 / 0.0035 = 0.765 m²/s³
// omega = 0.765 / (0.09 * 0.09375) = 90.7 s^-1

// 0/k
internalField   uniform 0.09375;

// 0/epsilon
internalField   uniform 0.765;

// 0/omega (for k-omega SST)
internalField   uniform 90.7;

Negative k values during simulation

Some k-epsilon and k-omega models can develop negative k values in regions of strong flow deceleration or adverse pressure gradient. While physically impossible (k is a variance and cannot be negative), the numerical scheme can produce it. OpenFOAM provides options to bound k and prevent FPE:

// In fvSolution, under SIMPLE or PIMPLE:
limitTurbulenceViscosity  yes;

The limitTurbulenceViscosity flag clips nut to prevent runaway turbulent viscosity amplification, which is a common secondary cause of FPE in turbulent simulations.

3. Negative cell volumes or degenerate mesh

A mesh with negative cell volumes (produced by a failed blockMesh or snappyHexMesh run) will immediately cause FPE when the solver tries to compute face fluxes. Run checkMesh before starting the solver:

checkMesh -latestTime

Look for: ***Zero or negative cell volumes or ***Zero or negative face areas. If present, the mesh must be rebuilt — there is no numerical fix for an invalid mesh.

Diagnosing snappyHexMesh failures

snappyHexMesh can produce degenerate cells near complex geometry, particularly at sharp feature edges or when refinement levels create very poor aspect ratio cells. The symptoms in checkMesh are:

***Zero or negative cell volumes detected, min volume = -2.34e-12
***Max aspect ratio = 45231
***Number of faces with skewness > 4: 23

To fix, tighten quality controls in snappyHexMeshDict:

meshQualityControls
{
    maxNonOrtho         65;
    maxBoundarySkewness 20;
    maxInternalSkewness 4;
    maxConcave          80;
    minVol              1e-13;
    nSmoothScale        4;
    errorReduction      0.75;
}

Stricter quality controls in snappyHexMesh prevent degenerate cells at the cost of slightly less surface conformance. This is almost always the right trade-off for solver stability.

4. Boundary condition inconsistency

A fixedValue with an unrealistic magnitude (e.g., velocity in the wrong units, or pressure in Pa when the solver expects m²/s²) can cause FPE after a few iterations when the solver amplifies the error. Double-check:

• Velocity U: in m/s, not km/h or mph
• Pressure p: in m²/s² (kinematic) for incompressible solvers, not Pa
• Temperature T: in Kelvin, not Celsius
• All patches in 0/ must be listed in the mesh boundary file

Unit conversion for incompressible pressure

This is one of the most common FPE causes for engineers coming from ANSYS Fluent, where pressure is always in Pa. In OpenFOAM's incompressible solvers (simpleFoam, pimpleFoam), pressure is kinematic: p = P/rho, where P is the physical pressure in Pa and rho is the density in kg/m³. For air at standard conditions:

// Physical pressure difference: 100 Pa
// rho_air = 1.225 kg/m³
// Kinematic pressure = 100 / 1.225 = 81.6 m²/s²

// 0/p — correct for incompressible solver
internalField   uniform 0;  // gauge pressure, reference = 0
inlet
{
    type    fixedValue;
    value   uniform 81.6;  // NOT 100 Pa
}

Setting uniform 100 in an incompressible case (thinking it's 100 Pa) tells the solver the kinematic pressure is 100 m²/s² — equivalent to a 122.5 Pa pressure difference for air. On a high-velocity case, this error can trigger FPE within a few iterations as the solver amplifies the inconsistency.

5. How to find which field diverged first

Enable field writing at every time step temporarily (writeInterval 1) and run for 5–10 iterations. Open the output in ParaView and look for the field with NaN values — that's where the divergence started.

Alternatively, add writeFormat ascii temporarily and look for nan or inf in the field files at the last written time step.

Using foamLog to extract residual history

The foamLog utility (available in all OpenFOAM distributions) parses the solver log and extracts residual histories into separate files that can be plotted:

foamLog log.pimpleFoam

This creates a logs/ directory with files like Ux_0, p_0, k_0, each containing the time and initial residual for that variable. Plot them with gnuplot or Python. The variable whose residual shows a spike immediately before the crash is the primary instability source.

6. Wrong turbulence model for the wall mesh resolution

A mismatch between the wall mesh resolution (y+) and the turbulence model wall treatment is a less obvious but frequent cause of FPE in RANS simulations. If y+ at the first wall cell is 0.5 (low-Re resolution) but the wall boundary uses high-Re wall functions (nutkWallFunction), the wall function extrapolation produces unphysical turbulent viscosity values that can overflow.

Conversely, if y+ = 50 (wall function resolution) but the model uses low-Re treatment, the turbulence equations attempt to resolve a boundary layer that is coarser than required, leading to instability in the viscous sublayer region.

Always verify wall treatment consistency:

// For y+ ≈ 1 (low-Re treatment):
wall
{
    type    kLowReWallFunction;    // 0/k
    value   uniform 0;
}
wall
{
    type    epsilonLowReWallFunction;  // 0/epsilon
    value   uniform 0;
}
wall
{
    type    nutLowReWallFunction;  // 0/nut
    value   uniform 0;
}

// For y+ 30-300 (high-Re wall functions):
wall
{
    type    kqRWallFunction;       // 0/k
    value   uniform 0;
}
wall
{
    type    epsilonWallFunction;   // 0/epsilon
    value   uniform 0;
}
wall
{
    type    nutkWallFunction;      // 0/nut
    value   uniform 0;
}