When the Mesh Skews, the Laplacian Leaks — Non-Orthogonal Diffusion Correction

A diffusion solver that is fine on cubes starts lying on skewed triangles. Same equation. Same code. The temperature field blurs and convergence slows. Not the boundary data. Not the material. The way the face gradient is measured.

Pushing the non-orthogonality slider from 20° to 65° stretches the over-relaxed cyan vector past the face area, and it stretches the minimum-correction pink vector even more.

When cyan grows shorter than pink, the ρ readout under the canvas crosses 1. At 45° only minimum correction should go red first.

The face gradient is the troublemaker#

In the finite-volume method, diffusion is a sum of face fluxes. The flux through one face $f$ is written as follows.

$\Gamma_f$ is the face diffusivity. $\nabla\phi_f$ is the face gradient. $\mathbf{S}_f$ is the outward area vector, area times normal.

The hard part is turning $\nabla\phi_f \cdot \mathbf{S}_f$ into cell values. On an orthogonal mesh the line $\mathbf{d}$ between cell centers lies along the normal. $(\phi_N-\phi_P)/|\mathbf{d}|$ is the normal gradient.

On an unstructured mesh that assumption breaks. $\mathbf{d}$ and $\mathbf{S}_f$ misalign. The angle is the non-orthogonality $\theta$. $(\phi_N-\phi_P)/|\mathbf{d}|$ is now the gradient along $\mathbf{d}$, not along the normal. Ignore the gap and diffusion leaks. The field smears.

Split S_f into two pieces#

Split the area vector into two parts. One along $\mathbf{d}$, one leftover.

$\mathbf{E}_f$ is the part along $\mathbf{d}$. $\mathbf{T}_f$ is the tangential remainder. The flux then splits.

The first term uses only $\phi_N-\phi_P$. It enters the matrix implicitly. Two cell centers, nothing else. The second term needs the face-average gradient $\overline{\nabla\phi}_f$, an average of neighbor gradients. That is a poor matrix coefficient, so it goes to the right-hand side as an explicit correction. Rebuilding that gradient by least squares is a separate job, written up in weights that keep a gradient alive on an unstructured mesh. When a deferred loop wobbles, the culprit is sometimes that gradient, not ρ.

If $\theta=0$, then $\mathbf{T}_f=0$. The old central difference comes back.

Minimum, orthogonal, over-relaxed: three splits#

The length of $\mathbf{E}_f$ along $\mathbf{d}$ is not unique. With $\hat{\mathbf{d}}=\mathbf{d}/|\mathbf{d}|$ and $\cos\theta=\hat{\mathbf{d}}\cdot\mathbf{S}_f/|\mathbf{S}_f|$, there are three choices.

• Minimum: $|\mathbf{E}_f|=|\mathbf{S}_f|\cos\theta$. Shortest $\mathbf{T}_f$.
• Orthogonal: $|\mathbf{E}_f|=|\mathbf{S}_f|$. The implicit magnitude stays at the face area.
• Over-relaxed: $|\mathbf{E}_f|=|\mathbf{S}_f|/\cos\theta$. The implicit part grows as the mesh skews.

On the figure above, θ at 60° makes the over-relaxed $\mathbf{E}_f$ longer than the area vector. The diagonal coefficient got bigger. Minimum correction does the opposite: $\mathbf{E}_f$ shrinks and $\mathbf{T}_f$ stretches.

Code: a deferred-correction Laplace solver#

The explicit term uses the gradient from the previous pass. One solve is not enough. The correction is refreshed over several sweeps. That is deferred correction.

import numpy as np
 
def decompose_face(S_f, d, scheme="over"):
    """Area vector S_f = E_f (implicit, along d) + T_f (explicit correction)."""
    d_hat = d / np.linalg.norm(d)
    Smag = np.linalg.norm(S_f)
    dS = d_hat @ S_f                      # = |S_f| cos(non-orthogonality)
    if scheme == "min":                   # minimum correction
        E = dS * d_hat
    elif scheme == "orth":                # orthogonal correction
        E = Smag * d_hat
    else:                                 # over-relaxed
        E = (Smag**2 / dS) * d_hat
    return E, S_f - E                     # E_f, T_f
 
