Pushing the cell Péclet number to 2 drove the solution to -0.33 — where the energy to minimize disappeared

Four methods gave the same answer on one bar#

Chapter 1 of a finite element textbook solves the same tapered-bar problem four times: direct stiffness, minimum total potential energy, weighted residuals, and Galerkin. All four produce the identical 5×5 stiffness matrix. The lecture notes put it plainly — whichever method you pick, you land on essentially the same accurate answer.

That sentence carries a condition, and the condition hides inside the bar problem where nobody trips over it. Flow equations do not satisfy it. The measurement is which guarantee disappears when the condition breaks, and which part stabilization can hand back. Raise advection a from zero and the pink dot's lap work leaves 0.000 for exactly $2\pi a$.

The amber ball that used to slide inward without cutting the grey level curves starts spiraling up them once a grows. Push diffusion s to the stop and the lap total stays put.

A minimum needs a symmetric matrix#

The minimum total potential energy formulation differentiates the total potential with respect to the nodal unknowns and sets the result to zero. In discrete form:

Here $\boldsymbol{\phi}$ is the vector of nodal unknowns, $\mathbf{K}$ the stiffness matrix, $\mathbf{f}$ the load vector. Now differentiate with respect to component $i$:

Note what appeared: $K_{ij} + K_{ji}$. For this to reduce to $\mathbf{K}\boldsymbol{\phi} = \mathbf{f}$, you need $K_{ij} = K_{ji}$. Without symmetry, minimization actually solves the symmetrized matrix $\tfrac{1}{2}(\mathbf{K}+\mathbf{K}^{T})$, which is a different equation from the one you wrote down.

There is a deeper version of the same statement. For the residual field $\mathbf{r}(\boldsymbol{\phi}) = \mathbf{f} - \mathbf{K}\boldsymbol{\phi}$ to be the gradient of some scalar function, its Jacobian must be symmetric. If it is not, no potential function exists at all. The fastest way to test existence is to walk a closed loop: a gradient field always returns zero work after a full lap. The climb over the level curves, once a leaves zero, is that lap failing to close at zero.

Advection breaks it by exactly $u$#

Write the weak form of the one-dimensional advection-diffusion equation $u\,\phi' = \epsilon\,\phi''$ and the two terms show their different characters:

Here $w$ is the test function, $\epsilon$ the diffusivity, $u$ the advection velocity. The first term is unchanged when you swap $w$ and $\phi$ — it is symmetric. The advection term flips sign under integration by parts. For test functions vanishing at the boundaries, $\int w\,u\,\phi'\,dx = -\int \phi\,u\,w'\,dx$, so advection is a purely skew-symmetric contribution.

Assemble with linear elements and that structure survives in the coefficients. With element size $h$, an interior row reads:

Diffusion deposits the same value on both neighbors; advection adds $+u/2$ on one side and $-u/2$ on the other. So the departure from symmetry is exactly one number:

Refine the mesh as much as you like and this stays $u$, because $h$ never entered it. As long as there is flow, the minimum total potential energy principle is not coming back. Folding a stiffness into a 6×6 Voigt matrix only worked while the matrix stayed symmetric. When advection makes a second-order scheme ring near a step, the fix lives on the TVD advection side, with a flux limiter.

Three schemes on the same mesh in Python#

Ten elements, $u = 1$, boundary conditions $\phi(0)=0$ and $\phi(1)=1$. A single artificial-diffusion coefficient $\beta$ generates all three schemes through $\epsilon_{\text{eff}} = \epsilon + \beta\,u\,h/2$: $\beta = 0$ is pure Galerkin, $\beta = 1$ is full upwind, and $\beta = \coth(Pe_h) - 1/Pe_h$ is SUPG.

import numpy as np
 
def assemble_ad(n, vel, eps, beta):
    """1D advection-diffusion on linear elements. beta = artificial diffusion."""
    h = 1.0 / n
    eps_eff = eps + beta * vel * h / 2.0
    kd = (eps_eff / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])   # diffusion: symmetric
    ka = (vel / 2.0) * np.array([[-1.0, 1.0], [-1.0, 1.0]])     # advection: skew
    K = np.zeros((n + 1, n + 1))
    for e in range(n):
        K[e:e + 2, e:e + 2] += kd + ka
    return K
 
def solve_bvp(K):
    n = K.shape[0] - 1
    A, b = K.copy(), np.zeros(n + 1)
    A[0, :], A[0, 0], b[0] = 0.0, 1.0, 0.0
    A[n, :], A[n, n], b[n] = 0.0, 1.0, 1.0
    return np.linalg.solve(A, b)
 
