[Paper Review] Fooling a Nonlinear Flux Into Behaving Linearly — Jin–Xin & Suliciu Relaxation

Write an exact Riemann solver (the device that resolves the wave structure at a cell interface) for long enough and the equation of state trips you. The Jacobian of the nonlinear flux $f(u)$ has to be diagonalized every step. A real-gas equation of state has no closed eigenstructure. In 1995 Shi Jin and Zhouping Xin flipped the direction. Don't march the nonlinear equation. Promote the flux to a new unknown, build a linear system, then relax it back. Thomann et al. (2019) pushed that skeleton out to all-speed Euler. The detour is sturdier in one place and brittle in another. That price is the whole note.

Paper — Andrea Thomann, Markus Zenk, Gabriella Puppo, Christian Klingenberg, An all speed second order IMEX relaxation scheme for the Euler equations, arXiv:1907.08398 (2019).

Drop the relaxation-speed slider a below 1 and the Burgers shock — left state 1, right state $-0.4$ — starts to ring.

Near 1.3 the front locks in. Raise eps alone and the wave speed stays put while only the thickness grows.

Why a nonlinear flux makes the Riemann problem hard#

Take a scalar conservation law $u_t + f(u)_x = 0$. Here $u$ is the conserved quantity and $f(u)$ is the flux. Solving a shock tube inside one cell at every face needs the wave speed $f'(u)$ first.

The trouble is real flow. In the compressible Euler equations the flux carries pressure, and pressure is a nonlinear function of density and internal energy. Add a real-gas equation of state and the eigenvalues and eigenvectors no longer come out analytically. Diagonalizing the Jacobian numerically every step is expensive. In the low-Mach regime the acoustic waves scale like $1/M$, so an explicit time step suffocates.

Hence the question: can we defer the nonlinearity instead of facing it head-on every step?

The Jin–Xin promotion — lifting the flux into an unknown#

Their answer: instead of computing the flux $f(u)$ as a value, introduce a new variable $v$ to hold it, and attach a source that slowly drags $v$ toward $f(u)$.

Here $u$ is the original conserved quantity, $v$ is the relaxation variable standing in for the flux, $a$ is a constant relaxation speed, and $\varepsilon$ is the relaxation time. The left-hand side is completely linear. The coefficient matrix has fixed eigenvalues $\pm a$ — the nonlinear diagonalization is gone.

As $\varepsilon \to 0$ the right-hand source forces $v = f(u)$, and the first equation collapses back to the original conservation law. So this $2\times 2$ linear system is a viscous approximation of the nonlinear equation. Because its wave speeds are the constants $\pm a$, the Riemann problem can be solved once, by hand. The figure above is this system on Burgers, $f(u)=u^2/2$.

The sub-characteristic condition — small a, exploding system#

You can't pick the relaxation speed $a$ freely. For a stable approximation it must obey the sub-characteristic condition (Whitham's condition):

Here $f'(u)$ is the true wave speed of the original equation. The condition says the frozen speeds $\pm a$ carried by the relaxation system must always bracket that true speed. For Burgers, $f'(u)=u$, so $a \ge \max|u|$ is enough.

In the chart below, vary the relaxation speed $a$ and the solution range $\max|u|$. The instant the orange curve $f'(u)$ leaves the $\pm a$ band is the threshold of instability.

Drop $a$ below the largest wave speed of the solution and the curve punches through the band. There the relaxation approximation produces negative diffusion instead of physical diffusion. Information flows the wrong way and the system explodes. The next section pins down exactly why.

The hidden diffusion, via Chapman–Enskog#

What happens when $\varepsilon$ is small but not zero? Expand $v = f(u) + \varepsilon v_1 + \cdots$ (a Chapman–Enskog expansion) and substitute into the relaxation system. An effective equation drops out:

The right-hand side is the diffusion the relaxation quietly injects. Its coefficient is $a^2 - f'(u)^2$. For that to stay non-negative you need exactly $a \ge |f'(u)|$. So the sub-characteristic condition was really the condition that the relaxation's artificial viscosity never goes negative.

The trade-off is now visible. A large $a$ gives comfortable stability but a big $a^2 - f'^2$ that smears the solution. Pinning $a$ right at the sub-characteristic limit is sharp but risky. And $\varepsilon$ sets the overall magnitude of the diffusion — which is why raising $\varepsilon$ in the earlier simulation thickens the shock.

The margin is the viscosity. The Riemann data in the code are $u=1$ on the left and $u=-0.4$ on the right, so $\max|u|$ starts at 1 and the floor is 1. At $a=1.3$, $a^2=1.69$. On the left, $f'=1$ and the diffusion coefficient is $1.69-1=0.69$. On the right, $f'=-0.4$, so $1.69-0.16=1.53$. Same $\varepsilon=10^{-4}$, and the right state still smears more. At $a=0.8$, $a^2=0.64$ and the left coefficient is $0.64-1=-0.36$. The sign flips. Information runs back through the front, and the amplitude does not stay at the initial value 1.

Suliciu — relax only the pressure, not the whole flux#

Jin–Xin relaxes every component of the flux. Elegant, but the diffusion is excessive. What practitioners reach for more often is Suliciu-type relaxation. In the Euler equations the only genuine nonlinearity is the pressure, so you relax only the pressure into a new variable $\pi$.

Here $\rho$ is density, $u$ velocity, $p$ the actual pressure, and $\pi$ the relaxed pressure. The sub-characteristic condition becomes $a > \rho\sqrt{\partial_\rho p}$, bracketing the sound speed. The relaxed system's eigenvalues are $u,\ u\pm a/\rho$, all linearly degenerate — only contact-like waves remain, which are easy to handle.

