NEON FLUID

# Overview

The screen is a tank of dark fluid. The cursor pushes the fluid and puts neon dye into it. The dye moves with the flow, rolls into vortices, and fades in approximately one second. A tone curve, a small light effect, and faint scanlines give the dye its neon look.

The effect is a stable-fluids solver (Jos Stam, 1999) in raw WebGL2, with no library. The velocity, the pressure, and the dye are half-float textures. Each step of the solver is a fragment shader that draws one full-screen quad into a texture. The module makes nine programs from one vertex shader and the nine fragment shaders of the FS object: clear, splat, advect, divergence, curl, vorticity, pressure, gradient, display. The shaders use GLSL ES 1.00 syntax, which a WebGL2 context accepts.

There are two grids. The simulation grid SIM holds the velocity, the curl, the divergence, and the pressure, with 160 cells on the short side of the canvas. The dye grid DYE has 1024 texels on the short side. The long side is the short side multiplied by the aspect ratio, so the cells are square.

The element count in the statistics (this.count) is the number of simulation cells, s.w × s.h. It is not the number of dye texels. A 1710 × 930 viewport gives 294 × 160 = 47,040 cells. initFBOs calculates the count again after each resize.

# Interaction

• Move. When the pointer moves more than 0.5 px between two frames, the module adds splats. A splat is a gaussian spot that is added to the velocity texture and then to the dye texture. The velocity of the splat is the movement as a fraction of the viewport, multiplied by FORCE = 6000. The y sign is inverted, because texture coordinates go up.
• Move fast. The module does not read pointer.vx, pointer.vy, or pointer.speed. It uses its own value moved, the distance in px since the last frame. A larger moved gives a larger force and brighter dye: the colour factor is 0.18 + 0.5 · min(1, moved / 40). A stroke longer than 12 px is divided into 2 to 4 splats along its line, so that a fast stroke is not a row of spots. Each of these splats adds the full force and the full colour.
• Colour. hue increases by 0.05 for each second and by 0.0004 for each px of travel. One full turn through the five PALETTE colours thus takes 20 s at rest, or 2,500 px of movement.
• Hold the button. While the button is down, the radius parameter of the stroke splats is multiplied by 2.2, which makes the spot approximately 1.5 times wider. A held button with no pointer movement adds nothing.
• Press edge. On the frame where pointer.down becomes true, the module adds 14 splats at the cursor position, with a radius scale of 1.6. Their velocities point outward at 14 equal angles with a magnitude of 2,500 texels/s, and their colours go one time around the palette. All 14 splats are at the same point, so the 14 velocity vectors add to zero. Only a rounding error of the half-float texture stays: in a test, approximately 26 texels/s at the centre, against 2,494 texels/s after one such splat. The dye does add, to a bright spot that is almost white. The result is a white flash that the flow then moves away, not an outward blast.
• Idle. After 2.5 s with no input, the shared pointer goes into autopilot. The module sees this as normal movement. The path is slow, so the splats are weak and dim.

init does not take the reduced flag from its argument, so the module does not store it and update does not read it. Reduced motion does not change this effect. The module reads no URL parameter. There is one mobile branch: a viewport narrower than 800 px on a device with a touch screen (navigator.maxTouchPoints > 0) gets SIM = 96 and DYE = 512. init makes this decision one time, and a resize does not make it again.

# Model

Stable fluids

The Navier-Stokes equations for an incompressible fluid have two parts: a momentum equation, and the condition that the divergence of the velocity is zero. Stam's method solves them in separate steps that are stable for all time steps: add forces, move the field along itself (advection), and remove the divergence (projection). This module has no viscosity step, because the bilinear filter of the advection already smooths the field.

The velocity unit is simulation texels for each second. In the formulas below, [L] [R] [T] [B] are the left, right, top, and bottom neighbour cells, and [C] is the cell itself. The grid distance is 1.

Neighbour coordinates

The vertex shader calculates the texture coordinates of the four neighbours, vL vR vT vB, from the uniform texelSize, and the GPU interpolates them as varyings. The offsets are the same for each fragment, so this work is done one time for each vertex and not one time for each fragment. All nine fragment shaders share the header H that declares the varyings.

Splat

