Catching Shocks Without Multiplying a Matrix at Every Face

At Mach 0.01 the polynomial split alone lets pressure chatter on the grid, so the Mp term has to pull that pressure difference into the mass flux.

In 1993, Meng-Sing Liou at NASA Lewis was unhappy with the matrix work inside the Roe solver. Shocks came out fine. But solving every single face meant decomposing the full eigenstructure of a 5×5 Jacobian. Liou took another road. He split the flux into "convective" and "pressure" parts and separated each with polynomials of the Mach number. The result, AUSM (Advection Upstream Splitting Method), never multiplies a matrix. This post follows that split through the math, then solves the Shu–Osher problem in Python to check it.

Pushing the Mach slider from near zero out past unity shows where the transonic polynomial joins the supersonic upwind branch.

M = 0.40β = 0.125

Cyan = forward split (M+ / P+), pink = backward split (M- / P-). Outside the sonic band |M|>1 the curves collapse onto the fully upwind branch; inside, the polynomial blend keeps them smooth and differentiable.

Outside $|M|>1$ the two curves snap onto the fully upwind branch. Inside, the polynomial joins them smoothly and keeps the function differentiable. Raise β to bulge the curvature in the transonic band. That smoothness is exactly what makes convergence possible without a Roe solver.

1993: Liou splits the flux in two#

Look at the face flux of the Euler equations. In one dimension,

\mathbf{F} = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ (E+p)u \end{pmatrix} = \underbrace{\rho u \begin{pmatrix} 1 \\ u \\ H \end{pmatrix}}_{\text{convective}} + \underbrace{\begin{pmatrix} 0 \\ p \\ 0 \end{pmatrix}}_{\text{pressure}}

Here $\rho$ is density, $u$ velocity, $p$ pressure, $E$ total energy, and $H=(E+p)/\rho$ the total enthalpy.

Liou's insight was simple. The convective term is decided by which direction the mass flux $\dot m = \rho u$ comes from. The pressure term rides acoustic waves in both directions. Different physics, so handle them separately. Where Roe resolved the difference of two states through a Jacobian (flux-difference splitting), AUSM divides each state's contribution to the face on its own (flux-vector splitting).

The AUSM flux at the interface is assembled like this.

\mathbf{F}_{1/2} = \dot m_{1/2}\,\boldsymbol{\Psi}_{L/R} + \mathbf{p}_{1/2}

$\boldsymbol{\Psi}=(1,\,u,\,H)^{T}$ is what the convection carries, and $\mathbf{p}_{1/2}=(0,\,p_{1/2},\,0)^{T}$ is the pressure piece. Everything hinges on how we choose the face mass flux $\dot m_{1/2}$ and the face pressure $p_{1/2}$ .

Splitting the Mach number into polynomials#

The face mass flux comes from the face Mach number $M_{1/2}$ .

M_{1/2} = \mathcal{M}^{+}_{(4)}(M_L) + \mathcal{M}^{-}_{(4)}(M_R)

$M_L=u_L/a_{1/2}$ and $M_R=u_R/a_{1/2}$ are the left and right cell Mach numbers, and $a_{1/2}$ is the face sound speed. The split functions $\mathcal{M}^{\pm}$ are the heart of the method.

The simplest first-degree split is plain upwinding.

\mathcal{M}^{\pm}_{(1)}(M) = \tfrac{1}{2}\left(M \pm |M|\right)

In supersonic flow ( $|M|\ge 1$ ) we use this directly. One side goes to zero, giving full upwinding. The trouble is near sonic. The first-degree function kinks at $M=0$ where its derivative is discontinuous, and solver convergence stalls there.

AUSM+ fills the transonic band smoothly with a fourth-degree polynomial ( $\beta=1/8$ ).

\mathcal{M}^{\pm}_{(4)}(M) = \begin{cases} \tfrac{1}{2}(M\pm|M|) & |M|\ge 1 \\ \pm\tfrac{1}{4}(M\pm 1)^2 \pm \beta(M^2-1)^2 & |M|<1 \end{cases}

Pressure is split the same way ( $\alpha=3/16$ ).

\mathcal{P}^{\pm}_{(5)}(M) = \begin{cases} \tfrac{1}{2}(1\pm\operatorname{sgn}M) & |M|\ge 1 \\ \tfrac{1}{4}(M\pm 1)^2(2\mp M) \pm \alpha M(M^2-1)^2 & |M|<1 \end{cases}

Filling the transonic band this way removes the kink at $M=0$ .

Pressure with pressure, convection with convection#

Once the split functions are ready, the face value gathers in one line. The mass flux is