def exact_ad(x, pe):
    return (np.exp(pe * (x - 1.0)) - np.exp(-pe)) / (1.0 - np.exp(-pe))
 
def skew_ratio(K):
    return np.linalg.norm(K - K.T) / np.linalg.norm(K + K.T)
 
def loop_work(K, m=20000):
    """Work of the residual field -K.phi once around the unit circle in DOF space."""
    t = np.linspace(0.0, 2.0 * np.pi, m, endpoint=False)
    path = np.stack([np.cos(t), np.sin(t)])          # position
    tang = np.stack([-np.sin(t), np.cos(t)])         # dl / dt
    return float(np.sum(np.sum(-(K @ path) * tang, axis=0)) * (2.0 * np.pi / m))
 
n, vel = 10, 1.0
x = np.linspace(0.0, 1.0, n + 1)
print(f"{'Pe_h':>5} {'scheme':>9} {'max|err|%':>10} {'min phi':>9} {'skew/sym':>9} {'loop W':>8}")
for pe_h in [0.5, 1.0, 2.0, 5.0]:
    eps = vel / (2.0 * n * pe_h)
    ex = exact_ad(x, vel / eps)
    for name, beta in [("Galerkin", 0.0), ("upwind", 1.0),
                       ("SUPG", 1.0 / np.tanh(pe_h) - 1.0 / pe_h)]:
        K = assemble_ad(n, vel, eps, beta)
        phi = solve_bvp(K)
        print(f"{pe_h:5.1f} {name:>9} {100 * np.max(np.abs(phi - ex)):10.2f} "
              f"{phi.min():9.4f} {skew_ratio(K):9.4f} {loop_work(K[4:6, 4:6]):8.4f}")
 
K0 = assemble_ad(n, 0.0, 0.1, 0.0)
print(f"\nvel = 0 : skew/sym = {skew_ratio(K0):.2e},  loop W = {loop_work(K0[4:6, 4:6]):.2e}")
print(f"pi * u  = {np.pi * vel:.4f}")

 Pe_h    scheme  max|err|%   min phi  skew/sym   loop W
  0.5  Galerkin       3.45    0.0000    0.2924   3.1416
  0.5    upwind      13.17    0.0000    0.1954   3.1416
  0.5      SUPG       0.00   -0.0000    0.2704   3.1416
  1.0  Galerkin      13.53    0.0000    0.5774   3.1416
  1.0    upwind      19.80    0.0000    0.2924   3.1416
  1.0      SUPG       0.00   -0.0000    0.4428   3.1416
  2.0  Galerkin      35.17   -0.3334    1.1010   3.1416
  2.0    upwind      18.17   -0.0000    0.3885   3.1416
  2.0      SUPG       0.00    0.0000    0.5572   3.1416
  5.0  Galerkin      69.61   -0.6961    2.1517   3.1416
  5.0    upwind       9.09    0.0000    0.4836   3.1416
  5.0      SUPG       0.00    0.0000    0.5773   3.1416
 
vel = 0 : skew/sym = 0.00e+00,  loop W = -4.29e-16
pi * u  = 3.1416

Push the Pe_h slider past 1 and the a_E bar crosses red while the nodes sink below zero.

Leave beta at zero and set Pe_h to 2. The nodal minimum lands near the table's -0.3334. Hit SUPG and the error sticks at 0 %, while the pink skew line stays on 3.1416.

The boundary values are pinned at 0 and 1, yet the Galerkin solution at $Pe_h = 2$ dips to $-0.3334$. At $Pe_h = 5$ it reaches $-0.6961$. Neither value can exist physically.

One coefficient flips sign at $Pe_h = 1$#

Define the cell Péclet number as $Pe_h = u h / (2\epsilon)$ and the right-hand coefficient of that row rewrites itself:

Once $Pe_h > 1$ the off-diagonal entry turns positive. At that instant the matrix stops being an M-matrix, and the discrete maximum principle — the guarantee that interior values stay inside the range set by the boundaries — leaves with it. That is why min phi in the table holds at zero through $Pe_h = 1.0$ and goes negative at 2.0. The moment a_E crossed red above is this sign change.

Stabilization does not restore symmetry#

A common misreading shows up here. When upwinding or SUPG is said to "restore stability", that does not mean symmetry and the minimization principle come back. Look at the loop W column: three schemes across four values of $Pe_h$, and all twelve rows read 3.1416. That number is $\pi u$, with neither $\epsilon$ nor $\beta$ anywhere in it.

