[Paper Review] The Same Point Read 0.544 Apart Along Two Paths — Curl Involutions on a Staggered Grid

The same point read as two different numbers#

After 400 steps of pushing an interface field $\mathbf{b}$ around, I reconstructed the color function at the middle cell of the grid. Integrating $\mathbf{b}$ from the origin should give the color function back. That is what $\mathbf{b} = \nabla c$ means. Going in $x$ first and turning in $y$ gave one number. Going in $y$ first and turning in $x$ gave another. They differed by 0.544. The color function is a quantity that lives between 0 and 1.

Push the shear slider and the Rusanov dissipation slider together and the round interface starts to color on only one of the two grids.

The left panel is the staggered layout, the right one collocated. Watch which side stays black, and when the other lights up. That color is where the 0.544 path gap is accumulating.

This post reads Chiocchetti and Dumbser (2023), a staggered semi-implicit scheme for a first order hyperbolic model of viscous two-phase flow with surface tension, and reproduces one piece of it: the curl involution. Why that number has to be zero, why it is not zero under an ordinary variable layout, and what has to change to drive it down to round-off.

Promoting the interface field to an independent variable#

A diffuse interface model represents the boundary between two phases as a smooth transition in a color function $c$. Surface tension needs $\nabla c$ and its divergence. The trouble is the number of derivatives. Curvature is a second derivative of $c$, so when the interface is only a few cells wide, noise gets amplified straight through.

A first order hyperbolic reformulation lowers that by one. It promotes $\mathbf{b} = \nabla c$ to an independent variable and solves a transport equation for $\mathbf{b}$ itself. Apply $\nabla$ to both sides of $\partial_t c + \mathbf{v}\cdot\nabla c = 0$ and you get this.

Here $b_k$ is the $k$ component of the interface field, $v_j$ is the velocity, and $v_j b_j = \mathbf{v}\cdot\mathbf{b}$ is a single scalar. The whole system now contains first derivatives only. Eigenvalues are defined, and any hyperbolic solver applies. Break the constraint and a pressure that should have stayed quiet jumps with it. A nonconservative term that lifted the pressure 21% is the same kind of ledger.

There is a price. $\mathbf{b}$ is now stored independently of $c$. The two can drift apart.

The constraint is part of the equation, not part of the initial data#

If $\mathbf{b}$ is the gradient of any scalar, its curl vanishes identically.

That relation is called an involution. The transport equation preserves it on its own, because $\partial_t (\nabla \times \mathbf{b}) = -\nabla \times \nabla(v_j b_j) = \mathbf{0}$. At the continuous level, zero at $t=0$ means zero forever.

So it gets treated as an initial-condition matter. That is exactly the trap. Nothing guarantees that $\nabla_h \times \nabla_h$ vanishes for the discretized operators. Without that guarantee the constraint leaks a little every step, and after 400 steps the leak is the 0.544 from the opening.

What breaks when the curl is nonzero? Stokes answers that. The integral along a closed loop equals the area integral of the curl the loop encloses. A nonzero curl means $\mathbf{b}$ is not the gradient of any color function. The interface that corresponds to the stored $\mathbf{b}$ does not exist.

Where the discrete operators share the Schwarz identity#

The paper's prescription is simple: lay the discrete gradient and the discrete curl out so they use the same four points.

Put the color function at cell centers, $b_x$ on the $x$ faces $(i+\frac12, j)$, and $b_y$ on the $y$ faces $(i, j+\frac12)$. Compute the curl at the node $(i+\frac12, j+\frac12)$. For any scalar $\psi$ it then expands like this.

The two numerators are the same four terms written in a different order. The difference is identically zero. This is the Schwarz identity for mixed partials (mixed partials commute) carried intact to the discrete level. It is not an approximation. In floating point the same four numbers are added and subtracted, so only round-off survives.

