[Paper Review] Not Killing Sound at Low Mach — An AP IMEX Wave Scheme

Put a low-Mach flow in a compressible code and a faster sound speed sounds like a sharper answer. It is not. Drop the Mach number and a standard upwind scheme wipes a perfectly good wave off the grid. Arun, Das Gupta, and Samantaray (2019) pin that collapse on the linear wave system and build an asymptotic-preserving (AP) scheme with IMEX-RK, stable without reference to the Mach number. The reason the acoustic operator has to be implicit is sitting in the numbers below.

Slide $\varepsilon$ to the left and the red upwind curve climbs with slope $-1$, while the cyan AP curve stays flat.

At $\varepsilon=10^{-3}$ the ratio of the two curves reaches the hundreds. What grows as Mach falls is not the wave. It is the upwind numerical viscosity.

When ε goes to zero, the grid eats the wave#

The scaled isentropic Euler equations carry one small parameter, the Mach number $\varepsilon$. The pressure term enters the momentum equation as $\nabla p / \varepsilon^2$.

$\rho$ is density, $\mathbf{u}$ is velocity, and $p=\rho^\gamma$ is pressure. $\varepsilon$ is reference velocity over reference sound speed.

As $\varepsilon \to 0$ the solution goes to the incompressible limit. That is a singular perturbation. The authors cut the analysis down to the linear wave system.

$\varrho$ is the linearized density perturbation and $a$ is the sound speed of that state. The effective sound speed is $c=a/\varepsilon$. Shrink the Mach number and that acoustic speed blows up.

The well-prepared space — the solutions that must survive#

In the limit $\varepsilon\to 0$ the solution does not wander. It collects in the well-prepared space: constant density, divergence-free velocity.

That space is the kernel of the acoustic operator $L(U)=(a\nabla\cdot\mathbf{u},\,a\nabla\varrho)$. The discrete solution has to land there too. Dellacherie (2010) showed that a scheme which breaks the invariance loses accuracy at low Mach. Keep the kernel and the scheme is asymptotically accurate (AA). Keep the stability limit independent of $\varepsilon$ and it is asymptotic preserving (AP).

Upwind viscosity grows in proportion to the sound speed#

The culprit is upwind numerical viscosity. Diagonalize the wave system into characteristic variables and it splits into two scalar advections at speeds $\pm c$. A Rusanov upwind step on each mode has this amplification.

$\nu$ is the acoustic CFL, $\theta$ is the phase per cell, and $k$ is the wavenumber.

One step barely hurts. The stack is the problem. Reaching a fixed physical time $T$ takes $M = T/\Delta t = T c /(\nu\Delta x)$ acoustic steps. For small $\theta$ the decay per step is $\tfrac{1}{2}\nu(1-\nu)\theta^2$, and $M$ of them produce this log-decay.

Drop the Mach number by a decade and the effective viscosity climbs a decade. The slope $-1$ on the red curve above is that proportion.

Operator splitting — convection explicit, acoustics implicit#

Split the equations by time scale and two operators show up.

$H$ is convection, time scale $O(1)$. $L/\varepsilon$ is acoustics, time scale $O(\varepsilon)$. All of the stiffness sits in $L/\varepsilon$. Integrate that part implicitly. Leave convection explicit. The mix is IMEX (Implicit-Explicit). Whether the same split still holds a time step all the way to Mach 0 is the question in an IMEX TVD reproduction.

Implicit acoustics make central differences stable. The upwind viscosity drops out. $\Delta t$ is set by slow convection, not by fast sound. Inside this model, a smaller $\varepsilon$ is not a reason to shrink the step.

IMEX-RK absorbs the stiffness#

Time integration is IMEX Runge-Kutta (IMEX-RK). The explicit tableau $(\tilde A,\tilde b)$ carries $H$. The implicit tableau $(A,b)$ carries $L/\varepsilon$. At each stage the acoustics condense into an elliptic equation for density, a pressure-Poisson lookalike.

$\beta$ is the diagonal coefficient of the IMEX tableau and $\nabla^2$ is the Laplacian. Invertibility of that elliptic problem comes from saddle-point theory. On the grid, circulant matrices give the same conclusion. A unique solution, uniform stability independent of $\varepsilon$, and invariance of the well-prepared space.

Python — reproducing the amplitude decay of one mode#

The code below computes the upwind amplification of one Fourier mode and prints the amplitude left after physical time $T$, one Mach number at a time. The implicit scheme keeps the mode, so its amplitude stays at 1.

import numpy as np
 
a  = 1.0        # linearization sound speed
k  = 2 * np.pi  # single Fourier mode on [0,1]
T  = 1.0        # fixed physical end time
nu = 0.5        # acoustic CFL number
 
