Filling an Interface with a Fluid That Isn't There — The Ghost Fluid Method

Into the cell where water meets air, you pour fake air that isn't actually there. This ghost fluid (a fictitious fluid extrapolated across the interface) exists nowhere on the grid. Yet that very ghost is what keeps the solver stable when a shock crosses the interface. What gets copied and what gets extrapolated are different jobs. The original GFM rings where that split is wrong.

Raise the density jump and the staircase gets steeper, while pressure still walks straight through the red dashed line. "show ghost" toggles the extension past the interface.

With the dashed lines on, each fluid carries its own entropy across the boundary. Slide the interface and the pressure curve stays one piece. Density is the quantity that breaks.

When two materials collide in one cell#

Gas and liquid, or explosive and metal — two materials with wildly different properties meet on the same grid. Use averaged properties and the interface smears numerically. Across a water-air interface with a 1000:1 density ratio, that smearing quickly turns into nonphysical oscillations.

GFM takes a different route. It does not blur the interface. It extends each fluid as a fake fluid past its own region, then runs a single-material solver twice, as if only one material existed. The interface is tracked as the zero of a level set (a signed-distance function that splits the domain) $\phi$.

Here $\phi$ is the signed-distance function, and its zero set $\phi=0$ is the interface.

What breaks and what continues across the interface#

To fill the ghost correctly, you must know exactly what is continuous across the interface. Pressure $p$ and normal velocity $u_n$ are continuous there. Entropy $s$ and tangential velocity $u_t$ are not.

The brackets $[\cdot]$ denote the jump across the interface. These four lines are the whole of GFM. Copy the continuous quantities, extrapolate the discontinuous ones. The figure above is that split: pressure stays smooth, density is the step.

Averaging both phases inside one cell fails a different way. The six-equation two-fluid model picks up complex eigenvalues once the density ratio is large, and refining the grid makes the blow-up arrive sooner. GFM refuses that average. A cell holds one material.

Treating water like a gas — the stiffened EOS#

Water is nearly incompressible, so the ideal-gas equation of state cannot handle it. To bind both materials into the same Mie–Grüneisen form, GFM uses the stiffened-gas EOS.

$\rho$ is density, $e$ is internal energy per unit mass, $\gamma$ is the specific-heat ratio, and $p_\infty$ is a constant that encodes material stiffness. An ideal gas has $p_\infty = 0$; water is roughly $\gamma \approx 4.4$, $p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$. The sound speed follows as:

Thanks to $p_\infty$, water's sound speed comes out far larger than air's at the same pressure. Unify the EOS into one line and both fluids fit into the same solver.

Two touches that fill a ghost cell#

Here is the core procedure. When solving fluid A's region, fill the cells occupied by B with A's ghost. The touch splits in two.

1. Copy: the continuous quantities — pressure and normal velocity — are taken directly from the real B cells across the interface.
2. Extrapolate: the discontinuous quantities — entropy (hence density) and tangential velocity — are extrapolated isentropically from the nearest real A cell across the interface.

Isentropic extrapolation preserves the entropy $s = p / \rho^\gamma$ and re-solves the density.

$s_A$ is the real A fluid's entropy and $p_\text{copy}$ is the copied pressure. The extrapolation is usually carried out with a fast-marching method that propagates values quickly along the interface normal. Once a single sweep completes A's ghost field, you forget B exists and solve A's single-material Euler equations. Repeat the same for B.

The initial data already fix the ghost densities. Left state: density 1, pressure 1, $\gamma = 1.4$, entropy 1. Copy the right-hand pressure 0.1 and the ghost density is $0.1^{1/1.4}$, about 0.193, not the real 0.125. The single-material solver sees a ratio near 5.2, not 8.

The right gas has $\gamma = 1.67$. Its entropy is $0.1 / 0.125^{1.67}$, about 3.22. Copy pressure 1 onto that entropy and the ghost density is about 0.496, half the real left state. Pressure crosses. Entropy stays with the fluid that owns the ghost.

With $p_\infty = 0$, sound speed is the square root of $\gamma p / \rho$. Left is about 1.183, right about 1.156, so the left state sets the CFL. Spacing on 400 points is $1/399$, about 0.00251. CFL 0.4 gives a step near 0.00085, and $0.14$ divided by that step is about 165 steps. A sound wave needs about 0.85 to cross the tube. The end time 0.14 is under a fifth of that. A density ratio of 8 is mild next to water-air at 1000, and 0.193 against 0.496 are already off the real states.

Where the original GFM rings#

The original GFM runs into trouble when one side of the interface is very stiff (like water, with large $p_\infty$). When a strong shock or detonation crosses the interface, the ghost state built by plain copy and extrapolation disagrees with the true Riemann solution. That mismatch shows up as pressure oscillations near the interface, and the oscillations often lead to divergence.

In the simulation below, send a pressure pulse through the stiff fluid (pink region). The higher you push the stiffness ratio, the larger the ripples just past the interface.

Push the stiffness ratio to 9 and a high-frequency oscillation appears clearly just past the interface. It does not go away when you refine the grid, because the cause is the boundary condition, not the grid.

The modified GFM — set values from an interface Riemann problem#

The fix is to solve for the ghost state rather than guess it. The modified GFM (MGFM or IGFM) takes the real states on each side of the interface as the left and right states of a Riemann problem and solves it.

$W_A, W_B$ are the real states on each side, $\mathcal{R}$ is the Riemann solver, and $p^*, u^*$ are the contact-discontinuity pressure and velocity. Assign the same $p^*, u^*$ to both ghosts and the pressure–velocity continuity conditions are exactly satisfied. Turn on "use modified GFM" above and the oscillation disappears at the same stiffness.

