[Paper Review] Drawing Shocks with tanh — Implementing TENO-THINC Reconstruction

Push the formal order to fifth, sixth, eighth — the contact discontinuity still smears as time marches on. WENO or TENO, it makes no difference. As long as you draw a jump with a polynomial, that fate is sealed. The 2023 paper by Takagi and colleagues sidesteps it. High-order TENO handles the smooth regions; only the discontinuous cells get swapped for a tanh curve (THINC). This post implements that core idea at the single-cell level and spins a top-hat once around the grid to compare the fate of polynomials and tanh directly.

The paper on one page#

• Title: High-Order Low-Dissipation Shock-Resolving TENO-THINC Schemes for Hyperbolic Conservation Laws
• Authors: Shinichi Takagi, Hiro Wakimura, Lin Fu, Feng Xiao
• Venue: Communications in Computational Physics, 2023
• One line: Couple even-point TENO (six- and eight-point) with THINC reconstruction to get low dissipation in smooth regions and sub-cell resolution at discontinuities at the same time.

The trick lies in deciding when to swap, for free. TENO already tags each candidate stencil as smooth or not via $\delta_k \in \{0,1\}$. Reusing those flags tells you whether a cell hides a discontinuity, and THINC fires only there. The only new parameter is THINC's slope $\beta$.

Why even high-order polynomials smear the contact#

Polynomial reconstruction draws a smooth curve by nature. When a curve must pass through cell-average points across a near-vertical jump, it faces a choice: allow oscillation (overshoot), or kill the slope and smear.

TENO and WENO choose to suppress oscillation. They discard stencils that cross the discontinuity, enforcing the ENO property. That buys stability, but a little numerical dissipation accumulates every step. A contact discontinuity (a surface where velocity and pressure match but density jumps) does not self-steepen through a pressure mechanism. So after long advection its width bleeds across dozens of cells. Raising the order barely reduces this smearing.

See it for yourself below. Same cell averages, only the reconstruction changes.

Polynomial resolves the jump as a smooth S-curve and slightly overshoots above and below. Switch to THINC and the same data stands nearly vertical inside a single cell. Hybrid (BVD) uses tanh only in the jump cell and keeps polynomials elsewhere.

THINC — tanh draws a sub-cell jump#

THINC (Tangent of Hyperbola for INterface Capturing) replaces the polynomial with a monotone tanh. The reconstruction over the cell-local coordinate $X \in [-0.5, 0.5]$ in cell $i$ is

where $f_a = (f_{i+1}+f_{i-1})/2$ is the neighbor mean, $f_d = (f_{i+1}-f_{i-1})/2$ is half the jump, $\beta$ is the tanh slope, and $d_i$ is the jump-center location.

The key is that $d_i$ is not free. The cell-average conservation constraint (the reconstruction must integrate back to $f_i$) pins it down. The solution is analytic:

with $T_1 = \tanh(\beta/2)$, $T_2 = \tanh(\alpha\beta/2)$, $\alpha = (f_i - f_a)/f_d$. Here $\alpha$ is the normalized position of the cell mean between its neighbors. If the cell is non-monotone ($(f_{i+1}-f_i)(f_i-f_{i-1}) \le 0$) the tanh is meaningless, so it falls back to first order ($\hat f_{i+1/2} = f_i$).

β: the knob that holds sharpness#

A single $\beta$ sets THINC's character: small is gentle, large is steep. The paper records a useful correspondence.

• $\beta \approx 1.1$ : TVD spectrum at the van Leer limiter level
• $\beta \approx 1.3$ : Superbee level (over-compressive tendency)
• $\beta = 1.6 \sim 2.0$ : sharp shock-capturing regime

The paper adopts $\beta = 1.8$. Push it too high and even smooth curves get squared into steps, producing artificial flattening. In the simulation above, slide $\beta$ from 1.0 up to 2.4 and watch how the jump stands up and how the overshoot readout changes — the balance becomes tangible.

TENO finds it, THINC fills it — the BVD hybrid#

Apply THINC everywhere and even smooth waves break. So it must fire only in discontinuous cells. The paper builds a per-cell indicator $\zeta_i$ from the $\delta_k$ that TENO weighting already computed.

$\zeta_i = 1$ signals a discontinuity on the right, $-1$ on the left. When the $\zeta$ signs of the neighbors face each other across cell $i$, that cell is judged to truly hold a discontinuity and gets swapped for THINC.

The practical BVD (Boundary Variation Diminishing) principle is simpler still. Compute both reconstructions, then keep whichever shrinks the difference between left and right reconstructed values at the cell interface (the boundary variation). The implementation below follows this BVD test.

Python — spinning a top-hat once around#

Advect a top-hat under linear advection $u_t + u_x = 0$ with periodic boundaries. Compare minmod TVD against THINC-BVD on the same grid for the same time.

import numpy as np
 
