[Paper Review] Pick the Reconstruction with Smaller Jumps — Deng (2018) MUSCL-THINC-BVD

In 2018 the Xiao group at Tokyo Tech measured a thickness nobody likes to publish. Compressible two-phase flow, 1000 steps. A jump that started one cell wide had spread across eight. WENO did not stop it. Dissipation can look small on a single step. Time still ruins the interface.

Their reconstruction is MUSCL-THINC-BVD. Every cell carries two candidates. The one with the smaller face jump stays.

Moving the β slider from 1.4 to 2.2 tightens the tanh inside a jump cell, and the Gaussian shoulders answer with teeth.

Watch whether cyan bars light only at the edges of the square wave. Also watch whether those bars stop multiplying near 1.6.

One-Page Summary#

• Authors / journal: Deng, Inaba, Xie, Shyue, Xiao. Journal of Computational Physics 371 (2018) 945–966.
• Problem: In the five-equation model for compressible two-phase flow, material interfaces smear as time runs.
• Proposal: Build MUSCL (smooth stretches) and THINC (jumps) in every cell, then keep the smaller face jump. That rule is BVD, boundary variation diminishing.
• What changes: No post-processed anti-diffusion and no artificial compression. The same rule hits volume fraction and the other conserved variables, so they stay consistent with each other.

Two Demands That Conflict#

Compressible two-phase numerics wants two things at once. Smooth regions need accuracy and little dissipation. Jumps, meaning interfaces and shocks, need monotonicity and a thin profile. One function rarely does both.

MUSCL is a straight line inside each cell. Monotonicity holds. The slope is only first-order. Each crossing of an interface shaves a little off. Over 1000 steps, one cell becomes eight. THINC fits a tanh inside the cell and holds a step in one or two cells. Drop that fit on a smooth region and fake stairs appear.

The paper starts with a decision. Do not blend the two into one function. Pick one per cell.

MUSCL — The Reliable Diffuser#

The baseline candidate is MUSCL with a minmod slope limiter.

$\bar q_i$ is the cell average. $\sigma_i$ is the slope at the center. minmod returns the smaller magnitude when the two differences share a sign, and zero when they do not.

Face values read like this.

No oscillations. In a jump, only the thickness grows, and it grows slowly.

THINC — A Monotone Function That Mimics Jumps#

THINC, tangent of hyperbola for interface capturing, fits a hyperbolic tangent inside the cell.

$\bar q_{\min}$ and $\bar q_{\max}$ are the minimum and the amplitude taken from the neighbors. $\theta$ is the jump direction. $\beta$ is the thickness. $\tilde x_i$ is the jump center, solved so the cell average stays put.

$\beta$ from 1.4 to 2.0 is the stable band. 1.6 is the usual pick. Larger values trap the jump in one cell, which looks clean, and they cut teeth into smooth regions. A record that bolts tanh onto a shock stencil, instead of picking a candidate, sits in the reconstruction that stacks TENO on THINC. There the two are blended by weights.

BVD — Measuring the Variation Between Two Candidates#

With two candidates, measure which face jump is smaller when they meet the neighbor reconstructions. Total boundary variation of cell $i$ is defined as follows.

$P$ is either MUSCL or THINC. Both neighbors stay on MUSCL. Only the middle cell switches.

The selection rule is this.

$C_i$ says where the jump sits inside the cell. $\delta$ is about $10^{-4}$.

The first gate multiplies the two one-sided differences and demands a strictly positive product. A zero factor fails it. A pure 0-to-1 step, one cell wide, is that shape: the cell before the jump has a right-hand difference and no left-hand one. MUSCL is returned before any boundary variation is built.

The second gate is location. $\delta$ sits near $10^{-4}$. A cell average inside that fraction of the neighbor amplitude is not eligible for THINC. Flat crests drop out.

Consistency on the Five-Equation Model#

The paper does not stop at a reconstruction contest. The five-equation model carries volume fraction $\alpha_1$, phasic densities $\alpha_1\rho_1$ and $\alpha_2\rho_2$, momentum $\rho u$, and total energy $E$. The same BVD rule is applied to all five.

Earlier methods sharpened only the volume fraction, then patched the rest so pressure would not ring at the interface. BVD forces the five variables to share one decision in each cell. A jump cell uses THINC for all five. A smooth cell uses MUSCL for all five. Consistency comes along. A separate anti-diffusion stage is not required.

The BVD Decision Tree in NumPy#

import numpy as np
 
def muscl_minmod_edges(q):
    """ Take (q[i-1], q[i], q[i+1]) and return MUSCL left/right edge values """
    qm, q0, qp = q
    a, b = q0 - qm, qp - q0
    if a * b <= 0.0:
        slope = 0.0
    else:
        slope = np.sign(a) * min(abs(a), abs(b))
    return q0 - 0.5 * slope, q0 + 0.5 * slope  # qL, qR
 
