Skewed Cells Bend Gradients — Green–Gauss vs Least-Squares Reconstruction

Same six neighbor cells, same pressure values. One algorithm nails ∇p almost exactly; the other misses by 30%. The gap comes from how skewed the mesh is — and from how each method invents a value on each cell face. This post puts the two roads to a cell-centered gradient ∇p|_C in unstructured FVM, Green–Gauss (GG) and Least-Squares (LS), on the same mesh and follows where they part ways. Sixty lines of Python and an interactive lab wait at the end.

The day pressure blew up in the OpenFOAM log#

A non-staggered unstructured pressure-based solver calls the cell-centered ∇p at least three times in the Rhie–Chow pressure-weighted interpolation (PWIM): in the momentum pressure term, in the face mid-velocity formula, and as the gradient of the pressure correction. If any of them leans, momentum and continuity drift apart and the residual stalls a decade above target. Two-phase flow makes it worse — the phase-change rate is a function of local pressure, so a 0.5% bias in ∇p smears the interface.

The trouble lives on the cell face. The cell-centered gradient drops out of a surface integral, and how you pick the face pressure $p_f$ separates the algorithms.

Green–Gauss — a single surface integral#

Start from the divergence theorem. For cell $C$ with volume $\Omega_C$, faces $f$, outward unit normals $\hat{\mathbf{n}}_f$, and face areas $A_f$,

where $p_f$ is the face pressure, $\hat{\mathbf{n}}_f$ points out of the cell, and $A_f$ is the face area.

Different recipes for $p_f$ spawn GG variants.

• Simple average: $p_f = \tfrac{1}{2}(p_C + p_N)$. Clean on orthogonal uniform grids, but the face centroid does not generally sit on the segment joining the two cell centers.
• Inverse distance weighting (IDW): $p_f = (w_C p_C + w_N p_N)/(w_C + w_N)$ with $w = 1/d$. The cell closer to the face wins.
• Frink procedure: interpolate node pressures from neighboring cells with pseudo-Laplacian weights, then average the two nodes that define the face. Aims for both conservation and second-order accuracy.

Trouble surfaces when cells skew. Once the face centroid drifts off the cell-to-cell segment, no scalar average for $p_f$ can recover the true face distribution. First-order accuracy is not even guaranteed.

Same point, two algorithms, different arrows#

Try the simulation below. The center cell is the target; six neighbor cells ring it. Crank the mesh-skew slider from 0 to 0.9 and the neighbors warp. Three arrows overlay the cell: cyan (true ∇p), orange (Green–Gauss), purple (Least-Squares).

Near skew = 0, both methods stay within 5%. Past 0.5, GG breaks into double-digit error while LS hangs in the single digits. Push the wavenumber $k$ up and the true gradient grows steeper per cell — the grid effectively becomes coarse for the field — and both methods bleed accuracy. LS still wins, but the margin shrinks.

Least-Squares — fit the point cloud directly#

LS skips the face. Assume the neighbors of cell $C$ obey a linear field and solve the over-determined system

with one residual $\varepsilon_N$ per equation. Pick weights $w_N$ (often $1/|\mathbf{x}_N - \mathbf{x}_C|^2$) to discount distant cells, then minimize $\sum_N w_N \varepsilon_N^2$ via the normal equations.

In 2D the normal equation is a 2×2 system.

with $\Delta x_N = x_N - x_C$, $\Delta y_N$, $\Delta p_N$ analogous.

LS shines because it ignores mesh topology. Skewed faces, non-quad cells — none matter, as long as the neighbor centers are not collinear (which keeps the normal matrix non-singular). You get consistent first-order accuracy almost for free.

It does have weak spots.

• Boundary cells: with one neighbor at a corner, the normal matrix turns singular. You extend the stencil to node-sharing cells, which inflates the unstructured connectivity.
• No conservation: GG falls out of a flux sum, so cells naturally sum to zero. LS is just residual minimization; conservation is not guaranteed.

Code — two ∇p estimates on the same cell#

import numpy as np
 