p      = (vUv − point) · (aspect, 1)
target = target + exp(−(p · p) / r) · color

The multiplication of p.x by the aspect ratio makes the spot round on a canvas that is not square. The gaussian is 1/e at a distance of √r of the canvas height. With RADIUS = 0.0022 on a 16:9 canvas, this is approximately 6% of the height. For the velocity target, color is (fx, fy, 0). For the dye target, it comes from color(h), which reads the five-entry PALETTE as a ring and mixes the two entries nearest to h linearly, so the hue moves with no steps.

Curl and vorticity confinement

curl = 0.5 · ((vy[R] − vy[L]) − (vx[T] − vx[B]))          ≈ ∂vy/∂x − ∂vx/∂y

f    = 0.5 · (|curl[T]| − |curl[B]|,  |curl[R]| − |curl[L]|)
f    = f / (|f| + 0.0001)
f    = f · CURL · curl[C]
f.y  = −f.y
v    = clamp(v + f · dt, −1000, 1000)

The advection step loses small vortices, because each bilinear read is a small blur. Vorticity confinement puts that rotation back. Before the multiplication, f is the gradient of the curl magnitude, turned by 90° and made a unit vector. It points along the rim of a vortex. The multiplication by the signed curl selects the direction in which the vortex already turns, so the force makes each vortex turn faster.

Projection: divergence, pressure, gradient

div = 0.5 · ((vx[R] − vx[L]) + (vy[T] − vy[B]))
    wall: a neighbour outside the texture counts as −v[C]

p   ← 0.8 · p                                        (clear pass, one time)
p   ← (p[L] + p[R] + p[B] + p[T] − div) / 4          (pressure pass, 20 times)

v   = v − (p[R] − p[L],  p[T] − p[B])

The divergence tells how much fluid a cell makes or loses. The pressure is the solution of the Poisson equation ∇²p = div, and the subtraction of its gradient from the velocity removes the divergence. This is the step that makes the fluid roll and not compress. The wall rule mirrors the velocity at the screen edge, so the pressure pushes the flow back from the walls.

One Jacobi iteration moves information by one cell, so 20 iterations from zero cannot reach across 160 cells. The clear pass is the answer to this. It does not set the pressure to zero. It keeps 80% of the pressure of the last frame as the start value. The pressure changes little between frames, so the iterations start near the solution. The factor is less than 1 so that an old error fades and does not collect.

Note that the gradient pass subtracts p[R] − p[L] and not 0.5 · (p[R] − p[L]), which is two times the central-difference gradient. The same form is in the WebGL Fluid Simulation of Pavel Dobryakov (see the references). An incomplete Jacobi solve gives a pressure that is too small at large scales, and the larger step removes more of that divergence on each frame. The result is not an exact projection, but it is sufficient for a visual effect.

Advection and dissipation

src = vUv − dt · v(vUv) · texelSize
q   = q(src) / (1 + dissipation · dt)

The shader does not push a value forward to a new cell. It goes back along the velocity to the point where the value came from, and reads the field there with the bilinear filter of the GPU. The result is always a mix of old values, so the step is stable for all time steps. The velocity is advected by itself first, then the dye by the new velocity. The dye pass also uses the texel size of the simulation grid, because the velocity unit is simulation texels.

The division is an implicit Euler step of an exponential decay. The half-life is ln 2 / dissipation: 0.77 s for the dye and 2.8 s for the velocity.

Display

dx  = |dye[R]| − |dye[L]|,   dy = |dye[T]| − |dye[B]|
n   = normalize(dx, dy, 2 · |texelSize|)
c   = dye · clamp(n.z + 0.7, 0.7, 1.0)
c   = 1 − exp(−1.8 · c)
c   = c + 0.6 · max(c − 0.6, 0)²
out = c · (0.97 + 0.03 · sin(1.7 · y)) + (0.02, 0.02, 0.03)

The dye brightness is used as a height map: steep edges become darker by up to 30%. The exponential curve compresses the additive dye into the range 0 to 1 and keeps the colours saturated. The last line adds scanlines with a period of 3.7 device pixels and a dark blue floor.

# Frame pipeline