\dot m_{1/2} = a_{1/2}\,M_{1/2}\, \begin{cases}\rho_L & M_{1/2}>0 \\ \rho_R & M_{1/2}\le 0\end{cases}

The sign is the upwind direction. If $M_{1/2}>0$ , the left cell's $\rho,u,H$ ride out. The pressure is a weighted sum of both sides.

p_{1/2} = \mathcal{P}^{+}_{(5)}(M_L)\,p_L + \mathcal{P}^{-}_{(5)}(M_R)\,p_R

Convection comes from one side only; pressure comes from both. That asymmetry is why AUSM keeps a contact discontinuity (density jumps but pressure stays flat) sharp. Roe's approximate solver treats the contact through one eigenvalue of the Jacobian, but AUSM separates it naturally at the point where the mass flux vanishes. No extra dissipation sneaks in to rattle the pressure.

AUSM+ and AUSM+-up — all the way down to low Mach#

The original AUSM had a weakness. In the nearly incompressible limit, with Mach as low as 0.01, the pressure field oscillated cell by cell. The acoustic term was dissipated too strongly relative to momentum. In 2006 Liou patched this band with AUSM+-up, adding a pressure-diffusion term to the mass flux and a velocity-diffusion term to the pressure.

M_{1/2} = \mathcal{M}^{+}_{(4)}(M_L) + \mathcal{M}^{-}_{(4)}(M_R) + M_p

M_p = -\frac{K_p}{f_a}\,\max\!\left(1-\sigma \bar{M}^2,\,0\right)\frac{p_R-p_L}{\rho_{1/2}\,a_{1/2}^2}

$M_p$ pulls the pressure difference into the mass flux and damps the low-Mach oscillation. $f_a$ is a scaling function that slides between 0 and 1 with the local Mach number. At high Mach, $f_a\to 1$ , the extra term vanishes, and the original shock-capturing returns untouched. One formula handles both the incompressible limit and hypersonic flow. That is why AUSM+-up is called an "all-speed" scheme.

Python — testing it on Shu–Osher#

For the verification problem we pick Shu–Osher. A Mach 3 shock charges into a sine-wave density disturbance — a one-dimensional Euler test that probes shock capturing and high-frequency preservation at once. The code below solves it with first-order AUSM+.

import numpy as np
 
GAMMA = 1.4
 
def mach_split(M, sign):
    # AUSM+ fourth-degree split Mach polynomial (beta = 1/8)
    beta = 1.0 / 8.0
    if abs(M) >= 1.0:
        return 0.5 * (M + abs(M)) if sign > 0 else 0.5 * (M - abs(M))
    if sign > 0:
        return 0.25 * (M + 1.0) ** 2 + beta * (M * M - 1.0) ** 2
    return -0.25 * (M - 1.0) ** 2 - beta * (M * M - 1.0) ** 2
 
def pressure_split(M, sign):
    # AUSM+ fifth-degree split pressure polynomial (alpha = 3/16)
    alpha = 3.0 / 16.0
    if abs(M) >= 1.0:
        return 0.5 * (1.0 + np.sign(M)) if sign > 0 else 0.5 * (1.0 - np.sign(M))
    if sign > 0:
        return 0.25 * (M + 1.0) ** 2 * (2.0 - M) + alpha * M * (M * M - 1.0) ** 2
    return 0.25 * (M - 1.0) ** 2 * (2.0 + M) - alpha * M * (M * M - 1.0) ** 2
 
def ausm_plus(wl, wr):
    rl, ul, pl = wl
    rr, ur, pr = wr
    al = np.sqrt(GAMMA * pl / rl)
    ar = np.sqrt(GAMMA * pr / rr)
    a12 = 0.5 * (al + ar)                    # face sound speed (simple mean)
    Ml, Mr = ul / a12, ur / a12
    M12 = mach_split(Ml, +1) + mach_split(Mr, -1)
    p12 = pressure_split(Ml, +1) * pl + pressure_split(Mr, -1) * pr
    Hl = (pl / (GAMMA - 1) + 0.5 * rl * ul * ul + pl) / rl
    Hr = (pr / (GAMMA - 1) + 0.5 * rr * ur * ur + pr) / rr
    mdot = a12 * M12 * (rl if M12 > 0 else rr)
    psi = np.array([1.0, ul, Hl]) if M12 > 0 else np.array([1.0, ur, Hr])
    return mdot * psi + np.array([0.0, p12, 0.0])
 