The reason is simple. All $\beta$ does is inflate $\epsilon_{\text{eff}}$, which lives in the symmetric part of the matrix. The skew part, $u$, is untouched. The skew/sym ratio dropping from 2.1517 to 0.4836 at $Pe_h = 5$ is not the skew part shrinking; it is the symmetric denominator growing.

At $Pe_h = 0.5$ Galerkin prints max|err|% 3.45, min phi 0.0000, and skew/sym 0.2924. At $Pe_h = 1.0$ the error rises to 13.53 while the minimum is still 0.0000. From $K_{i,i+1} = (\epsilon/h)(Pe_h - 1)$, the right-hand coefficient is zero at $Pe_h = 1$ and turns positive the moment the cell Péclet crosses 1. At $Pe_h = 2$ Galerkin's minimum is -0.3334, the error is 35.17 %, and skew/sym is 1.1010. On the same row upwind has minimum -0.0000, error 18.17 %, and skew/sym 0.3885. SUPG has error 0.00 and skew/sym 0.5572. All three loop W entries read 3.1416. With $u = 1$, $\pi u$ is 3.1416, so the lap work is the same whether $\beta$ is 0, 1, or $\coth(Pe_h) - 1/Pe_h$. At $Pe_h = 5$ Galerkin reaches minimum -0.6961, error 69.61 %, and skew/sym 2.1517. Upwind cuts the error to 9.09 % and skew/sym to 0.4836. What shrank is the symmetric denominator. The zero-velocity check prints skew/sym 0.00e+00 and loop W -4.29e-16. Ten elements and $u = 1$ give $h = 0.1$, and $Pe_h = uh/(2\epsilon)$ then puts $\epsilon = 0.025$ at $Pe_h = 2$. Inflating that $\epsilon$ with $\beta$ still leaves $K_{i,i+1} - K_{i+1,i} = u$ equal to 1.

So what stabilization actually buys is a weaker guarantee. You gave up the minimization principle — optimality in the energy norm — and in exchange you get the M-matrix property and the discrete maximum principle. You pay in accuracy. At $Pe_h = 2$, upwinding killed the oscillation but carries 18.17 % error. SUPG adds only as much $\beta$ as needed and stays exact at the nodes.

This goes to 0 as $Pe_h \to 0$ and to 1 as $Pe_h \to \infty$: turn stabilization off when diffusion dominates, go full upwind when advection does. Nodal exactness, though, is a privilege of the one-dimensional constant-coefficient problem. In two dimensions you need the original SUPG form, which adds artificial diffusion only along the streamline, and even then "exact" is gone.

The name finite volume uses for the same place#

The calculation above spoke finite element, but the conclusion is independent of the discretization. In finite volume, a central-differenced convection term produces exactly the same stencil and flips sign at exactly the same $Pe_h = 1$. That is what first-order upwinding answers, and the numerical diffusion it introduces equals the $u h / 2$ you get at $\beta = 1$.

Where the two worlds part company is conservation. Finite volume upwinding modifies the face flux, so the global balance survives intact. It is the same kind of split as conservative versus primitive forms disagreeing on shock speed. Finite element artificial diffusion adds a term to the stiffness matrix, so what it still conserves has to be checked on its own.

Choosing a different member of the weighted-residual family costs something else. Least squares minimizes the residual norm, so it always produces a symmetric positive definite matrix — the minimization principle returns. In exchange the condition number is squared, and on linear elements the second derivative vanishes inside each element, so the diffusion term drops out entirely. The basis and the integration rule quietly change the character of the formulation.

At $Pe_h = 2$, buy a bound before buying a finer mesh#

The textbook line about all methods agreeing holds only for self-adjoint operators. Diffusion, elasticity, and potential flow sit inside that region. Advection brings $K_{i,i+1} - K_{i+1,i} = u$, the symmetry breaks, and loop W stays at 3.1416 — $\pi u$ — no matter what $\epsilon$ does. The energy-norm optimum that CG promised is gone. A nonsymmetric matrix belongs with the Krylov side, from Jacobi through BiCGSTAB.

A negative min phi means the cell Péclet has passed 1. At $Pe_h = 2$ Galerkin reaches -0.3334 with 35.17 % error. Upwind kills the oscillation and leaves 18.17 % error. SUPG is exact at the nodes, yet skew/sym only falls from 1.1010 to 0.5572. The skew piece itself did not move. In three dimensions, taking $Pe_h$ from 5 down to 1 multiplies the cell count by 125. When that factor is unaffordable, what gets bought is boundedness, not accuracy.