Skip to content

Stabilised Navier–Stokes: equal-order P1/P1 (SUPG + PSPG)

Download .py All tutorials

Taylor–Hood (P2 velocity / P1 pressure) is inf-sup stable, and it is what the other flow tutorials use. It is also expensive: in 3-D a P2 velocity carries roughly four times the DOFs of P1 on the same mesh. Equal-order P1/P1 is not a stable pair — used naively its pressure carries a spurious mode — but residual-based stabilisation compensates, and in jNO that stabilisation is not a feature you enable. It is terms you write.

\[(\mathbf u\!\cdot\!\nabla)\mathbf u - \nu\,\Delta\mathbf u + \nabla p = 0,\qquad \nabla\!\cdot\mathbf u = 0\]

The exact solution is Kovasznay (1948), so every number here is checked against a closed form.

The two pieces the formula needs

Neither is a stabilisation feature; they are what the formula is made of.

A vector Laplacian. The momentum strong residual is \(\mathbf r_M = (\mathbf u\!\cdot\!\nabla)\mathbf u - \nu\Delta\mathbf u + \nabla p\), so jno.np.laplacian has to accept a vector field.

An element metric. \(\tau\) has to know the element's shape, not just its size — on a boundary-layer cell the streamwise and wall-normal directions are not interchangeable. dom.cell_size is \(|\det J|^{1/d}\), an isotropic scalar that cannot see stretch at all. dom.cell_metric is the covariant metric \(G = J^{-\mathsf T}J^{-1}\), per quadrature point:

G  = d.cell_metric
gG = lambda a: inner(a, inner(G, a, n_contract=1), n_contract=1)      # aᵀ G a
tau = jno.lag((gG(ub) + C_I * NU**2 * inner(G, G, n_contract=2)) ** -0.5)

trace(G) and inner(G, G, n_contract=2) read as written — trace takes the last two axes and inner aligns on the trailing ones.

jno.lag on tau

\(\tau\) is a coefficient, not part of the equation. Differentiating through \(\mathbf u^\mathsf T G\mathbf u\) hands Newton a tangent it gains nothing from. jno.lag freezes it within each linearisation.

SUPG and PSPG are two terms, and their signs differ

This is the part that is easy to get wrong, so it is worth stating plainly.

They are the same strong residual weighted by two different test perturbations — but they are two separate terms in the list, because an additive term carries exactly one test field: the test field is what names the equation block. And each term joins an equation that already has a sign convention, so it must carry that convention:

r_m = adv(gu, ub) - NU * lap(u, [xi, yi]) + gp        # the strong residual

momentum   = momentum   + tau * inner(adv(gv, ub), r_m, n_contract=1)   # SUPG → momentum   (+)
continuity = continuity - tau * inner(gq, r_m, n_contract=1)            # PSPG → continuity (−)

Momentum is written \(+(\mathbf u\!\cdot\!\nabla\mathbf u,\mathbf v)\), so SUPG is +. Continuity is written \(-(q,\nabla\!\cdot\mathbf u)\), so PSPG is −. Getting that sign wrong is not a small error: with + on PSPG, Newton diverges outright.

What it buys, measured

Kovasznay at \(\mathrm{Re}=20\), equal-order P1/P1, against the closed form:

velocity error pressure error
no stabilisation, ms=0.11 2.53e-02 4.09e-01
SUPG + PSPG, ms=0.11 3.07e-02 9.42e-02

A 4.3× better pressure — which is the whole point, since the pressure is what the unstable pair gets wrong — for a slightly worse velocity. That trade is expected and is why the term exists. Under refinement the observed rates are 1.78 (velocity) / 1.70 (pressure), in line with theory for a PSPG-stabilised equal-order pair.