Carrying that into the time integration takes one line. If every update has the form $-\Delta t\, \nabla_h \psi$, it disappears the moment a curl is applied. So $\nabla_h \times \mathbf{b}$ keeps its initial value regardless of how many steps are taken. If that value is zero, it stays zero.

Two discretizations under the same shear layer in Python#

One round interface on a $64^2$ periodic grid, pushed by a tanh shear layer. One run uses the staggered operator pair above; the other stores all of $\mathbf{b}$ at cell centers and uses a component-wise Rusanov flux. Same $\Delta t$, same velocity field.

import numpy as np
 
N, L = 64, 1.0
h = L / N
xc = (np.arange(N) + 0.5) * h
Xc, Yc = np.meshgrid(xc, xc, indexing="ij")
 
def shear_velocity(amp=1.0, delta=0.12):
    """Kelvin-Helmholtz style: tanh shear layer + sinusoidal perturbation."""
    return amp * np.tanh((Yc - 0.5) / delta), 0.35 * amp * np.sin(2 * np.pi * Xc)
 
def interface_potential(r0=0.22, w=0.05):
    r = np.sqrt((Xc - 0.5) ** 2 + (Yc - 0.5) ** 2)
    return 0.5 * (1.0 + np.tanh((r0 - r) / w))
 
def grad_stagger(phi):
    """b_x at x-faces (i+1/2, j), b_y at y-faces (i, j+1/2)."""
    return (np.roll(phi, -1, 0) - phi) / h, (np.roll(phi, -1, 1) - phi) / h
 
def curl_node(bx, by):
    """curl at nodes (i+1/2, j+1/2) -- the operator compatible with grad_stagger."""
    return (np.roll(by, -1, 0) - by) / h - (np.roll(bx, -1, 1) - bx) / h
 
def curl_centered(bx, by):
    """collocated curl: central differences of cell-centred components."""
    return ((np.roll(by, -1, 0) - np.roll(by, 1, 0))
            - (np.roll(bx, -1, 1) - np.roll(bx, 1, 1))) / (2 * h)
 
def psi_cell(bx, by, vx, vy, upwind):
    """the single scalar psi = v . b, collocated at cell centres."""
    if upwind:
        bxc = np.where(vx > 0, np.roll(bx, 1, 0), bx)
        byc = np.where(vy > 0, np.roll(by, 1, 1), by)
    else:
        bxc = 0.5 * (bx + np.roll(bx, 1, 0))
        byc = 0.5 * (by + np.roll(by, 1, 1))
    return vx * bxc + vy * byc
 
def update_compatible(bx, by, vx, vy, dt, upwind=True):
    """b <- b - dt * grad_h(psi): every added term is a discrete gradient."""
    psi = psi_cell(bx, by, vx, vy, upwind)
    return (bx - dt * (np.roll(psi, -1, 0) - psi) / h,
            by - dt * (np.roll(psi, -1, 1) - psi) / h)
 
def update_collocated(bx, by, vx, vy, dt, lam):
    """component-wise Rusanov: the dissipation is not the gradient of anything."""
    psi = vx * bx + vy * by
    fx = 0.5 * (psi + np.roll(psi, -1, 0)) - 0.5 * lam * (np.roll(bx, -1, 0) - bx)
    fy = 0.5 * (psi + np.roll(psi, -1, 1)) - 0.5 * lam * (np.roll(by, -1, 1) - by)
    return (bx - dt * (fx - np.roll(fx, 1, 0)) / h,
            by - dt * (fy - np.roll(fy, 1, 1)) / h)
 
def path_gap(bx, by):
    """integrate b along two paths to the same cell; the gap is the enclosed curl."""
    i1 = j1 = N // 2
    pa = bx[:i1, 0].sum() * h + by[i1, :j1].sum() * h   # x first, then y
    pb = by[0, :j1].sum() * h + bx[:i1, j1].sum() * h   # y first, then x
    return pa - pb
 
