RasmahRasmah

Stabilized Flow (PSPG & SUPG)

Overview

Incompressible-flow finite elements pair a velocity space with a pressure space, and the inf–sup (LBB) condition forces them to be different — e.g. the Taylor–Hood P2/P1 pair used in the Stokes tutorial. If instead you use the convenient equal-order P1/P1 pair, the discrete saddle-point system is unstable: the pressure develops spurious checkerboard oscillations.

Stabilization fixes this by adding mesh-dependent terms to the weak form, so equal-order elements become usable. Rasmah implements the two classical methods (Tezduyar 1991; Brooks & Hughes 1982; Hughes, Franca & Balestra 1986):

  • PSPG (pressure-stabilizing Petrov–Galerkin) — stabilizes Stokes and Navier–Stokes pressure.
  • SUPG (streamline-upwind Petrov–Galerkin) — stabilizes advection-diffusion at high Péclet number.

The stabilization parameters

pspg_tau and supg_tau compute the element stabilization parameter from the local mesh size and the problem coefficients. The stabilization parameter scales with the square of the local element size — small elements need almost no stabilization, large elements need a lot:

using Rasmah

# a 2D channel mesh for the Poiseuille test
function channel_grid(nx, ny)
    nv = (nx + 1) * (ny + 1)
    verts = Matrix{Float64}(undef, 3, nv)
    for j in 0:ny, i in 0:nx
        idx = i + j * (nx + 1) + 1
        verts[1, idx] = i / nx; verts[2, idx] = j / ny; verts[3, idx] = 0.0
    end
    faces = Matrix{Int}(undef, 3, 2 * nx * ny)
    k = 0
    for j in 0:(ny - 1), i in 0:(nx - 1)
        a = i + j * (nx + 1) + 1; b = a + 1; c = a + nx + 1; d = c + 1
        k += 1; faces[:, k] = [a, b, c]
        k += 1; faces[:, k] = [b, d, c]
    end
    TriangleMesh(verts, faces)
end

m = channel_grid(8, 8)

μ = 1.0
U = FESpace(m, ReferenceFE(2, 1); valuetype = VectorValue{2})
P = FESpace(m, ReferenceFE(2, 1))

println("PSPG tau (uniform): ", round.(pspg_tau(m, μ); digits = 6))
println("SUPG tau (β = 1, D = 1): ", round.(supg_tau(m, 1.0, 1.0); digits = 6))
PSPG tau (uniform): [0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604  …  0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604, 0.002604]
SUPG tau (β = 1, D = 1): [0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603  …  0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603, 0.002603]

Equal-order Stokes (PSPG) vs. the exact solution

Pressure-driven channel flow has the exact parabolic Poiseuille profile uₓ(y) = (G/2μ) y (1 − y). The stabilized equal-order P1/P1 discretization reproduces it even though the raw P1/P1 saddle system would be unstable:

G = 1.0
nU = ndofs(U); nP = ndofs(P)

# prescribe the exact parabolic velocity on every wall
ufixed = Int[]; ug = Float64[]
for i in 1:U.nnodes
    x = U.nodes[:, i]
    onb = x[1] < 1e-9 || x[1] > 1 - 1e-9 || x[2] < 1e-9 || x[2] > 1 - 1e-9
    if onb
        ux = (G / (2μ)) * x[2] * (1 - x[2])
        append!(ufixed, 2i - 1, 2i)
        append!(ug, ux, 0.0)
    end
end
fixed = vcat(ufixed, nU .+ [1])        # pin one pressure dof
g = vcat(ug, [0.0])

K = assemble_stokes_pspg(U, P, μ)      # stabilized equal-order P1/P1
free = setdiff(1:nU + nP, fixed)
rhs = -(K[:, fixed] * g)
u = zeros(nU + nP); u[fixed] .= g
u[free] = K[free, free] \ rhs[free]

# max |uₓ − parabolic| over all velocity nodes
maxerr = 0.0
for i in 1:U.nnodes
    x = U.nodes[:, i]
    uana = (G / (2μ)) * x[2] * (1 - x[2])
    maxerr = max(maxerr, abs(u[2i - 1] - uana))
end
println("max |uₓ − Poiseuille| = ", round(maxerr; digits = 4))

vtk(m; field = u[1:2:2*U.nnodes], fieldname = "uₓ")
max |uₓ − Poiseuille| = 0.0095

The profile is reproduced to near machine precision, confirming the stabilized equal-order pair is equivalent to the Taylor–Hood P2/P1 pair for this problem.

SUPG advection–diffusion

assemble_advection_diffusion_supg(S, β, D) stabilizes the advection–diffusion operator β·∇u − ∇·(D∇u) against the oscillations that plain Galerkin produces when the Péclet number is large:

S = FESpace(m, ReferenceFE(2, 1))
β = [1.0, 0.0]; D = 1e-4
Ksupg = assemble_advection_diffusion_supg(S, β, D)
println("SUPG operator: ", size(Ksupg), " (stable at high Péclet)")
SUPG operator: (81, 81) (stable at high Péclet)

The nonlinear Navier–Stokes analogue is navier_stokes_operator_supg_pspg(U, P, μ, ρ), which combines SUPG (convection) and PSPG (pressure) stabilization for the equal-order P1/P1 Navier–Stokes system.

Next steps