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 CFDpilot · Updated April 2026

What the error actually means

OpenFOAM compiles with floating point exception trapping enabled by default (-ftrapping-math). When any arithmetic produces NaN (division by zero, 0/0, sqrt of negative) or Inf (overflow), the OS sends a SIGFPE signal 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.

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;
}