vx, vy = shear_velocity()
phi0 = interface_potential()
lam = float(np.max(np.abs(vx) + np.abs(vy)))
dt = 0.4 * h / lam
bs = list(grad_stagger(phi0))
bc = [np.gradient(phi0, h, axis=0, edge_order=2),
      np.gradient(phi0, h, axis=1, edge_order=2)]
 
print(f"N={N}  h={h:.4e}  dt={dt:.4e}  lambda={lam:.3f}")
print(f"{'t':>6}{'|b|max stag':>13}{'|curl|max stag':>16}"
      f"{'|b|max coll':>13}{'|curl|max coll':>16}")
for k in range(5):
    if k:
        for _ in range(100):
            bs = update_compatible(*bs, vx, vy, dt)
            bc = update_collocated(*bc, vx, vy, dt, lam)
    print(f"{k*100*dt:>6.3f}{np.max(np.hypot(*bs)):>13.3f}"
          f"{np.max(np.abs(curl_node(*bs))):>16.3e}"
          f"{np.max(np.hypot(*bc)):>13.3f}"
          f"{np.max(np.abs(curl_centered(*bc))):>16.3e}")
 
print(f"\npath-integral gap of phi at (0.5, 0.5):"
      f"  staggered {path_gap(*bs):+.3e}   collocated {path_gap(*bc):+.3e}")
 
bn = list(grad_stagger(phi0))
for _ in range(400):
    bn = update_compatible(*bn, vx, vy, dt, upwind=False)
print(f"staggered with NO upwinding: |b|max={np.max(np.hypot(*bn)):.3e} (unstable)"
      f"  |curl|max={np.max(np.abs(curl_node(*bn))):.3e}")

N=64  h=1.5625e-02  dt=4.6330e-03  lambda=1.349
     t  |b|max stag  |curl|max stag  |b|max coll  |curl|max coll
 0.000       10.167       0.000e+00        9.818       3.286e-04
 0.463       12.318       1.421e-12        6.538       5.545e+01
 0.927       18.063       1.748e-12        4.703       1.687e+02
 1.390        9.854       1.705e-12        4.124       2.269e+02
 1.853        5.746       2.132e-12        3.342       1.998e+02
 
path-integral gap of phi at (0.5, 0.5):  staggered -3.331e-16   collocated -5.442e-01
staggered with NO upwinding: |b|max=3.622e+03 (unstable)  |curl|max=1.091e-10

Raise the $\kappa$ slider and the two path bars split, while the shaded curl integral grows with that gap.

At zero the two bars stop on the same number. Any lift and the path gap matches that area integral.

The staggered curl sits at $10^{-12}$. That is double precision round-off multiplied by $1/h^2 \approx 4096$. The collocated one reaches 55 by $t = 0.46$ and 227 by $t = 1.39$. Since the largest $|\mathbf{b}|$ is about 4, the curl is more than fifty times bigger than the field carrying it. What is left is rotational noise, not an interface.

At $t=1.390$ the collocated curl peaks at 226.9. The largest $|\mathbf{b}|$ on that row is 4.124. Divide them and the ratio clears 50. Rotation outweighs the field. Integrate that curl along both paths to $(0.5, 0.5)$ and the gap is $-5.442\times 10^{-1}$. A color function lives between 0 and 1. The path gap is more than half of that range.

The time column is the step count times $dt=4.6330\times 10^{-3}$. One hundred steps land on 0.463. Four hundred land on 1.853, the last row. Between those rows the collocated curl climbs from $3.286\times 10^{-4}$ through 55.45, 168.7, and 226.9, then eases to 199.8. $|\mathbf{b}|$ has fallen to 3.342, so the absolute curl dips with it. The curl did not leave. Over the same 400 steps the staggered curl only moves from $1.421\times 10^{-12}$ to $2.132\times 10^{-12}$. Round-off times $1/h^2 \approx 4096$.

