[Paper Review] Second-Order Accuracy in a Single Stage — A GRP Solver for the Kapila Multiphase Model

A time derivative at the interface buys second-order time accuracy without another Runge-Kutta stage.

Second-order time accuracy was supposed to mean several Runge-Kutta stages. This paper gets the order from a single stage. The tool is a GRP (Generalized Riemann Problem) solver that returns the interface value and its time derivative together.

Paper: Tuowei Chen, Zhifang Du, Generalized Riemann Problem Method for the Kapila Model of Compressible Multiphase Flows, arXiv:2310.08241 (2023).

Turn the GRP check off and the slope slider locks, while the interface value freezes in time. Raise CFL and only the half-step mark that enters the flux moves forward.

advection speed a: 1.00left slope uₓ: 0.50CFL: 0.80GRP (off = Godunov)

Set the slope to zero and the cyan dot stops, sitting on top of Godunov. Flip the sign of the advection speed $a$ and the cell that supplies the value switches. The amber ring is the state at $\Delta t/2$ .

Second-order in one stage — the problem GRP solves#

One finite-volume step is fixed by the interface flux $\mathbf{F}_{i+1/2}$ . Godunov solves a Riemann problem and gets one number. That number is constant in time. Second order in time then means repeating the same spatial operator, Runge-Kutta style.

GRP evaluates the flux at the half time $t^{n+1/2}$ .

\mathbf{U}(\mathbf{x}_{i+1/2}, t^{n+1/2}) = \mathbf{U}^{n,*}_{i+1/2} + \frac{\Delta t}{2}\left(\frac{\partial \mathbf{U}}{\partial t}\right)^{n,*}_{i+1/2}

$\mathbf{U}^{n,*}$ is the Riemann state at the interface. $(\partial \mathbf{U}/\partial t)^{n,*}$ is its instantaneous time derivative. Added as a Taylor step, the interface state moves $\Delta t/2$ into the future. One evaluation.

From piecewise-constant to piecewise-linear — Riemann vs GRP#

A Riemann problem has piecewise-constant data on each side. A GRP problem has piecewise-linear data. Waves leaving the jump are no longer straight, and the solver returns a time derivative with the Riemann state.

For scalar advection $u_t + a u_x = 0$ , call the slope in cell $i$ by $\sigma_i$ . If $a>0$ , the interface Riemann state is the extrapolation from the left cell.

U^{*}_{i+1/2} = u_i + \sigma_i \frac{\Delta x}{2}

$u_i$ is the cell average, $\sigma_i$ the in-cell slope, $\Delta x$ the cell width. The cyan dot above is this expression. Setting the slope to zero, and falling back to Godunov, is the same expression.

The Lax-Wendroff procedure — turning a time derivative into a spatial one#

The time derivative comes out of the equation itself.

\frac{\partial u}{\partial t} = -a \frac{\partial u}{\partial x} = -a\,\sigma_i

The left-hand side becomes the spatial derivative on the right. The reconstruction already knows $\sigma_i$ . The half-step interface value is

U^{n+1/2}_{i+1/2} = U^{*}_{i+1/2} - a\,\sigma_i \frac{\Delta t}{2}

and the flux is $F = a\,U^{n+1/2}$ . For Euler or Kapila the same step is $\partial u/\partial t = -A\,\partial u/\partial x$ , with Jacobian $A$ . The paper linearizes the acoustic waves this way and solves the GRP approximately.

Python — measuring the convergence order of a scalar GRP solver#

A smooth initial field $u_0=\sin(2\pi x)$ is advected for one period. The $L_1$ orders of Godunov and GRP are compared on periodic linear advection.

import numpy as np
 
# 1D linear advection  u_t + a u_x = 0,  periodic domain [0,1]
# First-order Godunov vs single-stage second-order GRP (Lax-Wendroff)
 
def cell_gradient(u, dx):
    # central (unlimited) slope so the smooth-solution order is clean
    return (np.roll(u, -1) - np.roll(u, 1)) / (2.0 * dx)
 
def godunov_step(u, a, dx, dt):
    # piecewise constant: upwind interface value (a > 0)
    flux = a * u                       # F_{i+1/2} = a u_i
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def grp_step(u, a, dx, dt):
    g = cell_gradient(u, dx)           # spatial slope sigma_i
    u_star = u + g * dx / 2.0          # interface Riemann solution U*
    u_half = u_star - a * g * dt / 2.0 # Lax-Wendroff half-step U^{n+1/2}
    flux = a * u_half                  # F_{i+1/2}^{n+1/2}
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def l1_error(num, exact, dx):
    return np.sum(np.abs(num - exact)) * dx
 
