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. This post measures what the condition is, and exactly which guarantee you lose when it breaks. Then it checks, with numbers, what stabilization schemes hand back — and what they can never hand back.

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.

Try it below on a toy matrix with two degrees of freedom.

Set advection a to zero and the pink dot completes its lap carrying 0.000 back. The grey ellipses are energy level curves, and the amber ball slides straight inward without cutting across them. Raise a and the lap work becomes exactly $2\pi a$, while the ball starts spiraling and climbing back over the level curves. Now push diffusion s all the way up. The lap total does not move at all.

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 second 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. That is a different situation from transforming the constitutive tensor in curvilinear coordinates, where symmetry was what let the stiffness fold into a 6×6 Voigt matrix.

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

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 simulation below marches the same system in time, so you can watch where the oscillation grows as the steady state assembles itself.

Leave beta at zero and push Pe_h past 1. The a_E bar crosses to the left and turns red, and on the next time steps the nodal values sink below zero. Press SUPG and the error drops to 0 %. But the pink line in the lower right — the skew part of the matrix — does not move under any of the three buttons.

Stabilization does not restore symmetry#

A common misreading is worth clearing 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.

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. The story of how conservative and primitive forms split on shock speed repeats here. Finite element artificial diffusion adds a term to the stiffness matrix, so you have to check separately what it still conserves.

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. As in the post on the quadrature floor for discontinuous Galerkin, this is a place where the basis and the integration rule quietly change the character of the formulation.

When you inherit a non-symmetric matrix#

The textbook line about all methods agreeing is valid over self-adjoint operators. Diffusion, elasticity, and potential flow live inside that region. The moment advection appears, you are outside it.

In practice, check three things in order. First, is the assembled matrix symmetric? If yes, CG-family solvers apply and optimality in the energy norm comes along. If not, you are on GMRES and that guarantee is gone. Second, does the cell Péclet number exceed 1? If it does, the oscillation is not a bug — it is the defined behavior of the scheme. Third, if stabilization is on, did it buy accuracy or boundedness? Almost always boundedness, with accuracy as the currency.

When the solution escapes the boundary range, refining the mesh is not a workaround but the direct fix, because shrinking $h$ shrinks $Pe_h$ with it. Just note that in three dimensions, bringing $Pe_h$ from 5 down to 1 multiplies the cell count by 125. Run that arithmetic before you pick a stabilization term.