1. Stroke splatsreads pointer.x, pointer.y, velocity, dyewrites velocity, dye

   Only when moved > 0.5. 1 to 4 splats. Each splat is two draws: one into velocity.write and one into dye.write, each followed by a swap.

2. Press burstreads pointer.down, wasDown, velocity, dyewrites velocity, dye

   Only on the press edge. 14 splats at the cursor position.

3. curlreads velocitywrites curl

   The scalar curl of the velocity, one value for each cell.

4. vorticityreads velocity, curlwrites velocity

   Adds the confinement force f · dt and limits the velocity to ±1000 texels/s.

5. divergencereads velocitywrites divergence

   Central differences, with the mirror rule at the walls.

6. clearreads pressurewrites pressure

   Multiplies the pressure of the last frame by 0.8.

7. pressure, 20 timesreads divergence, pressurewrites pressure

   One Jacobi iteration for each draw. The pair of pressure textures swaps after each draw.

8. gradientreads pressure, velocitywrites velocity

   Subtracts the pressure gradient. The velocity now has much less divergence.

9. advect the velocityreads velocitywrites velocity

   The velocity texture is the velocity input and also the source input. Dissipation 0.25.

10. advect the dyereads velocity, dyewrites dye

    The same program with the dye as the source. Dissipation 0.9. This is the only solver pass on the dye grid.

11. displayreads dyewrites canvas

    blit(null) draws into the default framebuffer with the viewport of the full canvas.

Before the first step, update disables blending and limits the time step of the solver passes: sdt = min(dt, 1 / 40).

The header comment of the module gives the same order but does not list the clear pass. The projection runs before the advection, so the field that moves the dye has low divergence.

# Code

Vertex shader and shared header

Look at the four neighbour varyings. texelSize is a uniform, so the same vertex shader serves the simulation grid and the dye grid. H is put in front of each fragment shader.

const VS = `precision highp float; attribute vec2 aPosition; varying vec2 vUv, vL, vR, vT, vB; uniform vec2 texelSize;
void main(){ vUv = aPosition*0.5+0.5; vL = vUv-vec2(texelSize.x,0.); vR = vUv+vec2(texelSize.x,0.); vT = vUv+vec2(0.,texelSize.y); vB = vUv-vec2(0.,texelSize.y); gl_Position = vec4(aPosition,0.,1.); }`;
const H = `precision highp float; precision highp sampler2D; varying vec2 vUv, vL, vR, vT, vB;`;

Field shaders

These are the seven shaders of the solver. Compare each main with its formula in the model section. In divergence, the four if statements are the wall rule. In gradient, there is no factor 0.5.

splat: H + `uniform sampler2D uTarget; uniform float aspectRatio, radius; uniform vec3 color; uniform vec2 point;
  void main(){ vec2 p = vUv-point; p.x*=aspectRatio; vec3 s = exp(-dot(p,p)/radius)*color; vec3 b = texture2D(uTarget,vUv).xyz; gl_FragColor = vec4(b+s,1.); }`,
advect: H + `uniform sampler2D uVelocity, uSource; uniform vec2 texelSize; uniform float dt, dissipation;
  void main(){ vec2 c = vUv - dt*texture2D(uVelocity,vUv).xy*texelSize; vec4 r = texture2D(uSource,c); gl_FragColor = r/(1.0+dissipation*dt); }`,
divergence: H + `uniform sampler2D uVelocity;
  void main(){ float L=texture2D(uVelocity,vL).x, R=texture2D(uVelocity,vR).x, T=texture2D(uVelocity,vT).y, B=texture2D(uVelocity,vB).y; vec2 C=texture2D(uVelocity,vUv).xy;
    if(vL.x<0.)L=-C.x; if(vR.x>1.)R=-C.x; if(vT.y>1.)T=-C.y; if(vB.y<0.)B=-C.y; gl_FragColor = vec4(0.5*(R-L+T-B),0.,0.,1.); }`,
curl: H + `uniform sampler2D uVelocity;
  void main(){ float L=texture2D(uVelocity,vL).y, R=texture2D(uVelocity,vR).y, T=texture2D(uVelocity,vT).x, B=texture2D(uVelocity,vB).x; gl_FragColor = vec4(0.5*(R-L-T+B),0.,0.,1.); }`,
