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

Same six neighbors, same pressures. One arrow sits on $\nabla p$. The other leans well off it. The gap is how skewed the cell is, and how the face pressure was built.

Push the mesh-skew slider from 0 to 0.9 and the orange Green–Gauss arrow leaves the true value, while the purple least-squares arrow stays shorter.

Near skew 0 both stay in single digits. Past 0.5 only the orange one climbs into double digits. The table sits after the code.

The day pressure blew up in the OpenFOAM log#

A non-staggered unstructured pressure-based solver calls the cell-centered $\nabla p$ at least three times inside Rhie–Chow pressure-weighted interpolation. The momentum pressure term, the face mid-velocity, and the gradient of the pressure correction. If one of them leans, momentum and continuity part and the residual stalls a decade up. In two-phase flow the phase-change rate is a function of local pressure, so a 0.5% bias in $\nabla p$ smears the interface.

The trouble is the face. The cell gradient falls out of a surface integral, and the recipe for the face pressure $p_f$ splits the road.

Green–Gauss — a single surface integral#

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

$p_f$ is the face pressure, $\hat{\mathbf{n}}_f$ the outward normal, $A_f$ the area.

The recipe for $p_f$ is the Green–Gauss variant.

• Simple average: $p_f=\tfrac12(p_C+p_N)$. Clean on an orthogonal uniform mesh. On an unstructured one the face centroid is not the midpoint of the two cell centers.
• Inverse distance: $p_f=(w_C p_C+w_N p_N)/(w_C+w_N)$, $w=1/d$. The cell closer to the face gets the larger weight.
• Frink: node pressures come from neighboring cells with pseudo-Laplacian weights, then the two nodes of the face are averaged again. Conservation and second order, both aimed at.

When the cell skews, the face centroid $\mathbf{x}_f$ leaves the line of centers. No average for $p_f$ recovers the distribution on the face. First order is often not even guaranteed. If the volume $\Omega_C$ in the denominator is a fake area, conservation is only formal. Closing the cell area with a diagonal cross product is where that area has to agree with $\partial x/\partial x=1$.

Same point, two algorithms, different arrows#

The ring in the figure is one center cell and six neighbors. At skew 0 both arrows sit on the cyan truth. Past 0.5 the orange one leaves first. Raise the wavenumber $k$ and the true $\nabla p$ gets steep relative to the cell. Both methods lose accuracy, and the purple lead narrows.

Least-Squares — fit the point cloud directly#

Least squares never visits the face. The neighbors are assumed linear, and the over-determined system is solved.

Each row leaves a residual $\varepsilon_N$. The weight $w_N$ is usually $1/|\mathbf{x}_N-\mathbf{x}_C|^2$. The $\nabla p$ that minimizes $\sum w_N\varepsilon_N^2$ is the normal equation.

In 2D that normal equation is 2×2.

$\Delta x_N=x_N-x_C$, and $\Delta y_N$, $\Delta p_N$ are the same kind of difference.

If the neighbor centers are not collinear the normal matrix is nonsingular, and a skewed face still keeps first order. That is why the method is dull to topology.

The weaknesses are real. A boundary corner with a single neighbor makes the normal matrix singular. The stencil has to reach vertex neighbors, and on an unstructured mesh that index is heavy. There is also no conservation. Green–Gauss is a sum of face fluxes, so cells telescope toward zero. Least squares only minimizes a residual.

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%}")

The print, with the caveat that it moves with the face shape, is this.

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%

The 1.6% least-squares error at skew 0 is not skew. The code sets $k=1.2$ and a ring radius $R=1.2$, so $kR=1.44$ rad. The field is $\sin(kx)\cos(ky)$, and at the origin the true gradient is $(k,0)=(1.2,0)$. A linear fit cannot hold that curvature. At skew 0 the normal matrix keeps only the angles. With $w=1/|\mathbf{d}|^2$, the sum of $w\Delta x^2$ is the sum of $\cos^2$. Six samples at 60 degrees sum to 3, the cross term is 0, and $M=3I$. The condition number is 1. The 1.6% is not a sick matrix.