def upwind_amp_factor(eps, N):
    """|g| of one characteristic mode under Rusanov on the wave system."""
    c  = a / eps              # acoustic speed c = a/ε
    dx = 1.0 / N
    th = k * dx               # phase per cell θ = kΔx
    re = 1.0 - nu * (1.0 - np.cos(th))
    im = nu * np.sin(th)
    return np.hypot(re, im)   # |g| <= 1
 
def retained_amplitude(eps, N):
    """Mode amplitude after physical time T (upwind vs implicit)."""
    c  = a / eps
    dx = 1.0 / N
    dt = nu * dx / c          # acoustic CFL step
    M  = T / dt               # number of explicit steps
    explicit = upwind_amp_factor(eps, N) ** M
    implicit = 1.0            # central-implicit preserves the mode
    return M, explicit, implicit
 
if __name__ == "__main__":
    N = 64
    print(f"{'eps':>8} {'steps M':>10} {'explicit':>12} {'AP':>6}")
    for eps in [0.4, 0.2, 0.1, 0.05, 0.02, 0.01]:
        M, ex, ap = retained_amplitude(eps, N)
        print(f"{eps:8.3f} {M:10.0f} {ex:12.3e} {ap:6.2f}")

The output looks like this.

     eps    steps M     explicit     AP
   0.400        320    6.804e-01   1.00
   0.200        640    4.629e-01   1.00
   0.100       1280    2.143e-01   1.00
   0.050       2560    4.591e-02   1.00
   0.020       6400    4.430e-04   1.00
   0.010      12800    1.980e-07   1.00

Drop the Mach $\varepsilon$ slider and the red amplitude hits the floor before physical time $T=1$ is over. Push CFL $\nu$ toward 1 and $\nu(1-\nu)$ shrinks, so the decay eases for a moment. Drop $\varepsilon$ further and that gain disappears under the $1/\varepsilon$ stack. Raising $N$ only partly helps: $\theta^2\propto N^{-2}$ and $\Delta x\propto N^{-1}$ do not cancel.

Plug the code's numbers into the step-count formula and it closes. With $a=1$, $T=1$, $\nu=0.5$, and $N=64$, $\Delta x=1/64$ and $c=1/\varepsilon$. Then $\Delta t=\nu\Delta x/c$, so $M=T/\Delta t=128/\varepsilon$. The printed steps 320, 640, 1280, 2560, 6400, and 12800 are that formula at $\varepsilon=0.4,0.2,0.1,0.05,0.02,0.01$. Halve $\varepsilon$ and $M$ doubles. Cut $\varepsilon$ by ten and $M$ grows by ten.

The log-decay $\tfrac12 T c(1-\nu)k^2\Delta x$ with the same numbers is the exponent $0.154/\varepsilon$. The coefficient comes from $k=2\pi$. At $\varepsilon=0.4$ the exponent is 0.386 and $\mathrm{e}^{-0.386}=0.680$, matching the printed $6.804\times 10^{-1}$ to three digits. At $\varepsilon=0.1$, $\mathrm{e}^{-1.542}=0.214$ against the printed $2.143\times 10^{-1}$. At $\varepsilon=0.01$ the exponent is 15.4. $\mathrm{e}^{-15.4}$ sits near $2.0\times 10^{-7}$, and the table prints $1.980\times 10^{-7}$. The phase $\theta=k\Delta x=2\pi/64\approx 0.098$ is not tiny, so the small-angle step is a little kind at the bottom of the table. And $\nu=0.5$ maximizes $\nu(1-\nu)$ at 0.25. This column is the decay at the CFL where one step hurts most.

AP versus AA#

They are not the same question. AP asks whether the $\varepsilon\to 0$ limit of the discrete scheme is a valid discretization of the limit equation, with a stability restriction independent of $\varepsilon$. AA asks something narrower. On a finite grid, with $\varepsilon$ small but not zero, does the solution stay near the well-prepared space? Upwind fails both. IMEX implicit passes both. The fork is the time discretization, not the spatial order.

Even at the worst ν, the AP column stays 1.00#

Every AP entry in the table is 1.00. The implicit step removes the upwind viscosity, so the exponent $0.154/\varepsilon$ never multiplies the amplitude. $\nu=0.5$ is the CFL where one step decays the most, and after $M=128/\varepsilon$ steps the AP value is still 1.00. For this linear wave system, treating only $L/\varepsilon$ implicitly is the right call.

That call does not travel to every low-Mach system. There are equations where the acoustics are already implicit and the time step shrinks again anyway. The place a second $1/M^2$ hides in the relative velocity is the stiffness left in the enthalpy.