def euler_step(U, dx, cfl):
    r = U[0]
    u = U[1] / r
    p = (GAMMA - 1.0) * (U[2] - 0.5 * r * u * u)
    a = np.sqrt(GAMMA * p / r)
    dt = cfl * dx / np.max(np.abs(u) + a)
    n = U.shape[1]
    F = np.zeros((3, n + 1))
    for f in range(1, n):                    # interior face fluxes
        F[:, f] = ausm_plus((r[f-1], u[f-1], p[f-1]), (r[f], u[f], p[f]))
    F[:, 0] = F[:, 1]                         # transmissive boundary
    F[:, n] = F[:, n - 1]
    return U - (dt / dx) * (F[:, 1:] - F[:, :-1]), dt
 
def shu_osher(n=400, t_end=1.8):
    x = np.linspace(-5, 5, n, endpoint=False) + 5.0 / n
    dx = 10.0 / n
    r = np.where(x < -4, 3.857143, 1.0 + 0.2 * np.sin(5.0 * x))
    u = np.where(x < -4, 2.629369, 0.0)
    p = np.where(x < -4, 10.33333, 1.0)
    U = np.zeros((3, n))
    U[0], U[1], U[2] = r, r * u, p / (GAMMA - 1.0) + 0.5 * r * u * u
    t = 0.0
    while t < t_end:
        U, dt = euler_step(U, dx, cfl=0.4)
        t += dt
    return x, U[0]
 
if __name__ == "__main__":
    x, rho = shu_osher()
    print(f"density min/max : {rho.min():.3f} / {rho.max():.3f}")
    print(f"density at x=2.0 : {rho[np.argmin(np.abs(x - 2.0))]:.3f}")

Running it gives a density field where the sine wave is compressed and survives behind the shock. After the Mach 3 shock passes, the peak density climbs past 4. The split functions alone capture the shock stably, while the pressure stays smooth and free of oscillation.

The left state is the post-shock side of Mach 3 at gamma 1.4. Pressure ratio: $2 \times 1.4 / 2.4 \times 8 = 9.333$ , plus 1 is 10.333, the code's 10.33333. Density ratio: $2.4 \times 9 / 5.6 = 3.857$ , the code's 3.857143. Velocity 2.629369 follows that shock. Right state: sine amplitude 0.2, wave number 5, pressure 1.

beta is 1/8. At Mach 0 the branches are $0.375$ and $-0.375$ , sum 0. At absolute Mach 1 the polynomial meets the upwind branch. 400 cells over length 10 give spacing 0.025. CFL 0.4, stop time 1.8. alpha is 3/16.

The same AUSM+ flux applied to a 1D Sod shock tube runs live below.

CFL = 0.60

Sod shock tube solved live with AUSM+. Watch three structures appear from one diaphragm: a rarefaction fan moving left, a contact discontinuity in the middle (sharp in density, flat in pressure), and a shock on the right. Push CFL toward 0.95 to see the first-order scheme stay stable but smear the contact further.

Three structures split out of one diaphragm: a rarefaction fan spreading left, a contact discontinuity in the middle (density bends, pressure stays flat), and a shock on the right. Even pushed to CFL 0.95 the first-order scheme never blows up. It just smears the contact further.

Standing it next to Roe#

AUSM is not free. The first-order version smears the contact a little more than Roe. Using the raw cell values, as in the code above, stays first-order accurate. In practice you apply MUSCL or WENO reconstruction to the left and right states first to lift it to higher order. The split functions stay the same; only the inputs change.

Item	Roe (FDS)	AUSM+ (FVS)
Cost per face	Jacobian eigendecomposition	Polynomial evaluation
Contact discontinuity	Very sharp	Sharp (no entropy fix)
Carbuncle phenomenon	Fragile	Robust
Low-Mach behavior	Needs separate preconditioning	Built into AUSM+-up

When the grid aligns normal to a hypersonic shock, Roe produces an unphysical bulge called the carbuncle. The AUSM family is far more robust against this illness. The very structure of separating convection and pressure feeds less energy into the transverse instability of the shock front.

Drop Mp at Mach 0.01 and the pressure chatter comes back#

If the shock is the whole question and the Mach number is not small, AUSM+ with beta 1/8 and alpha 3/16 is enough. There is no Jacobian. The contact splits where the mass flux hits zero.

At a near-incompressible Mach of 0.01, the polynomials alone let the pressure field chatter at the grid spacing. $M_p$ has to pull the left-right pressure difference into the mass flux, and that term has to switch off as the Mach number rises. Another cut between pressure and convection is the Toro–Vázquez split, which moves the sound speed into the pressure subsystem.

Catching Shocks Without Multiplying a Matrix at Every Face — AUSM Flux Splitting