MacCormack Advection on Staggered Grids in GLSL

This article provides a technical overview of implementing the MacCormack advection method on staggered (Marker-and-Cell) grids using GLSL compute shaders. It details the coordinate handling required for staggered velocity components, the multi-step forward-backward advection calculations, the error-correction formulation, and the essential clamping mechanisms needed to prevent numerical instability and overshoot artifacts in real-time fluid simulations.

The Staggered Grid Layout

In a staggered Marker-and-Cell (MAC) grid configuration, physical quantities are defined at distinct spatial offsets within each grid cell:

  • Cell Centers \((i, j)\): Store scalar quantities such as pressure, density, and temperature.
  • Horizontal Faces \((i \pm \frac{1}{2}, j)\): Store the horizontal velocity component \(u\).
  • Vertical Faces \((i, j \pm \frac{1}{2})\): Store the vertical velocity component \(v\).

Because \(u\) and \(v\) do not share spatial locations, advecting a velocity component requires reconstructing the full velocity vector at that component's specific face location before performing characteristic tracing.

The MacCormack Advection Scheme

Standard semi-Lagrangian advection suffers from severe numerical dissipation (blurring of vortices and fine details). The MacCormack method mitigates this by evaluating a forward advection step, reversing it with a backward advection step to estimate the numerical error, and applying that error as a correction.

The scheme proceeds in four discrete mathematical steps:

  1. Forward Step: Compute an intermediate candidate field \(\phi^*\) by standard semi-Lagrangian back-tracing: \[\phi^* = \text{Advect}(\phi^n, \mathbf{u}, \Delta t)\]
  2. Backward Step: Advect the intermediate field forward in time (reversing the velocity vector): \[\phi^{**} = \text{Advect}(\phi^*, -\mathbf{u}, \Delta t)\]
  3. Error Estimation and Correction: Compute the difference between the original and reconstructed fields to cancel first-order diffusion: \[\phi^{n+1}_{\text{raw}} = \phi^* + \frac{1}{2}(\phi^n - \phi^{**})\]
  4. BPS Clamping: Prevent dispersion artifacts and non-physical oscillations by clamping the corrected value to the local bounding box of the forward-advection lookup: \[\phi^{n+1} = \text{clamp}(\phi^{n+1}_{\text{raw}}, \phi_{\min}, \phi_{\max})\]

GLSL Compute Shader Evaluation Pipeline

Evaluating MacCormack advection in compute shaders typically requires two dispatch passes due to data dependencies between \(\phi^*\) and \(\phi^{**}\).

1. Velocity Interpolation at Staggered Faces

To advect the horizontal component \(u\) located at index \((x, y)\), you must sample the co-located full velocity vector \(\mathbf{v}_u = (u, \bar{v})\):

// Reconstructing full velocity at the u-component face
float get_v_at_u_face(sampler2D vTex, vec2 uv, vec2 texelSize) {
    // Average four neighboring v-values surrounding the u-face
    float v0 = texture(vTex, uv + vec2(-0.5, -0.5) * texelSize).r;
    float v1 = texture(vTex, uv + vec2( 0.5, -0.5) * texelSize).r;
    float v2 = texture(vTex, uv + vec2(-0.5,  0.5) * texelSize).r;
    float v3 = texture(vTex, uv + vec2( 0.5,  0.5) * texelSize).r;
    return 0.25 * (v0 + v1 + v2 + v3);
}

An analogous interpolation is performed when advecting the vertical component \(v\), averaging four neighboring \(u\) values.

2. Pass 1: Forward Advection (\(\phi^*\))

A compute shader samples the current field \(\phi^n\) along the backward characteristic \(\mathbf{x} - \mathbf{u}(\mathbf{x})\Delta t\) (often integrated using second-order Runge-Kutta). The result is written to an intermediate texture texPhiHat.

layout (local_size_x = 16, local_size_y = 16) in;