vorticity: H + `uniform sampler2D uVelocity, uCurl; uniform float curl, dt;
  void main(){ float L=texture2D(uCurl,vL).x, R=texture2D(uCurl,vR).x, T=texture2D(uCurl,vT).x, B=texture2D(uCurl,vB).x, C=texture2D(uCurl,vUv).x;
    vec2 f = 0.5*vec2(abs(T)-abs(B), abs(R)-abs(L)); f /= length(f)+1e-4; f *= curl*C; f.y*=-1.;
    vec2 v = texture2D(uVelocity,vUv).xy + f*dt; gl_FragColor = vec4(clamp(v,-1000.,1000.),0.,1.); }`,
pressure: H + `uniform sampler2D uPressure, uDivergence;
  void main(){ float L=texture2D(uPressure,vL).x, R=texture2D(uPressure,vR).x, T=texture2D(uPressure,vT).x, B=texture2D(uPressure,vB).x; float d=texture2D(uDivergence,vUv).x;
    gl_FragColor = vec4((L+R+B+T-d)*0.25,0.,0.,1.); }`,
gradient: H + `uniform sampler2D uPressure, uVelocity;
  void main(){ float L=texture2D(uPressure,vL).x, R=texture2D(uPressure,vR).x, T=texture2D(uPressure,vT).x, B=texture2D(uPressure,vB).x;
    vec2 v = texture2D(uVelocity,vUv).xy - vec2(R-L,T-B); gl_FragColor = vec4(v,0.,1.); }`,

Framebuffers and texture formats

A shader cannot read the texture that it writes. dfbo thus makes two equal targets, read and write, and swap() exchanges them after each draw. fbo makes one texture with its framebuffer. Its attach(id) helper binds the texture to a texture unit and returns the unit number for the sampler uniform.

initFBOs gives each field the smallest format that holds it: RGBA16F for the dye, RG16F for the velocity, R16F for the scalars. The dye and the velocity use LINEAR, because the advection depends on bilinear reads. The scalar fields use NEAREST, because the shaders read them only at cell centres. The divergence and the curl are single targets, because no pass reads and writes them at the same time. A render to a half-float texture needs the extension EXT_color_buffer_float.

fbo(w, h, internal, format, type, filter) {
  const gl = this.gl;
  const tex = gl.createTexture(); gl.bindTexture(gl.TEXTURE_2D, tex);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, filter); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, filter);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE);
  gl.texImage2D(gl.TEXTURE_2D, 0, internal, w, h, 0, format, type, null);
  const fb = gl.createFramebuffer(); gl.bindFramebuffer(gl.FRAMEBUFFER, fb);
  gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, tex, 0);
  gl.viewport(0, 0, w, h); gl.clear(gl.COLOR_BUFFER_BIT);
  return { tex, fb, w, h, tx: 1 / w, ty: 1 / h, attach(id) { gl.activeTexture(gl.TEXTURE0 + id); gl.bindTexture(gl.TEXTURE_2D, tex); return id; } };
},
dfbo(w, h, i, f, t, filter) {
  const a = this.fbo(w, h, i, f, t, filter), b = this.fbo(w, h, i, f, t, filter);
  return { read: a, write: b, w, h, tx: 1 / w, ty: 1 / h, swap() { const t = this.read; this.read = this.write; this.write = t; } };
},

initFBOs() {
  const gl = this.gl, s = this.res(this.cfg.SIM), d = this.res(this.cfg.DYE);
  const HF = gl.HALF_FLOAT;
  this.dye = this.dfbo(d.w, d.h, gl.RGBA16F, gl.RGBA, HF, gl.LINEAR);
  this.velocity = this.dfbo(s.w, s.h, gl.RG16F, gl.RG, HF, gl.LINEAR);
  this.divergence = this.fbo(s.w, s.h, gl.R16F, gl.RED, HF, gl.NEAREST);
  this.curl = this.fbo(s.w, s.h, gl.R16F, gl.RED, HF, gl.NEAREST);
  this.pressure = this.dfbo(s.w, s.h, gl.R16F, gl.RED, HF, gl.NEAREST);
  this.count = s.w * s.h;
},