Python — a gas-gas GFM shock tube#

Solve a shock tube where two ideal gases ($\gamma_A = 1.4$, $\gamma_B = 1.67$) meet, using isentropic-extrapolation GFM. It uses an HLL flux and, each step, two ghost sweeps merged by the level set.

import numpy as np
 
GAMMA = (1.4, 1.67)   # specific-heat ratios of left (air) and right (monatomic gas)
 
def primitive(U, g):
    rho = U[0]; u = U[1] / rho
    e = U[2] / rho - 0.5 * u * u
    p = (g - 1.0) * rho * e          # stiffened EOS, p_inf = 0 (ideal gas)
    return rho, u, p
 
def conserved(rho, u, p, g):
    E = p / (g - 1.0) + 0.5 * rho * u * u
    return np.array([rho, rho * u, E])
 
def hll(UL, UR, g):
    rL, uL, pL = primitive(UL, g); rR, uR, pR = primitive(UR, g)
    aL = np.sqrt(g * pL / rL); aR = np.sqrt(g * pR / rR)
    sL = min(uL - aL, uR - aR); sR = max(uL + aL, uR + aR)
    FL = np.array([rL * uL, rL * uL * uL + pL, uL * (UL[2] + pL)])
    FR = np.array([rR * uR, rR * uR * uR + pR, uR * (UR[2] + pR)])
    if sL >= 0: return FL
    if sR <= 0: return FR
    return (sR * FL - sL * FR + sL * sR * (UR - UL)) / (sR - sL)
 
def extrapolate_ghost(rho, u, p, mat, side):
    # fill the whole domain with `side` fluid: copy p,u; extrapolate density isentropically
    rg, ug, pg = rho.copy(), u.copy(), p.copy()
    g = GAMMA[side]
    idx = np.where(mat == side)[0]
    lo, hi = idx[0], idx[-1]
    ref = lo if side == 1 else hi      # real cell closest to the interface
    s_ref = p[ref] / rho[ref] ** g     # entropy to preserve
    for i in range(len(rho)):
        if mat[i] != side:
            j = min(max(i, lo), hi)            # nearest real cell
            ug[i] = u[j]; pg[i] = p[j]         # copy
            rg[i] = (pg[i] / s_ref) ** (1.0 / g)   # isentropic extrapolation
    return rg, ug, pg
 
def run_ghost_fluid_tube(N=400, t_end=0.14, cfl=0.4):
    x = np.linspace(0, 1, N); dx = x[1] - x[0]
    x0 = 0.5                                   # initial interface position
    mat = np.where(x < x0, 0, 1)
    rho = np.where(x < x0, 1.0, 0.125)
    u   = np.zeros(N)
    p   = np.where(x < x0, 1.0, 0.1)
    t = 0.0
    while t < t_end:
        a = np.sqrt(np.where(mat == 0, GAMMA[0], GAMMA[1]) * p / rho)
        dt = min(cfl * dx / np.max(np.abs(u) + a), t_end - t)
        new = {}
        for side in (0, 1):                    # two ghost sweeps
            rg, ug, pg = extrapolate_ghost(rho, u, p, mat, side)
            g = GAMMA[side]
            U = np.array([conserved(rg[i], ug[i], pg[i], g) for i in range(N)]).T
            F = np.zeros((3, N + 1))
            for i in range(1, N):
                F[:, i] = hll(U[:, i - 1], U[:, i], g)
            F[:, 0] = F[:, 1]; F[:, N] = F[:, N - 1]
            new[side] = U - dt / dx * (F[:, 1:] - F[:, :-1])
        U = np.where(mat[None, :] == 0, new[0], new[1])   # level-set merge
        for i in range(N):
            rho[i], u[i], p[i] = primitive(U[:, i], GAMMA[mat[i]])
        x0 += dt * u[np.argmin(np.abs(x - x0))]           # move the interface
        mat = np.where(x < x0, 0, 1)
        t += dt
    return x, rho, u, p
 
if __name__ == "__main__":
    x, rho, u, p = run_ghost_fluid_tube()
    print(f"t=0.14:  max p = {p.max():.3f}, "
          f"density near contact {rho[180]:.3f} -> {rho[220]:.3f}")

The output shows the density jump surviving across the interface while the pressure stays smooth. Two single-material sweeps are the backbone of GFM.

Traps to avoid when you write the code#

Extrapolate along the grid axes and entropy mixes on a slanted interface. The sweep direction is the interface normal, nothing else. A high stiffness ratio, water against air, makes the original copy-and-extrapolate step ring in the pressure. In that case the ghost is filled from the interface Riemann state from the start. Skip the level-set reinitialization and the zero contour can linger while the distance values drift, so the interface speed goes wrong.

When a gas-gas extrapolation is dropped onto water#

Two ideal gases with $p_\infty = 0$ ($\gamma = 1.4$ and $1.67$) keep a density step at the contact under plain isentropic extrapolation. Water is $\gamma \approx 4.4$, $p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$. Add that constant inside the sound-speed numerator at the same pressure, and the ripples seen when the demo stiffness ratio hits 9 are a smaller cousin of the same mismatch. Refining the grid leaves them in the boundary condition.

Decide whether the ghost pressure and velocity are a copy, or the interface Riemann state $p^*, u^*$, before counting the density jump. A ghost density that falls through zero is a separate limiter problem. The positivity-preserving flux limiter is the next place to look: mass can be right while density in a cell is not.