def convergence_study(a=1.0, cfl=0.4, t_end=1.0):
    print(f"{'N':>5} {'Godunov':>16} {'GRP':>16}")
    prev = {}
    for N in [40, 80, 160, 320]:
        dx = 1.0 / N
        x = (np.arange(N) + 0.5) * dx
        u0 = np.sin(2 * np.pi * x)     # smooth initial condition
        dt = cfl * dx / abs(a)
        nsteps = int(round(t_end / dt))
        dt = t_end / nsteps            # land exactly on t_end
        ug = u0.copy(); ur = u0.copy()
        for _ in range(nsteps):
            ug = godunov_step(ug, a, dx, dt)
            ur = grp_step(ur, a, dx, dt)
        exact = np.sin(2 * np.pi * (x - a * t_end))
        eg = l1_error(ug, exact, dx)
        er = l1_error(ur, exact, dx)
        og = f"{np.log2(prev['g']/eg):.2f}" if 'g' in prev else "  - "
        orr = f"{np.log2(prev['r']/er):.2f}" if 'r' in prev else "  - "
        print(f"{N:>5} {eg:.2e}({og}) {er:.2e}({orr})")
        prev = {'g': eg, 'r': er}
 
convergence_study()

The run prints this.

    N          Godunov              GRP
   40 1.63e-01(  - ) 1.31e-03(  - )
   80 8.76e-02(0.90) 2.69e-04(2.28)
  160 4.54e-02(0.95) 6.32e-05(2.09)
  320 2.31e-02(0.97) 1.55e-05(2.03)

Push the $\lambda\Delta t$ slider past 1.6 and explicit Euler either stays inside the band $[0,1]$ or breaks it. Moving the initial volume fraction and the equilibrium value shows whether Crank-Nicolson remains inside that band.

stiffness λΔt: 1.60α₀ (initial): 0.85α_eq (target): 0.10

Up to $\lambda\Delta t=3$ the cyan curve settles to equilibrium and the red one oscillates. The printed order 2.03 belongs to smooth advection. The slider is what happens when that half-step derivative is handed to a source integral.

At N=40 the error ratio is $1.63\times 10^{-1}/1.31\times 10^{-3}$ , about 124. The line about a hundred times smaller names this grid. At N=320 the ratio is $2.31\times 10^{-2}/1.55\times 10^{-5}$ , about 1490. An eightfold refinement, 40 to 320, cuts Godunov by 7.1 and GRP by 85. Eight squared is 64. The orders 2.28, 2.09, 2.03 overshoot 64 a little as they settle, while 0.90, 0.95, 0.97 climb toward 1.

CFL in the code is 0.4. The half-step correction has size $|a\,\sigma_i\,\Delta t/2|$ and the extrapolation $|\sigma_i\,\Delta x/2|$ . Their ratio is the CFL. The flux time pulls back 40 percent of the slope's push. At CFL 0.1 that share falls to 10 percent. The order proof stays. Less time derivative enters a single flux.

At N=80 the ratio is about 326, and at N=160 about 718. Doubling the grid steps the ratio through 124, 326, 718, 1490.

The stiff volume-fraction source of the Kapila model#

The volume-fraction equation in the Kapila five-equation model, a reduction of the BN seven-equation model, carries a stiff source.

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla\alpha_1 = \alpha_1 K\,\nabla\cdot\mathbf{u}

$\alpha_1$ is the volume fraction of phase 1. $K$ is the compaction/expansion coefficient taken from the Wood sound speed. How that speed bends with volume fraction is the same question as the Wood speed dropping inside a mixture. The source can push $\alpha_1$ outside $[0,1]$ . Explicit integration crosses the bound and diverges once the stiffness grows.

The GRP time derivative makes the interface state known at the new time $t^{n+1}$ . A Gauss-Green estimate then gives the cell average of $\nabla\cdot\mathbf{u}$ at that time. The source is integrated with Crank-Nicolson, semi-implicit, and a step-size limit can be written to keep the volume fraction inside the bound. The $\lambda\Delta t$ curves above are a model of that limit.

What felt shaky while reproducing it#

The scalar GRP came out second-order. The Kapila system is rougher. An unlimited central slope oscillates at a shock. The order showed up because the test was smooth. At a jump a limiter such as minmod is required, and the local order falls back to one.

The Wood sound speed is non-monotone in the volume fraction. That is why the Mach number wobbles inside the numerical smearing of an interface. The GRP time derivative follows compaction and expansion. It does not remove the non-monotonicity. The paper steps around that with robustness instead of a direct fix.

Production codes usually pair Strang splitting with Runge-Kutta, so a single-stage GRP does not drop in. Keeping the advantage means rewriting the flux in Lax-Wendroff form. The port still costs. The linearized acoustic step traces back to the gas-dynamics GRP of Ben-Artzi and Falcovitz.

Expecting the order while the source stays outside the flux#

The printed 2.03 is a number for a smooth slope. Turn a limiter on at a shock and that slot falls back to first order locally. Before deleting stages because the table says 2, the budget to rewrite the flux has to be there.

Leave the source outside the flux and some other time integrator inherits the same stiffness. Relaxation that carries two pressures and then collapses them to one either keeps the volume-fraction bound on that split side, or it cuts $\Delta t$ fine again. The printed 2.03 does not answer that.