Scope — stated up front

  • tau here is the steady form. A transient run adds a \((2/\Delta t)^2\) term inside the root. That term assumes the fixed grid from domain(time=…), so it is wrong under jno.solve.adaptive(), where the step size is chosen per step.
  • grad-div / LSIC is not included. The usual coefficient \(\tau_C = 1/(\mathrm{tr}(G)\,\tau_M)\) made both errors worse at this Reynolds number (velocity 6.14e-02 → 8.36e-02, pressure 1.78e-01 → 2.06e-01), so its calibration is left open rather than shipped as though it were verified. The term is one line if you want it: tau_c * div(gu) * div(gv).
  • Use the sparse-direct Newton, jno.solve.newton(direct=True). The matrix-free Newton (newton(direct=False), the default when this tutorial was written) goes NaN on a cold start from rest here.
  • dom.cell_metric is native 2-D/3-D volume terms only — like dom.cell_size. A 1-D form or a non-nodal element family refuses by name.
  • This tutorial does not claim a turbulence model, a free surface, or a compressible path. jNO's fluid scope is still laminar incompressible; see Formulations.

Full script

"""**Stabilised equal-order Navier-Stokes** -- P1 velocity and P1 pressure, which is *not* an
inf-sup-stable pair, made to work by residual-based stabilisation written as ordinary weak-form terms.

    (u.grad)u - nu lap u + grad p = 0,   div u = 0

The exact solution is Kovasznay (1948), so every number below is checked against a closed form.

Two things make this possible, and neither is a "stabilisation feature" -- they are the two pieces the
formula is built out of:

  * `jno.np.laplacian` on a **vector** field, so the momentum strong residual can carry `nu*lap(u)`;
  * `dom.cell_metric`, the element metric `G = J^-T J^-1`, which is what `tau` needs in order to be
    direction-aware. `dom.cell_size` is an isotropic scalar and cannot see a stretched cell.

Everything else is the FEM contract: the term list is the whole problem.

Why bother with equal order? Taylor-Hood needs P2 velocity. In 3-D that is roughly four times the
velocity DOFs of P1 for the same mesh. PSPG buys the cheaper pair.

References:
  tau, and the metric it is built on -- Tezduyar & Osawa, *CMAME* **190** (2000), Sec. 3
  PSPG                              -- Hughes, Franca & Balestra, *CMAME* **59** (1986)
  the exact solution                -- Kovasznay, *Proc. Camb. Phil. Soc.* **44** (1948) 58
"""

import os

os.environ.setdefault("XLA_PYTHON_CLIENT_MEM_FRACTION", "0.6")

import jax

jax.config.update("jax_enable_x64", True)

import numpy as np  # noqa: E402
from shapely.geometry import box  # noqa: E402

import jno  # noqa: E402

inner, grad, trace, lap = jno.np.inner, jno.np.grad, jno.np.trace, jno.np.laplacian

NU = 0.05
C_I = 36.0  # inverse-estimate constant for linear elements (Tezduyar & Osawa, Sec. 3)
RE = 1.0 / NU
LAM = RE / 2.0 - np.sqrt(RE**2 / 4.0 + 4.0 * np.pi**2)

# The exact Kovasznay fields, for the boundary data and the final check.
u_ex = lambda x, y: 1.0 - np.exp(LAM * x) * np.cos(2 * np.pi * y)  # noqa: E731
v_ex = lambda x, y: LAM / (2 * np.pi) * np.exp(LAM * x) * np.sin(2 * np.pi * y)  # noqa: E731
p_ex = lambda x, y: 0.5 * (1.0 - np.exp(2 * LAM * x))  # noqa: E731

d = jno.domain(box(-0.5, -0.5, 1.0, 1.5), mesh_size=0.075)
d.point_region("ppin", (-0.5, -0.5))  # one pressure gauge DOF

u, v = d.fem_symbols(value_shape=(2,), names=("u", "v"), order=1)  # P1 velocity  -- equal order,
p, q = d.fem_symbols(names=("p", "q"), order=1)  # P1 pressure  -- NOT inf-sup stable
xi, yi = d.variable("interior", split=True)[:2]
xb, yb = d.variable("boundary", split=True)[:2]
xpn, ypn = d.variable("ppin", split=True)[:2]