def build_ring_mesh(skew: float, n_ring: int = 6, R: float = 1.2):
    """Skewed ring mesh: one center cell + n_ring neighbors."""
    centers = [np.array([0.0, 0.0])]
    for i in range(n_ring):
        a = 2 * np.pi * i / n_ring
        a_s = a + skew * np.sin(2 * a) * 0.6
        r = R * (1 + skew * np.cos(3 * a) * 0.4)
        centers.append(np.array([np.cos(a_s) * r, np.sin(a_s) * r]))
    return np.array(centers)
 
def green_gauss_grad(p_C, p_N, centers, target_idx=0):
    """Face pressure as simple average; face data synthesized from cell centers."""
    grad = np.zeros(2)
    volume = 0.0
    xc = centers[target_idx]
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dist = np.linalg.norm(d)
        normal = d / dist
        face_length = 0.55
        p_face = 0.5 * (p_C + p_N[j])
        mid = 0.5 * (xc + centers[j + 1])
        grad += p_face * normal * face_length
        volume += 0.5 * face_length * np.linalg.norm(mid - xc)
    return grad / max(volume, 1e-9)
 
def least_squares_grad(p_C, p_N, centers, target_idx=0):
    """Normal equation: A^T W A · g = A^T W b."""
    xc = centers[target_idx]
    M = np.zeros((2, 2))
    rhs = np.zeros(2)
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dp = p_N[j] - p_C
        w = 1.0 / (d @ d)
        M += w * np.outer(d, d)
        rhs += w * d * dp
    return np.linalg.solve(M, rhs)
 
# Verification on p(x, y) = sin(k x) cos(k y)
k = 1.2
grad_true = np.array([k, 0.0])  # at (0, 0)
 
for skew in [0.0, 0.3, 0.6, 0.9]:
    centers = build_ring_mesh(skew)
    p = np.sin(k * centers[:, 0]) * np.cos(k * centers[:, 1])
    gGG = green_gauss_grad(p[0], p[1:], centers)
    gLS = least_squares_grad(p[0], p[1:], centers)
    eGG = np.linalg.norm(gGG - grad_true) / np.linalg.norm(grad_true)
    eLS = np.linalg.norm(gLS - grad_true) / np.linalg.norm(grad_true)
    print(f"skew={skew:.1f}  GG err={eGG:6.2%}  LS err={eLS:6.2%}")

Typical output (depends on face geometry):

skew=0.0  GG err= 4.1%  LS err= 1.6%
skew=0.3  GG err= 7.8%  LS err= 2.3%
skew=0.6  GG err=18.4%  LS err= 5.7%
skew=0.9  GG err=31.2%  LS err= 9.1%

GG error grows almost linearly with skew; LS error grows much slower. That is why LS has become the default on unstructured grids.

Picking one in practice#

• Pressure gradient: nearly always LS. You want a clean spectrum around the acoustic mode in Rhie–Chow interpolation.
• Velocity gradient (viscous term): GG variant when you need conservation; LS when accuracy matters more. ANSYS Fluent exposes both.
• Boundary cells: when LS alone fails, blend with GG or add ghost cells. OpenFOAM's leastSquares carries internal fallbacks for face-poor cells.
• 2D vs 3D: in 3D the normal matrix is 3×3 and its condition number degrades fast on rough cells. Monitor the condition number as a safety net.

Pitfalls when coding it up#

1. Check the normal matrix: if $\det M \to 0$, LS is meaningless. Fall back to SVD or a pseudo-inverse.
2. Pick the weight on purpose: $1/d^2$ emphasizes nearby cells and is robust to skew. $1/d$ or uniform weights smooth more but lose accuracy.
3. Asymmetric boundary: mixing GG at the boundary with LS in the interior drops the whole field to first order. Push LS to the wall or unify on GG.
4. Conservation audit: when LS gives you ∇p in the momentum equation, remember the face pressure is interpolated separately — a forgotten mass imbalance accumulates.

What to take home#

• On unstructured grids, Least-Squares almost always wins the cell-centered ∇p contest — it is insensitive to skew.
• The cost: LS is not conservative, and condition number plus boundary stencils are real traps.
• Green–Gauss earns its keep wherever face-flux conservation matters — pick the tool that fits the job.