def thinc_jump_edges(q, beta=1.6, eps=1e-20):
    """ tanh-fit jump: left/right edge values of THINC """
    qm, q0, qp = q
    qmin = min(qm, qp)
    qmax = max(qm, qp) - qmin
    if qmax < 1e-12:
        return q0, q0
    theta = np.sign(qp - qm)
    C = (q0 - qmin + eps) / (qmax + eps)
    if C <= 1e-6 or C >= 1 - 1e-6:
        return q0, q0
    B = np.exp(theta * beta * (2.0 * C - 1.0))
    A = (B / np.cosh(beta) - 1.0) / np.tanh(beta)
    qR = qmin + 0.5 * qmax * (1.0 + theta * A)
    num = 1.0 + theta * A * np.tanh(beta) + theta * np.tanh(beta)
    den = 1.0 + A * np.tanh(beta)
    qL = qmin + 0.5 * qmax * (num / den)
    return qL, qR
 
def bvd_decide(q_window, beta=1.6, delta=1e-4):
    """
    q_window: length-5 array [q[i-2], q[i-1], q[i], q[i+1], q[i+2]]
    returns: (qL_i, qR_i, picked_thinc)
    """
    qm2, qm1, q0, qp1, qp2 = q_window
    monotone = (qp1 - q0) * (q0 - qm1) > 0.0
    qmin = min(qm1, qp1)
    qmax = max(qm1, qp1) - qmin
    C = 0.5 if qmax < 1e-12 else (q0 - qmin) / qmax
    eligible = monotone and (delta < C < 1.0 - delta)
 
    qL_M, qR_M = muscl_minmod_edges([qm1, q0, qp1])
    if not eligible:
        return qL_M, qR_M, False
 
    qL_T, qR_T = thinc_jump_edges([qm1, q0, qp1], beta=beta)
    # TBV with both neighbors held at MUSCL
    _, qR_left = muscl_minmod_edges([qm2, qm1, q0])
    qL_right, _ = muscl_minmod_edges([q0, qp1, qp2])
    tbv_M = abs(qR_left - qL_M) + abs(qR_M - qL_right)
    tbv_T = abs(qR_left - qL_T) + abs(qR_T - qL_right)
    if tbv_T < tbv_M:
        return qL_T, qR_T, True
    return qL_M, qR_M, False
 
# Test: a smooth Gaussian + a square jump
N = 80
x = (np.arange(N) + 0.5) / N
q = np.where((x >= 0.6) & (x <= 0.8), 1.0, 0.0)
q += 0.5 * np.exp(-((x - 0.25) / 0.04) ** 2)
 
picks = 0
for i in range(2, N - 2):
    _, _, used_thinc = bvd_decide(q[i-2:i+3])
    picks += int(used_thinc)
print(f"THINC pick rate: {picks}/{N-4}")
# Expected: a few cells around the jump only

The snippet comment says 4 of 76 cells adopt THINC. The gate above it does not describe that profile. The coded square wave is 0 or 1, so one difference is zero and the product test fails. Those four cells need a same-sign ramp across two or three cells. The eight-cell smear after 1000 steps is already such a ramp. At one-cell thickness this branch never calls THINC.

On a straight line the neighbor differences match, so the minmod slope equals them and MUSCL boundary variation is zero. THINC's tanh opens a face gap, and the inequality keeps MUSCL. A Gaussian wider than a cell is locally that line, which is why smooth stretches refuse the swap.

A smooth Gaussian is locally close to a line, so the swap is refused. A square wave that is still one cell thick never opens the gate either.

Who Wins Each Cell — See It Directly#

The advection above pushes three reconstructions together and paints a bar on every cell where BVD picks THINC. Raising $\beta$ from 1.4 to 2.2 narrows only those cells. Too large, and teeth show up on the Gaussian shoulders. The bar count stops growing near 1.6, which is why that value stuck. Cyan appears at the edges only after the square wave has smeared a little. While the thickness is one cell, the bars stay empty.

How Far You Can Trust It#

The observed order stays below two. One candidate is MUSCL, so smooth accuracy stops where MUSCL stops. Swapping in a high-order candidate such as WENO can reach fifth order. This paper chose the plain version.

TBV assumes the neighbors are MUSCL. A case where both sides want THINC together is ignored inside one sweep. Section 4 of the paper still looks acceptable under that shortcut. A more general form was left for later work.

In multiple dimensions, $\beta$ splits by face normal. The interface normal is estimated with Young's algorithm. At corners and kinks the normal is wrong, and the jump leans to one side.

Time integration is separate. The paper uses SSPRK3. The BVD decision has to be rebuilt at every Runge–Kutta stage. If the decision flickers between stages, the interface jitters. Later work went back to that jitter.

Where the picker refuses a one-cell jump#

While the interface is still one cell thick, this branch never calls THINC. The product of differences turns positive only after the thickness has spread over two or three cells. Keeping a one-cell interface and being allowed to choose only after it has already spread are two jobs that miss each other.

Whether that miss leaks into the time-step limit is the next question. A thin interface still has an advection ceiling. The interface advection that stayed at CFL 0.05 after the capillary term was removed asks where that ceiling sits. If BVD really holds the jump in one cell, how large a step can move that cell?