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

Push the formal order to fifth, sixth, eighth — the contact still smears as time marches on. WENO or TENO, same fate. A polynomial drawing a jump does that. The 2023 paper by Takagi and colleagues steps off that path. High-order TENO keeps the smooth regions. Only discontinuous cells swap in a tanh (THINC).

Flip the reconstruction among Polynomial, THINC, and Hybrid and the same cell averages split into an S-curve, a vertical jump, or a mix.

Polynomial overshoots a little above and below the jump. THINC stands nearly vertical inside one cell. Hybrid keeps tanh only in the jump cell.

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.

What comes free is the decision of when to swap. 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.

The three reconstructions above are that choice. Same cell averages, different curves. Where TENO kills small eddies less is the subject of low-dissipation TENO. A contact bleeding across dozens of cells does not stop just because dissipation dropped.

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.

$d_i$ is not a free variable. 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. On the figure above, slide $\beta$ from 1.0 to 2.4 and the jump stands up while the overshoot readout moves with it.

At $\beta = 1.8$, $\beta/2 = 0.9$. A jump centered in the cell has $\alpha = 0$, so $d_i = 0$. At the right edge $X = 0.5$, $\tanh(\beta(X - d_i)) = \tanh(0.9) \approx 0.716$. Only about 72% of the neighbor jump is recovered inside one cell. At $\beta = 1.1$, $\tanh(0.55) \approx 0.50$, as gentle as van Leer. At $\beta = 2.4$, $\tanh(1.2) \approx 0.83$. The value 1.8 sits between a vertical wall and squaring.

$\alpha = (f_i - f_a)/f_d$. If the cell mean sits midway between the neighbors, $\alpha = 0$ and the jump center is 0. If the mean sticks to the right neighbor, $\alpha$ approaches 1, $(1 - T_2/T_1)$ inside the log shrinks, and $d_i$ shifts positive. If $|\alpha| \ge 1$ or $(f_{i+1} - f_i)(f_i - f_{i-1}) \le 0$, the code drops tanh and falls back to $\hat f = f_i$. Peak cells stay with the polynomial.

The initial top-hat lies between $0.35$ and $0.65$. Width $0.30$ on $N = 100$ is 30 cells. $dx = 0.01$ and $\Delta t = 0.45 \times dx = 0.0045$. Five revolutions are $5 / 0.0045 \approx 1111$ steps, at CFL 0.45. After that, minmod width is 23 cells and the peak is 0.84. Seven cells fall through 0.5 and the height loses 0.16. THINC-BVD keeps peak 1.00 and width 30. The two comment lines in the code are that arithmetic.

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#

Push revolutions and the red minmod curve and the cyan THINC curve run apart for the same time.

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.

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.

$\beta$ moves with the problem. 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.

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.

Freeze beta at 1.8?#

After five revolutions THINC-BVD shows peak 1.00 and width 30 cells. Minmod under the same conditions shows 0.84 and 23 cells. If smooth extrema get pressed into steps, that fixed value is too much. The choice that doubts one knob per wave is splitting THINC from MP5.

Drawing the jump with tanh, or sweeping up a smear afterward, is a different decision. If the reconstruction already holds the contact, a later sharpener is idle. If it does not, look at the detector before raising $\beta$ again.