ub, vv = u.bind(x=xi, y=yi), v.bind(x=xi, y=yi)
gu, gv = grad(u, [xi, yi]), grad(v, [xi, yi])
gp, gq = grad(p, [xi, yi]), grad(q, [xi, yi])
pp, qq = p.bind(x=xi, y=yi), q.bind(x=xi, y=yi)

# Notation first, so the weak form below reads like the paper.
div = lambda gw: trace(gw)  # noqa: E731  -- div w
adv = lambda gw, w: inner(gw, w, n_contract=1)  # noqa: E731  -- (w.grad)w

# --- Galerkin ------------------------------------------------------------------------------------
momentum = inner(adv(gu, ub), vv, n_contract=1) + NU * inner(gu, gv, n_contract=2) - pp * div(gv)
continuity = -qq * div(gu)

# --- stabilisation -------------------------------------------------------------------------------
G = d.cell_metric  # the element metric, per quadrature point
gG = lambda a: inner(a, inner(G, a, n_contract=1), n_contract=1)  # noqa: E731  -- a^T G a

# tau is LAGGED: it is a coefficient, not part of the equation, and differentiating through `u.G u`
# gives Newton a tangent it does not need to follow.
tau = jno.lag((gG(ub) + C_I * NU**2 * inner(G, G, n_contract=2)) ** -0.5)

r_m = adv(gu, ub) - NU * lap(u, [xi, yi]) + gp  # the momentum STRONG residual

# SUPG and PSPG are the same strong residual against two different test perturbations. They are two
# TERMS, not one: an additive term carries exactly one test field, because the test field is what
# names the equation block. And each carries the sign of the equation it joins -- momentum is written
# `+(u.grad u, v)` so SUPG is `+`; continuity is written `-(q, div u)` so PSPG is `-`.
momentum = momentum + tau * inner(adv(gv, ub), r_m, n_contract=1)  # SUPG  -> momentum block
continuity = continuity - tau * inner(gq, r_m, n_contract=1)  # PSPG  -> continuity block

bx = 1.0 - jno.np.exp(LAM * xb) * jno.np.cos(2 * np.pi * yb)
by = LAM / (2 * np.pi) * jno.np.exp(LAM * xb) * jno.np.sin(2 * np.pi * yb)

fem = jno.fem(
    [
        momentum,
        continuity,
        u(xb, yb)[0] - bx,
        u(xb, yb)[1] - by,
        p(xpn, ypn) - float(p_ex(-0.5, -0.5)),
    ]
)
assert not fem.is_linear, "the convective term keeps this nonlinear"
print(f"\nStabilised P1/P1 Kovasznay flow (Re={RE:.0f}): dofs={fem.dofs}")

# Sparse-direct Newton: the matrix-free Newton (newton(direct=False)) goes NaN on this cold start from rest.
sol = np.asarray(fem.solve(nonlinear=jno.solve.newton(direct=True, rtol=1e-8, atol=1e-8)))

off = fem.offsets
uv = sol[off[0] : off[1]].reshape(-1, 2)
ph = sol[off[1] :].reshape(-1)
pts = np.asarray(d.mesh.points)[:, :2]
x, y = pts[:, 0], pts[:, 1]

e_u = np.sqrt(np.mean((uv[:, 0] - u_ex(x, y)) ** 2 + (uv[:, 1] - v_ex(x, y)) ** 2))
pe = p_ex(x, y)
e_p = np.sqrt(np.mean(((ph - ph.mean()) - (pe - pe.mean())) ** 2))  # pressure is gauge-free

print(f"  velocity RMS error vs Kovasznay: {e_u:.3e}")
print(f"  pressure RMS error vs Kovasznay: {e_p:.3e}")
assert e_u < 5e-2, f"velocity error {e_u:.3e} is too large -- the stabilised pair should resolve this"
assert e_p < 1.5e-1, f"pressure error {e_p:.3e} is too large -- PSPG should be controlling it"
print("\nEqual-order P1/P1 converged against a closed-form solution.")