The ratios in the table are 2.6, 3.4, 3.2, 3.4. That is $4.1/1.6$, $7.8/2.3$, $18.4/5.7$, $31.2/9.1$. From skew 0 to 0.9, Green–Gauss grows by $31.2/4.1\approx 7.6$ and least squares by $9.1/1.6\approx 5.7$. The Green–Gauss increments are $+3.7$, $+10.6$, $+12.8$. Not a straight line. The later steps are steeper.

At skew 0.9 the code's radius is $R(1\pm 0.36)$ and the angular kick has amplitude $0.54$ rad. Those are $0.9\times 0.4$ and $0.9\times 0.6$. The face length $0.55$ is hard-coded. When all six faces share that length it cancels between the numerator and the denominator. Replacing 0.55 with another constant leaves the printed 4.1% and 31% where they are. Those percentages belong to the ring angles and to the sine field.

At skew 0 the neighbor distance is $R=1.2$, so the face midpoint sits at half of that, 0.6. After the length drops out, each face contributes $0.5\times 0.6$ and six faces sum to 1.8. The volume the code accumulates, 0.99, multiplies by 0.55 again and never enters the quotient. A regular hexagon with neighbor spacing $R$ has area $(\sqrt{3}/2)R^2$. With $R^2=1.44$ that area is about 1.25. The canceled denominator 1.8 is 44% larger than that area. If the numerator matched the true surface integral, the gradient would fall to $1/1.44$ of the truth and the relative error would be $1-1/1.44\approx 0.31$. The printed Green–Gauss error is 4.1%, far below that 0.31. The face-pressure sum is large in the same proportion as the denominator, so the two geometry errors cancel most of the way and 4.1% is what remains. The least-squares 1.6% does not use this area.

The bend is readable from $k$ and $R$. $kR=1.44$, and the Taylor cubic $(kR)^3/6$ is about 0.50, roughly a third of the linear term $kR$. Samples on the ring are already bent by that much. A leftover gradient error of 1.6% means most of that bend leaves the gradient on a symmetric stencil. Raising skew breaks the cancellation. The factor 7.6 from 4.1% to 31.2% is not new curvature. It is the geometry error that had been canceling, now showing through the angle kick and the radius wobble.

Picking one in practice#

Pressure gradients are almost always least squares. Rhie–Chow has to keep the acoustic mode clean.

Velocity gradients, the viscous term, split. If conservation is required, a Green–Gauss variant. If accuracy comes first, least squares. Fluent ships both.

A boundary cell that cannot solve least squares alone gets a Green–Gauss patch or a ghost. OpenFOAM leastSquares also falls back when face neighbors run out.

In 3D the normal matrix is 3×3, and a rough cell drives the condition number up fast. Watching that number is the safety check.

Pitfalls when coding it up#

If $\det M$ sits on zero, least squares means nothing. Switch to an SVD or a pseudoinverse.

The weight $1/d^2$ boosts the near cells and stays dull to skew. $1/d$ or a uniform weight is smoother and less accurate.

Drop only the boundary to Green–Gauss and keep least squares inside, and the whole system falls to first order. Push least squares through the boundary, or use Green–Gauss everywhere.

Least squares can own the pressure gradient while the face pressure is still interpolated on its own. Forget that and mass imbalance stacks up.

Cut least squares at the boundary and the 9% falls with it#

The skew 0.9 row makes the choice a number. Pick conservation and eat the 31% of Green–Gauss. Pick accuracy and keep the 9% of least squares. Cut that least squares off at the boundary, though, and the single digit the interior was holding drops in order with it. The regular hexagon whose condition number was 1 does not exist at a boundary corner.

A correct gradient arrow still loses if the face velocity remembers the under-relaxation factor. When a log diverges, the open question is whether the arrow was wrong or the face remembered the relaxation. Where Rhie–Chow remembers under-relaxation is that question.