Colour and splat

color mixes two neighbour palette entries, and k scales the brightness. splat uses one program for two targets and changes only the uTarget and color uniforms between the two draws. The aspect correction of the radius is on the third line of splat.

color(h, k = 1) {
  // blend between neighbouring palette entries → smooth neon hue drift
  const n = PALETTE.length, f = ((h % 1) + 1) % 1 * n, i = Math.floor(f), t = f - i;
  const a = PALETTE[i % n], b = PALETTE[(i + 1) % n];
  return [(a[0] + (b[0] - a[0]) * t) * k, (a[1] + (b[1] - a[1]) * t) * k, (a[2] + (b[2] - a[2]) * t) * k];
},

splat(x, y, dx, dy, color, rScale = 1) {
  const gl = this.gl, u = this.use('splat');
  const aspect = this.canvas.width / this.canvas.height;
  let r = this.cfg.RADIUS * rScale; if (aspect > 1) r *= aspect;
  gl.uniform1i(u.uTarget, this.velocity.read.attach(0));
  gl.uniform1f(u.aspectRatio, aspect); gl.uniform2f(u.point, x, y); gl.uniform3f(u.color, dx, dy, 0); gl.uniform1f(u.radius, r);
  this.blit(this.velocity.write); this.velocity.swap();
  gl.uniform1i(u.uTarget, this.dye.read.attach(0));
  gl.uniform3f(u.color, color[0], color[1], color[2]);
  this.blit(this.dye.write); this.dye.swap();
},

Pointer input

In the sub-step loop, k goes from 1 / steps to 1, and the splat position goes from near the last pointer position to the current position. The burst loop calls splat 14 times with the same nx, ny, which is the reason that the velocities cancel.

// ---- input ----
const dxp = p.x - this.lastX, dyp = p.y - this.lastY;
this.lastX = p.x; this.lastY = p.y;
const moved = Math.hypot(dxp, dyp);
this.hue += dt * 0.05 + moved * 0.0004;
if (moved > 0.5) {
  const nx = p.x / this.width, ny = 1 - p.y / this.height;
  const fx = dxp / this.width * c.FORCE, fy = -dyp / this.height * c.FORCE;
  const strength = Math.min(1, moved / 40);
  const col = this.color(this.hue, 0.18 + 0.5 * strength);
  // sub-steps along the stroke keep fast strokes continuous
  const steps = Math.min(4, Math.ceil(moved / 12));
  for (let i = 1; i <= steps; i++) {
    const k = i / steps;
    this.splat(nx - dxp / this.width * (1 - k), ny + dyp / this.height * (1 - k), fx, fy, col, p.down ? 2.2 : 1);
  }
}
if (p.down && !this.wasDown) {
  // burst: a ring of splats blasting outward
  const nx = p.x / this.width, ny = 1 - p.y / this.height;
  for (let i = 0; i < 14; i++) { const a = (i / 14) * Math.PI * 2; this.splat(nx, ny, Math.cos(a) * 2500, Math.sin(a) * 2500, this.color(this.hue + i / 14, 0.5), 1.6); }
}
this.wasDown = p.down;

The solver passes

Each pass has the same pattern: use the program, set texelSize, attach the input textures, blit into the target, swap if the target is a pair. The pressure loop sets only the uPressure sampler again on each iteration, because the read texture changes with each swap.

// ---- solve ----
const V = this.velocity;
let u = this.use('curl');
gl.uniform2f(u.texelSize, V.tx, V.ty); gl.uniform1i(u.uVelocity, V.read.attach(0)); this.blit(this.curl);

u = this.use('vorticity');
gl.uniform2f(u.texelSize, V.tx, V.ty); gl.uniform1i(u.uVelocity, V.read.attach(0)); gl.uniform1i(u.uCurl, this.curl.attach(1));
gl.uniform1f(u.curl, c.CURL); gl.uniform1f(u.dt, sdt); this.blit(V.write); V.swap();

u = this.use('divergence');
gl.uniform2f(u.texelSize, V.tx, V.ty); gl.uniform1i(u.uVelocity, V.read.attach(0)); this.blit(this.divergence);