def minmod(p, q):
    return np.where(p * q <= 0, 0.0, np.where(np.abs(p) < np.abs(q), p, q))
 
def thinc_right_edge(u, beta=1.8):
    # Eqs. 2.19-2.20: left reconstructed value at the right face (i+1/2) of cell i
    um, up = np.roll(u, 1), np.roll(u, -1)
    fa, fd = 0.5 * (up + um), 0.5 * (up - um)
    mono = ((up - u) * (u - um) > 0) & (np.abs(fd) > 1e-12)
    alpha = np.where(mono, (u - fa) / np.where(mono, fd, 1.0), 0.0)
    mono = mono & (np.abs(alpha) < 1.0)
    T1, T2 = np.tanh(beta / 2), np.tanh(alpha * beta / 2)
    d = np.where(mono, (1 / (2 * beta)) * np.log((1 - T2 / T1) / (1 + T2 / T1)), 0.0)
    vR = fa + fd * np.tanh(beta * (0.5 - d))
    return np.where(mono, vR, u)  # non-monotone -> first order
 
def reconstruct(u, beta, use_thinc):
    um, up = np.roll(u, 1), np.roll(u, -1)
    s = minmod(u - um, up - u)
    poly_vR = u + 0.5 * s
    if not use_thinc:
        return poly_vR
    poly_vL = u - 0.5 * s
    t = thinc_right_edge(u, beta)
    # BVD: keep the reconstruction that shrinks the right-face jump
    jump_poly = np.abs(poly_vR - np.roll(poly_vL, -1))
    jump_thinc = np.abs(t - np.roll(poly_vL, -1))
    return np.where(jump_thinc < jump_poly, t, poly_vR)
 
def rhs(u, dx, beta, use_thinc):
    vR = reconstruct(u, beta, use_thinc)        # a=1, so F_{i+1/2} = vR
    return -(vR - np.roll(vR, 1)) / dx
 
def advect_tophat(N=100, revolutions=5.0, beta=1.8, use_thinc=True):
    dx, x = 1.0 / N, (np.arange(N) + 0.5) / N
    dt = 0.45 * dx
    u = ((x > 0.35) & (x < 0.65)).astype(float)   # initial top-hat
    for _ in range(int(revolutions / dt)):
        u1 = u + dt * rhs(u, dx, beta, use_thinc)            # SSP-RK2
        u = 0.5 * u + 0.5 * (u1 + dt * rhs(u1, dx, beta, use_thinc))
    return x, u
 
x, u_tvd = advect_tophat(use_thinc=False)
_, u_thinc = advect_tophat(use_thinc=True)
print("TVD   peak=%.3f  width(>0.5)=%d" % (u_tvd.max(), (u_tvd > 0.5).sum()))
print("THINC peak=%.3f  width(>0.5)=%d" % (u_thinc.max(), (u_thinc > 0.5).sum()))
# Example output:
# TVD   peak=0.84  width(>0.5)=23
# THINC peak=1.00  width(>0.5)=30

After five revolutions, minmod TVD sags to a peak of 0.84 and narrows. THINC-BVD holds peak 1.00 and keeps nearly its initial width (30 cells). Even a higher-order TENO cannot fully stop this smearing, because of the polynomial nature itself — that is the paper's starting point.

Reconstruction and long-time advection, hands-on#

Play with the simulation below. It shows in real time how the two methods diverge as the same top-hat circles the grid.

As revolutions grows, the red minmod curve collapses at both ends and spreads into a trapezoid. The cyan THINC-BVD clings to the gray exact solution and keeps its rectangle. Drop $\beta$ to 1.1 and THINC softens like TVD; near 2.0 it is sharpest.

What felt shaky in the reproduction#

The promise is clean, but three things tripped me up while coding it.

First, the BVD decision is sensitive to detector quality. Mistake a smooth extremum for a discontinuity and that peak flattens unnaturally. The simple BVD above guards against it with the monotonicity check, but multidimensional systems need characteristic decomposition.

Second, $\beta$ is problem-dependent. The value $\beta = 1.8$ is tuned to 1D benchmarks. In strong compressible turbulence over-compression can create artificial steps, so an adaptive $\beta$ is worth considering.

Third, the cost claim depends on context. "THINC is nearly free because discontinuous cells are few" is true when you look only at reconstruction. But a branch-heavy BVD test can trigger warp divergence on GPUs. Porting it to an unstructured FV code like OpenFOAM means folding per-cell detection into the face-flux loop, which is not as clean as a 1D structured grid.

What to read next#

If you want THINC's roots, read Xiao's original VOF THINC; for the BVD principle itself, Deng's P$n$T$m$-BVD. Adaptive $\beta$ and multiphase extensions are carried on by the THINC-with-PE line of Shyue & Xiao. Today's one line: do not draw a jump with a polynomial — draw it with a jump.