def sweep_residual(S_f, d, scheme, sweeps=25):
    """E_f implicit, T_f lagged. Residual contracts at rho = |T|/|E|."""
    E, T = decompose_face(S_f, d, scheme)
    rho = np.linalg.norm(T) / np.linalg.norm(E)
    r, hist = 1.0, [1.0]
    for _ in range(sweeps):
        r *= rho                          # one deferred-correction sweep = residual x rho
        hist.append(r)
    return rho, hist
 
S_f = np.array([1.0, 0.0])                # face normal (area vector, |S|=1)
for ang in (20, 45, 65):
    d = np.array([np.cos(np.radians(ang)), np.sin(np.radians(ang))])
    for s in ("min", "orth", "over"):
        rho, hist = sweep_residual(S_f, d, s)
        tag = "converges" if rho < 1 else "diverges"
        print(f"theta={ang:2d}  {s:5s}  rho={rho:.3f}  {tag}  (after 25: {hist[-1]:.1e})")

The run says the three schemes pull apart as the angle grows. At 20° all three are mild. At 65° minimum correction diverges with $\rho=\tan 65°\approx 2.14$, while over-relaxed still contracts with $\rho=\sin 65°\approx 0.91$. Same mesh, same equation. The split decides convergence.

The three ratios close on the angle alone. With area 1, minimum correction has length of E equal to cosθ and length of T equal to sinθ, so ρ is tanθ. Tan of 45° is 1, and the deferred loop multiplies the residual by 1. Twenty-five sweeps change nothing. Tan of 20° is about 0.364. Tan of 65° is about 2.14.

Orthogonal correction has length of E equal to 1 and length of T equal to the square root of 2 minus 2cosθ. ρ hits 1 when cosθ is one half, at 60°. At 65° it is already about 1.07. Calling it unsafe only at 70° is late.

Over-relaxed ρ is sinθ, about 0.906 at 65°. The snippet multiplies the residual by that ratio 25 times. Raised to the 25th, 0.906 is about 0.085, so about 8 percent of the explicit error remains. At 20°, where sin is about 0.342, the same power falls near $10^{-12}$. Two or three correctors at 65° are 0.906 cubed, about 0.74. Three quarters of the error stays in that step. The 45° wall is exact, not an empirical threshold.

Why is ρ the length of T over the length of E? The implicit coefficient scales with E. Each sweep feeds back an error that scales with T. The residual shrinks or grows by that ratio. The loop closes only when the ratio is below 1.

What happens as non-orthogonality grows#

The three residuals run on one canvas. The angle slider sets ρ, and each sweep drops the curve by that ratio.

Below 30° the three curves fall together. At 45° only minimum correction goes flat. Past 60° orthogonal correction also crosses ρ of 1, and only over-relaxed correction, sitting on sinθ, remains.

There is a cost. Over-relaxation enlarges the implicit part and the explicit correction together. Reaching a tight residual takes more sweeps. nNonOrthogonalCorrectors is that sweep count. A clean mesh wants 0 or 1. A badly skewed mesh is pushed to 2 or 3. At 65° those two or three sweeps still leave about three quarters of the error. Adding more makes the time step expensive. Cutting the angle is almost always cheaper.

Why minimum correction is already dead at 45°#

At 45° minimum correction stops reducing the residual, because ρ is exactly 1. Orthogonal correction hits the same wall at 60° and is already past it at 65°. Over-relaxed correction stays under 1 until the angle is nearly 90°, because ρ is sinθ. The error left inside one step is the price.

Add more sweeps, or rebuild the mesh before the angle passes 60°? The order in which triangles cover the points sets that angle. Two ways of covering a point cloud with triangles asks where that order leaves non-orthogonality behind. Is a longer corrector loop the cheap fix, or is it cheaper to push the front again?