layout (binding = 0) uniform sampler2D uTex;
layout (binding = 1) uniform sampler2D vTex;
layout (binding = 2) uniform sampler2D phiSrc;
layout (binding = 3, r32f) writeonly uniform image2D phiHatOut;

uniform float dt;
uniform vec2 texelSize;

void main() {
    ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
    vec2 uv = (vec2(coord) + 0.5) * texelSize;

    // Velocity reconstruction at current location
    vec2 vel = vec2(texture(uTex, uv).r, get_v_at_u_face(vTex, uv, texelSize));

    // Midpoint Runge-Kutta (RK2) back-trace
    vec2 midUV = uv - 0.5 * dt * vel * texelSize;
    vec2 midVel = vec2(texture(uTex, midUV).r, get_v_at_u_face(vTex, midUV, texelSize));
    vec2 traceUV = uv - dt * midVel * texelSize;

    float phiStar = texture(phiSrc, traceUV).r;
    imageStore(phiHatOut, coord, vec4(phiStar, 0.0, 0.0, 0.0));
}

3. Pass 2: Backward Evaluation, Correction, and Clamping

The second compute shader carries out the backward trace, reads from the original and intermediate buffers, and computes the clamped correction.

layout (local_size_x = 16, local_size_y = 16) in;

layout (binding = 0) uniform sampler2D uTex;
layout (binding = 1) uniform sampler2D vTex;
layout (binding = 2) uniform sampler2D phiSrc;     // phi^n
layout (binding = 3) uniform sampler2D phiHatTex;  // phi^*
layout (binding = 4, r32f) writeonly uniform image2D phiFinalOut;

uniform float dt;
uniform vec2 texelSize;

void main() {
    ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
    vec2 uv = (vec2(coord) + 0.5) * texelSize;

    vec2 vel = vec2(texture(uTex, uv).r, get_v_at_u_face(vTex, uv, texelSize));

    // Forward trace UV (for clamping bounds)
    vec2 fwdUV = uv - dt * vel * texelSize;
    
    // Backward trace UV from current location (reversing advection)
    vec2 bwdUV = uv + dt * vel * texelSize;

    float phiStar = texture(phiHatTex, uv).r;
    float phiStarStar = texture(phiHatTex, bwdUV).r;
    float phiOriginal = texture(phiSrc, uv).r;

    // MacCormack error correction
    float phiCorrected = phiStar + 0.5 * (phiOriginal - phiStarStar);

    // Compute neighborhood min/max at the forward-trace point
    vec2 base = fwdUV / texelSize - 0.5;
    ivec2 iBase = ivec2(floor(base));

    float c00 = texelFetch(phiSrc, iBase + ivec2(0, 0), 0).r;
    float c10 = texelFetch(phiSrc, iBase + ivec2(1, 0), 0).r;
    float c01 = texelFetch(phiSrc, iBase + ivec2(0, 1), 0).r;
    float c11 = texelFetch(phiSrc, iBase + ivec2(1, 1), 0).r;

    float phiMin = min(min(c00, c10), min(c01, c11));
    float phiMax = max(max(c00, c10), max(c01, c11));

    // Clamping to suppress overshoot/undershoot
    float phiClamped = clamp(phiCorrected, phiMin, phiMax);

    imageStore(phiFinalOut, coord, vec4(phiClamped, 0.0, 0.0, 0.0));
}

Boundary Considerations on Staggered Compute Domains

When dispatching compute workloads over staggered domains:

  1. Resolution Mismatch: A grid with \(N \times N\) scalar cells contains \((N+1) \times N\) cells for \(u\) and \(N \times (N+1)\) cells for \(v\). Compute dispatch dimensions must account for the additional boundary faces.
  2. Boundary Wall Checks: During the forward and backward characteristic traces, trajectory points must not exit the fluid domain. UV lookups must be clamped to valid interior bounds to avoid sampling uninitialized or invalid domain edges.