The no-upwind line is a different failure. $|\mathbf{b}|$ reaches $3.622\times 10^{3}$, about 360 times the initial 10.167, while the curl stays at $1.091\times 10^{-10}$. The printed $\lambda=1.349$ and $h=1.5625\times 10^{-2}$ set the size of one step of leak, whenever the update is not the gradient of a single scalar. Four hundred of those steps build the path gap of 0.544. The identity does not bound the amplitude. Dissipation is a separate knob.

Why more dissipation does not add curl#

The second run of update_compatible uses upwinding, which means dissipation is in there. The curl is still $10^{-12}$. This is the practical advantage of the scheme.

The reason is where the upwinding sits. The upwind decision is used only when the scalar $\psi = \mathbf{v}\cdot\mathbf{b}$ is gathered at cell centers. The update after that is still the single object $-\Delta t\, \nabla_h \psi$. However $\psi$ gets built, the curl of a gradient is zero.

The collocated layout is the opposite. The Rusanov dissipation $-\frac{\lambda}{2}(b_{x,i+1} - b_{x,i})$ attaches to $b_x$ in $x$ only and to $b_y$ in $y$ only. The two components receive different operators. That pair is not the gradient of any scalar, so applying a curl leaves $\frac{\lambda h}{2}(\partial_x \partial_y^2 b_y - \partial_y \partial_x^2 b_x)$ behind.

The restriction is on form, not on accuracy. As long as the entire update is the discrete gradient of one scalar, pile on limiters or high order reconstruction freely. Count what the grid can represent, then pick the scheme.

Zero curl does not stop a blow-up#

The last line of output makes that point. Turn upwinding off and keep only the central average, and $|\mathbf{b}|$ climbs to $3.6\times 10^3$. It started at 10, so that is a factor of 360. The run diverged. The curl, meanwhile, is $1.09\times 10^{-10}$, still at round-off.

Involution preservation and stability are separate properties. The staggered layout guarantees only that the field is still the gradient of some interface. It promises nothing about that interface staying finite. Dissipation is a separate ingredient. Conflating the two leads to a bad diagnosis: seeing a blow-up and suspecting the variable layout.

Where this changes numbers in surface tension#

The system in the paper builds the surface tension stress directly from $\mathbf{b}$. The stress tensor looks roughly like $\sigma(|\mathbf{b}|\,\mathbf{I} - \mathbf{b}\otimes\mathbf{b}/|\mathbf{b}|)$. Mix rotational content into $\mathbf{b}$ and that stress produces tension for an interface that does not exist.

The symptom shows up first in a static bubble. Spurious currents grow in a problem whose velocity should be exactly zero. The curl component enters the stress through the curvature, the stress makes velocity, and the velocity rotates $\mathbf{b}$ further. Once the feedback loop closes, the amplitude grows exponentially.

The paper solves pressure semi-implicitly on the same grid. Only the acoustic waves go implicit. The rest stays explicit. The staggered layout does two jobs here. It hands over the curl identity, and it makes the pressure system symmetric positive definite so a conjugate gradient solver applies. One layout buys both properties.

A zero curl that still blows up is not a grid bug#

$b_x$ and $b_y$ live on different faces, so every $\mathbf{v}\cdot\mathbf{b}$ pays for an average. That average is psi_cell. Averaging shaves high frequencies. On an interface three or four cells wide, ignoring that damping and then blaming the layout kills the interface width first. Face and node variables also want separate boundary stencils. On an unstructured mesh the compatible pair has to be redesigned, which is why the paper stays on Cartesian grids.

The decision splits in two. If the model carries an interface field or a distortion field as its own variable, the curl identity is not optional. It leaks every step, and the leak rate scales with $\lambda h$. If the model rebuilds $\mathbf{b}$ from $c$ every step, the leak never starts. Either way, ten lines of path_gap tell you whether a scalar exists for the stored vector field. The cost of putting pressure implicit on that same grid sits on the same ledger as a semi-implicit pressure step 1001 times larger.