u = this.use('clear');
gl.uniform1i(u.uTexture, this.pressure.read.attach(0)); gl.uniform1f(u.value, c.PRESSURE); this.blit(this.pressure.write); this.pressure.swap();

u = this.use('pressure');
gl.uniform2f(u.texelSize, V.tx, V.ty); gl.uniform1i(u.uDivergence, this.divergence.attach(0));
for (let i = 0; i < c.ITERS; i++) { gl.uniform1i(u.uPressure, this.pressure.read.attach(1)); this.blit(this.pressure.write); this.pressure.swap(); }

u = this.use('gradient');
gl.uniform2f(u.texelSize, V.tx, V.ty); gl.uniform1i(u.uPressure, this.pressure.read.attach(0)); gl.uniform1i(u.uVelocity, V.read.attach(1)); this.blit(V.write); V.swap();

u = this.use('advect');
gl.uniform2f(u.texelSize, V.tx, V.ty);
gl.uniform1i(u.uVelocity, V.read.attach(0)); gl.uniform1i(u.uSource, V.read.attach(0));
gl.uniform1f(u.dt, sdt); gl.uniform1f(u.dissipation, c.VEL_DISS); this.blit(V.write); V.swap();
gl.uniform1i(u.uVelocity, V.read.attach(0)); gl.uniform1i(u.uSource, this.dye.read.attach(1));
gl.uniform1f(u.dissipation, c.DYE_DISS); this.blit(this.dye.write); this.dye.swap();

Resize and dispose

resize sets the canvas size, deletes all textures and framebuffers, and calls initFBOs. The new textures are empty, so a resize removes all dye and all motion. dispose deletes the nine programs and the vertex buffer, calls loseContext() of the WEBGL_lose_context extension, and removes the canvas. A browser permits only a limited number of live WebGL contexts, and a context stays alive until the garbage collector finds it. The explicit loss releases the GPU memory immediately when the visitor changes to a different effect.

# Tuned constants

# Performance

All field work is on the GPU. On the CPU, update() only sets uniforms and sends draw calls: 28 for a frame with no input, and 2 more for each splat, of which there are at most 18 on one frame. This is the reason that the measured CPU time is near zero. The GPU time is not in the measured values.

The GPU cost has three parts:

• Solver passes. 26 passes on the simulation grid. The cell count depends only on SIM and the aspect ratio, not on the viewport size. At 47,040 cells, this is approximately 1.2 million fragments for each frame.
• Dye advection. One pass on the dye grid: 1883 × 1024, approximately 1.9 million fragments, for a 1710 × 930 viewport.
• Display. One pass on the canvas with five texture reads for each pixel. This is the only part that scales with the viewport and the pixel ratio: 2565 × 1395, approximately 3.6 million pixels, for a 1710 × 930 viewport.

Limits that the module applies:

• The device pixel ratio is limited to 1.5. On the display of the measurements below (ratio 2), the canvas thus has 1.5 device pixels for each CSS pixel.
• The grids have a fixed size on the short side, so a large monitor does not make the solver slower.
• The mobile branch uses SIM = 96 and DYE = 512: 36% of the cells and 25% of the dye texels.
• All textures are half-float. They use half the memory of 32-bit float textures, and WebGL2 can filter them linearly with no extension.
• The context has no alpha channel, no antialiasing, no depth buffer, and no stencil buffer.
• There is no fallback. If WebGL2 or EXT_color_buffer_float is not available, init throws an error, and the controller reports it.

The cell count is different in the first row because the aspect ratio is different: 213 × 160 against 294 × 160.

# References

• Jos Stam, Stable Fluids (SIGGRAPH 1999). The source of the method: semi-Lagrangian advection and projection.
• Jos Stam, Real-Time Fluid Dynamics for Games (GDC 2003). A shorter description of the same solver, with full C code.
• Mark J. Harris, Fast Fluid Dynamics Simulation on the GPU (GPU Gems, chapter 38). The method as fragment shader passes, with vorticity confinement.
• Pavel Dobryakov, WebGL Fluid Simulation. A WebGL solver with the same pass structure, the same pressure carry-over of 0.8, and the same gradient form. The repository has a link to the live demonstration.
• Wikipedia, Vorticity confinement. The derivation of the confinement force.