This is where Thomann et al. (2019) begin. Aiming at the low-Mach regime, they split the pressure into a slow part and a fast acoustic part. The slow part is treated explicitly; the fast acoustic part is solved implicitly on top of the relaxation system. They add a new velocity variable $\hat u$ to secure Mach-number-independent diffusion, and the result is an asymptotic-preserving, second-order IMEX scheme that converges to incompressible Euler in the low-Mach limit. Thanks to the linearly degenerate structure, all of this runs without a single nonlinear diagonalization.

Python — solving Burgers through the relaxation system#

Let's reproduce the skeleton of the relaxation idea in numpy: Burgers' equation through the Jin–Xin system. Splitting into characteristic variables $r,s$ turns the left side into two simple advections.

import numpy as np
 
def burgers_flux(u):
    return 0.5 * u * u                     # f(u) = u^2 / 2
 
def relaxed_step(u, v, a, dx, dt, eps):
    # characteristic split: r moves at +a, s moves at -a
    r = 0.5 * (u + v / a)
    s = 0.5 * (u - v / a)
    nu = a * dt / dx                       # CFL number a*dt/dx
    # first-order upwind (periodic, via np.roll)
    r_new = r - nu * (r - np.roll(r, 1))   # right-moving wave
    s_new = s + nu * (np.roll(s, -1) - s)  # left-moving wave
    u_new = r_new + s_new
    v_new = a * (r_new - s_new)
    # relaxation source: pull v toward f(u) (exponential integrator)
    kappa = 1.0 - np.exp(-dt / eps)
    v_new += (burgers_flux(u_new) - v_new) * kappa
    return u_new, v_new
 
def run_relaxation(u0, a, eps, cfl=0.9, t_end=0.3):
    n = u0.size
    dx = 2.0 / n
    u = u0.copy()
    v = burgers_flux(u)                    # start on the equilibrium manifold
    dt = cfl * dx / a                      # sub-characteristic sets the CFL
    t = 0.0
    while t < t_end:
        u, v = relaxed_step(u, v, a, dx, dt, eps)
        t += dt
    return u
 
# Riemann initial data: left 1.0, right -0.4 (a shock)
n = 400
x = np.linspace(-1, 1, n, endpoint=False) + 1.0 / n
u0 = np.where(x < 0.0, 1.0, -0.4)
 
u_ok = run_relaxation(u0, a=1.3, eps=1e-4)   # a >= max|u|=1.0  -> stable
u_bad = run_relaxation(u0, a=0.8, eps=1e-4)  # a <  max|u|      -> violated
 
print("a=1.3  peak |u| =", round(float(np.max(np.abs(u_ok))), 3))   # ~1.0
print("a=0.8  peak |u| =", round(float(np.max(np.abs(u_bad))), 3))  # blows up

The output testifies directly to the sub-characteristic condition. At $a=1.3$ the peak amplitude stays near the initial value and a clean shock forms. At $a=0.8$ the amplitude jumps several-fold. No nonlinear Jacobian was diagonalized. Only linear advection at the constant speeds $\pm a$, repeated.

The step count is tied to $a$ as well. With $n=400$ on an interval of length 2, $\Delta x=0.005$. At CFL 0.9, $\Delta t=0.9\times 0.005/a$. For $a=1.3$, $\Delta t\approx 3.46\times 10^{-3}$, and $t=0.3$ is about 87 of those steps. For $a=0.8$, $\Delta t\approx 5.63\times 10^{-3}$, a coarser march. Lowering $a$ loosens the time step. With $\varepsilon=10^{-4}$, $\Delta t/\varepsilon$ exceeds 30, so the exponential in $\kappa=1-\exp(-\Delta t/\varepsilon)$ vanishes. $\kappa=1$. Every step pins $v$ to $f(u)$. Whatever thickness remains is only $\varepsilon(a^2-f'^2)$.

A critical look#

Relaxation isn't free. Artificial diffusion $\varepsilon(a^2-f'^2)$ always tags along, and the safer (larger) you make $a$, the more the solution smears. Enforcing the sub-characteristic condition needs a global bound on $\max|f'(u)|$, and near strong shocks or vacuum a conservative bound overshoots the diffusion. The logic that was clean for scalars demands, for real-gas Euler, a local estimate of $\rho\sqrt{\partial_\rho p}$ plus extra machinery to keep pressure and density positive. Thomann et al.'s IMEX coupling has to solve an elliptic equation for the acoustic part, so there's a crossover where the low-Mach gain is offset by the linear-solver cost. Reproducing it, the joint tuning of $\varepsilon$ and $a$ turns out to be touchier than expected.

OpenFOAM and Fluent's density-based solvers don't ship this relaxation Riemann solver directly, but the HLLC and AUSM families of approximate Riemann solvers share the same philosophy — simplify the wave structure to avoid diagonalization. The Suliciu solver is implemented in open-source codes like SU2 as an option that captures contact discontinuities exactly.

Where an a pinned just above the floor breaks#

The floor for this data is 1. The code picks 1.3. That 0.3 of headroom is the 0.69 diffusion coefficient on the left state. Shave the margin to zero, set $a=\max|u|$, and the coefficient on the front is 0. The shock looks sharp. One rounding nudge below the floor and the sign flips.

The miss that shows up later is freezing the bound on the first step. Once $|u|$ climbs past 1, the $a=1.3$ that looked safe is already outside the sub-characteristic condition. Suliciu moves the same test to $a>\rho\sqrt{\partial_\rho p}$. If the sound speed jumps from cell to cell and the bound lags by one update, the relaxation diffusion is already negative, implicit acoustics or not. Choosing to touch only the pressure sits in the same place as splitting advection off the pressure flux.