When Gradients Go Wrong on Skewed Meshes — Green–Gauss vs. Weighted Least-Squares Reconstruction

The residuals leaned over in a strange way, then blew up. Digging through the logs, I found one cell near the wall where the temperature gradient pointed almost perpendicular to the real one. I looked at the mesh: that one cell was oddly flat. The values were fine. The method that builds the gradient had tripped over the shape of the mesh.

In the finite volume method (FVM), the cell gradient feeds the diffusion flux, second-order reconstruction, and the limiter. Cell-center values alone are not enough. You have to gather neighbor information to estimate $\nabla\phi_C$. Two estimates sit on the table: Green–Gauss and least-squares. The run below shows where distortion swings one of them, and what the distance weight pulls back.

Drop anisotropy to about 0.2 and the stencil flattens. Raise the weight exponent from 0 to 1 and the orange arrow swings back toward the green truth.

Raise curvature κ at the same time and the gap between the arrows opens. On a flat stencil, watch whether the weight restores the direction.

Green–Gauss: wrap the cell in faces and average#

Apply the divergence theorem directly to the cell. Turn the volume integral into a face integral.

$V_C$ is the cell volume, $\mathbf{n}_f$ is the outward face normal, $A_f$ is the face area, and $\phi_f$ is the value on the face.

You cannot know the face value $\phi_f$ directly. You interpolate it from cell-center values. The most common linear interpolation looks like this.

$g_f$ is a weight set by the distance ratio from the face to the neighbor center $N$. It is cheap and intuitive. On structured meshes it gives second-order accuracy.

Skewed meshes shake Green–Gauss#

The trouble hides in $g_f$. Linear interpolation assumes the face center and the two cell centers all lie on one straight line. On unstructured meshes that assumption almost always breaks.

When the face center strays off the line joining the two cell centers ($\text{skewness}$), the interpolated $\phi_f$ becomes biased. If cells are flat or neighbor sizes differ a lot, the bias grows. This error accumulates over iterations. In the end the gradient points the wrong way.

Adding a correction term eases it. But the correction needs a gradient again. That creates a loop. Iterative correction is expensive, and if the distortion is severe, convergence slows down.

Least-squares: fit a plane to the neighbors#

Change the idea. Skip the face entirely. Fit a single plane straight to the value distribution around cell $C$.

If neighbor $k$ sits at offset $\mathbf{d}_k = \mathbf{x}_k - \mathbf{x}_C$, the linear approximation is $\phi_k \approx \phi_C + \nabla\phi_C \cdot \mathbf{d}_k$. Minimize the squared sum of this error over all neighbors.

$w_k$ is a per-neighbor weight. Differentiate and set to zero, and you get the normal equations.

In two dimensions this is a $2\times2$ system. Written out, it looks like this.

$\Delta\phi = \phi_k - \phi_C$. The left-hand matrix depends only on neighbor placement. You can compute it once per cell and store it. The upside is large. There is no face-interpolation assumption. It is insensitive to distortion. With three or more neighbors it always has a solution.

The face interpolation breaks in the same place as Green–Gauss versus least squares on a skewed cell. Three neighbors are enough to solve the normal equations, but if those three sit on a line the transverse answer is empty.

Weight by distance#

The weight $w_k$ decides the accuracy. A common choice is the inverse distance.

With $n=0$, every neighbor counts equally. A far neighbor has a large value difference, and with it a large curvature error. That neighbor pulls the plane fit. Raise $n$ to 1 and the far neighbor loses its voice. Weight shifts to nearby neighbors, where the local linear assumption holds well. For inviscid problems, $n=1$ is the usual choice. It cuts the far neighbor's vote and trims the curvature error.

Python — exact on linear fields, weighted on curved ones#

We confirm the power of least-squares in two parts. First, reconstruct the gradient of the linear field $\phi = 1.3x - 0.7y + 2$ on a randomly scattered stencil. In theory the error should be zero. Next, on the curved surface $\phi = e^{0.7x}\cos(0.9y)$, we pit the unweighted case against $n=1$ weighting.

import numpy as np
 
def weighted_lsq_gradient(offsets, dphi, n):
    # offsets: (K,2) neighbor offsets d_k; dphi: (K,) diffs phi_k - phi_C
    w = 1.0 / np.maximum(np.linalg.norm(offsets, axis=1), 1e-12) ** n
    M = (offsets[:, :, None] * offsets[:, None, :] * w[:, None, None]).sum(axis=0)
    b = (offsets * (dphi * w)[:, None]).sum(axis=0)
    return np.linalg.solve(M, b)
 
xc = np.array([0.4, -0.2])
 
# (A) linear field — least-squares is exact for any stencil
lin = lambda x, y: 1.30 * x - 0.70 * y + 2.0
scatter = np.array([[0.45, 0.30], [-0.40, 0.12], [0.05, -0.48], [-0.33, -0.29],
                    [0.90, 0.15], [-0.11, 0.44], [0.28, -0.09]])   # irregular stencil
dlin = np.array([lin(*(xc + d)) - lin(*xc) for d in scatter])
gA = weighted_lsq_gradient(scatter, dlin, 1.0)
print(f"linear field   LSQ = ({gA[0]:+.6f}, {gA[1]:+.6f})   exact = (+1.300000, -0.700000)")
 
# (B) curved surface — distance weighting trims the reconstruction error
phi  = lambda x, y: np.exp(0.7*x) * np.cos(0.9*y)
grad = lambda x, y: np.array([0.7*np.exp(0.7*x)*np.cos(0.9*y),
                             -0.9*np.exp(0.7*x)*np.sin(0.9*y)])
near = np.array([[0.10, 0.05], [-0.08, 0.09], [0.06, -0.10], [-0.09, -0.06], [0.11, 0.10]])
offs = np.vstack([near, [0.90, 0.70]])          # one far-away neighbor
dphi = np.array([phi(*(xc + d)) - phi(*xc) for d in offs])
gt   = grad(*xc)
for n in (0.0, 1.0):
    g = weighted_lsq_gradient(offs, dphi, n)
    err = np.linalg.norm(g - gt)
    print(f"curved field   n={n:.0f}  LSQ = ({g[0]:+.3f}, {g[1]:+.3f})   |error| = {err:.4f}")
print(f"curved field   true = ({gt[0]:+.3f}, {gt[1]:+.3f})")

The output comes out like this.

linear field   LSQ = (+1.300000, -0.700000)   exact = (+1.300000, -0.700000)
curved field   n=0  LSQ = (+0.885, +0.200)   |error| = 0.0293
curved field   n=1  LSQ = (+0.891, +0.205)   |error| = 0.0222
curved field   true = (+0.911, +0.213)

On a linear field the stencil can be as twisted as you like and the gradient is still good to six decimals. The printout (+1.300000, -0.700000) is that value. A linear function matches the dot product of offset and gradient at every neighbor, so the irregular point (0.90, 0.15) adds nothing to the error. Distortion does not shove the answer on this field.

A curved field does. The true gradient is (+0.911, +0.213). At n=0 the reconstruction (+0.885, +0.200) is short by 0.026 in x and 0.013 in y. At n=1, (+0.891, +0.205) is short by 0.020 and 0.008. The norm drops from 0.0293 to 0.0222. The gap, 0.0071, is about 24 percent, the one-quarter cut. The far neighbor (0.90, 0.70) has components much larger than the 0.1 scale of the five near points. At n=0 that point votes with the same weight as the rest. At n=1 distance cuts its vote.

When the neighbors line up, do not trust the inverse#

Least-squares being robust to distortion does not make it a cure-all. The left-hand matrix $M = \sum w_k \mathbf{d}_k \mathbf{d}_k^{\mathsf T}$ can become ill-conditioned. When the neighbors line up along a single straight line, $M$ approaches singular. The gradient is well determined along the stencil direction, but almost free in the transverse direction. The condition number explodes.

Push Collinearity toward 0.9 and the neighbors bunch into one line. cond(M) jumps from the tens to infinity.

The long thin axis of the ellipse is the axis the gradient cannot pin down. Boundary-layer cells pressed to the wall look like that. The diagonal entries are sums of squared offsets, so a stencil stretched one way leaves the other diagonal empty. On a boundary cell the neighbors sit on the fluid side. Adding the neighbor across the wall, or a face center, fills a little of that empty axis and the condition number drops. If the mesh is fixed, $M$ does not care that the velocity field changed. Store the inverse, or the QR factorization, instead of rebuilding it every step. The exponent $n$ is a problem choice. When the viscous term is large, $n=2$ is sometimes the better one.

Moving least squares rebuilds that weight neighbor by neighbor. Green–Gauss is cheap and weak to skew. Weighted least squares handles skew, until the stencil collapses onto a line and the matrix goes bad first. Before picking a method, check whether that cell's neighbors actually span the plane.