Fluid Dynamics, One Term at a Time
Contents
KO 한국어
u(x,t)DuDt=∂u∂t+(u⋅∇)u=−∇pρ+ν∇2u+f,∇⋅u=0\vel{\mathbf{u}(\mathbf{x},t)}\qquad \mat{\frac{D\mathbf{u}}{Dt}} = \frac{\partial\mathbf{u}}{\partial t} + \adv{(\mathbf{u}\cdot\nabla)\mathbf{u}} = \prs{-\frac{\nabla p}{\rho}} + \vis{\nu\nabla^{2}\mathbf{u}} + \frc{\mathbf{f}},\qquad \cont{\nabla\cdot\mathbf{u} = 0}

Chapter 6

The Full Lab

This is chapter 0’s lab again. When you first stirred the water in that box, you did not know what the five chips at the top left of the experiment window meant. Now you know what each one does. Chapter 0 promised that by the last chapter the terms you build up would add up to a fluid simulator, and this is that chapter. First it covers the one term no chapter has taken on yet, the force from outside. Then you switch the lab’s terms off one at a time and see again what each one does, and you read one flow in which all of them work together: the flow past a wing. Last, you open the code of the engine that runs the lab, stage by stage. The equation bar highlights every term.

The force from outside

The f\frc{\mathbf{f}} at the right end of the equation is a force from outside the water, written per unit mass. Chapter 0 said that the hand that stirred the water is this term. The engine adds velocity to the water around the hand’s path, in the direction the hand moved. It adds it in a bell shape that falls off quickly away from the hand, so the water right under the hand gets the strongest push. Some outside forces act on all the water at once. Gravity does, and so does buoyancy, which makes warmed water rise through the colder water around it because it is lighter. The lab in this chapter has only the force of the stirring hand.

A force that pushes all the water evenly does not always start a flow. Water in a glass feels gravity and stays still. The pressure rises with depth, and that difference in pressure pushes the water up just as hard as gravity pulls it down. In chapter 3, pressure was the field that adjusts everywhere so that the water neither crowds together nor thins out. The same pressure holds up water that gravity would crowd downwards. The engine does the same. Give the water in chapter 0’s closed box a push downwards everywhere, and the pressure settles into layers that rise with depth while the water stays still.

Side note

The engine ran chapter 0’s box (a 480×270 grid with one cylinder, steps of 1/120 s) for 10 s with a body force of 1 downwards. The pressure p/ρp/\rho rose by 1 for every 1 of depth (to within 0.001 %), and the water stayed nearly still: the fastest water never passed 0.008 in the 10 s. That is close to the velocity the force adds in one step (the force 1 times the step 1/120 s). With steps of 1/60 s it was 0.017, and with 1/240 s 0.004, so it changed in proportion to the step. With the pressure projection off, the same box sent its fastest water rushing to the floor at a speed of 1.4 within 1.5 s, and it got no faster after that. Buoyancy is a force set by the water’s temperature or saltiness. This engine’s ink is only a marker that does not change the flow (chapter 0), so it has no buoyancy.

Switching off one term at a time

The experiment below takes out chapter 0’s lab again, with two differences. The water comes in on the left and leaves on the right, so every obstacle you put in the water leaves a wake. And every term has a switch. Four switches, for force, advection, viscosity and pressure, turn the terms off and on, and the chips at the top left of the experiment window always show the terms that are on right now. The script stirs in front of the cylinder a few times with the same stroke and switches the terms off one by one. Then it is your turn. Drag with the mouse (swipe with a finger on a phone) to stir, or set the hand to Draw and the drag draws an obstacle along its path. You can also pick a cylinder, a plate or a bean shape under Obstacle, and drag an obstacle with the hand on Stir. Picking another obstacle keeps the drawn lines.

Side note

All values are dimensionless. The water comes in at a speed of 0.3, and the cylinder’s diameter is 0.16. At the starting viscosity ν\nu of 0.0005 the cylinder’s Reynolds number is about 100, and the Viscosity ν slider takes ν\nu from 0.0001 to 0.002 (a Reynolds number from 480 down to 24). One pressure switch turns −∇p/ρ\prs{-\nabla p/\rho} and ∇⋅u=0\cont{\nabla\cdot\mathbf{u} = 0} off and on together, because keeping the water from crowding or thinning is exactly the pressure’s job (chapter 3). The ink rides the flow whichever term is off, because the engine moves the ink every step apart from the advection term. A browser without WebGPU shows a recording instead, and in the recording you cannot change the switches or stir the water.

At first every term is on. The script stirs upwards in front of the cylinder, the outside force pushes the water into a swirl, and the flow carries the swirl round the cylinder and away.

Switch the force off, and the same stroke does nothing. The stroke pushes no water and puts in no ink, because in the engine this stage also puts in the hand’s ink. The hand that stirred the water in chapter 0 was this term.

The script switches advection off and turns the vorticity layer on. Its orange and sky blue are the directions the water turns, and have nothing to do with the chips of the same colours. The stirred swirl stands where the hand left it instead of riding the flow. In 3 s its centre moved only 0.013. With advection on, the same stroke’s swirl drifts 0.51 downstream in 2 s. The ink and the particles do not stop, though: they flow past the standing swirl. It is chapter 5’s vortex pair again, which stopped where it was when advection went off. The flow still carries other things along; only its carrying of its own velocity is off.

Switch advection back on and raise the viscosity four times, from 0.0005 to 0.002, and the cylinder’s Reynolds number falls from about 100 to 24. A swirl from the same stroke soon dies away. Two seconds after the stroke, the strongest vorticity was 10 % of what it had been. At the starting viscosity of 0.0005, 26 % is left. Viscosity shares velocity evenly between neighbours, so it hands the swirl’s turning out to the water around it and spreads it (chapter 4).

Switch the viscosity off, and the swirl lasts longer: 37 % is left after 2 s. The layers that leave the cylinder roll up more finely too and shed swirls. Yet even with the viscosity off, the swirl slowly weakens. In a two-dimensional flow without viscosity, every parcel of water keeps its vorticity as it moves, so even the strongest point should not weaken. What you see is chapter 5’s numerical error in the engine, the error that comes from splitting each step into stages (side note).

Side note

Stir water with the viscosity off using the script’s stroke, and 2 s later the strongest vorticity keeps 37 % with steps of 1/60 s (this lab), 50 % with steps of 1/120 s (chapter 0’s lab) and 62 % with steps of 1/240 s. Halving the grid in each direction gives 40 %, and doubling it gives 33 %. The length of the step changes the share left by far more than the grid does, and a finer grid does not raise it. So, as in chapter 5, the loss reads as the error that comes from splitting each step into stages.

Switch the pressure off, and the water empties out behind the cylinder, and particles crowd into bright lines along the edges of the empty patch. As in chapter 2 with the projection off, the water crowds together in one place and thins out in another. With pressure, the pressure rises in front of the cylinder and turns the water aside in advance, and drops behind it and draws the water around into the gap (chapter 3). Without pressure, the cylinder only blocks the water, which does not go round it, and no water fills the space behind it. Switch the pressure back on, and the divergence soon drops back near zero; a few seconds later a wake forms behind the cylinder again, as it did at first.

Side note

After 3.5 s with the pressure off, the divergence times the cell size hh over the inflow speed UU ran from −0.68 to 0.09; 1.5 s after the pressure came back on, it was within 0.000005 of zero. Each of the 16 combinations of the switches ran for 30 s of stirring, and none blew up. The fastest was the force alone, at 26 times the inflow speed. The engine treats a flow as blown up beyond 50 times.

Now try it yourself. Switch the terms off and on and stir, then set the hand to Draw and draw a line across the flow. The cells along the line turn solid, and the water goes round it. With the cylinder removed and a line drawn straight across the flow, 3 s later the water right behind the line runs backwards (at about 0.4 times the inflow speed), and the water going round its end speeds up to 1.8 times the inflow speed.

Side note

A drawn line is 0.03 thick, and one stroke takes up to 48 points. Up to six lines stay; draw a seventh and the first one goes. Restart clears them all. With only a keyboard, pick a shape under Obstacle instead. The recording cannot draw obstacles.

How a wing lifts

An aeroplane’s wing only cuts forwards through the air, yet it holds up a heavy aeroplane. That force is lift. In the experiment below, water flows from the left towards a wing. Holding the wing still and letting the water flow past it shows the same thing as a wing moving through still air. The wing’s cross-section is NACA 0006, one of the shapes that NACA (the US National Advisory Committee for Aeronautics) defined by formula in the 1930s. It is symmetric top to bottom, its thickest point lies 30 % of the chord (the length from the leading edge to the trailing edge) back from the leading edge, and it is 6 % of the chord thick. The angle at which the wing meets the oncoming flow is its angle of attack, and Angle of attack picks 0°, 4°, 8° or 16°. It starts at 8°. Every 3 s, the experiment lays a white line of ink across the flow in front of the wing. It marks water that set off at the same moment.

Side note

All values are dimensionless. The water comes in at a speed UU of 0.5, the wing’s chord cc is 0.4 and the viscosity ν\nu is 0.0002. The Reynolds number on the chord, Uc/νUc/\nu, is 1000. Real aeroplane wings fly at Reynolds numbers of about 10510^5 to 10710^7, a hundred to ten thousand times more. At this low Reynolds number, the layer of water slowed down on the wing’s surface is thick. The flow is two-dimensional: one cross-section of an endlessly long wing, with no swirls rolling off its tips. On the desktop’s 480×270 grid the chord is 108 cells and the thickness 6.5 cells; on a phone’s 240×135 grid, half that. The engine draws the wing by marking cells solid, so its surface steps like a staircase. The slip walls above and below are 0.5 from the wing’s pivot (a quarter of the chord back), so the wing is not quite a wing in endless open air. The white line is sixteen dots of ink laid at once, from a height of 0.3 to 0.7, about 0.38 in front of the wing, and it pushes no water. You cannot drag the wing or stir the water. A browser without WebGPU shows a recording instead, and in the recording you cannot change the angle.

As the water starts to flow, a swirl rolls up at the trailing edge, comes off and drifts away. In the vorticity layer it is orange, a swirl turning anticlockwise. As it goes, it leaves behind a flow turning the other way, clockwise, round the wing. Half a second after the water started, the circulation measured round a box about the wing was 0.038 clockwise, and round a box about the water behind the wing it was 0.038 anticlockwise: the same size, the opposite way. This clockwise flow round the wing adds to the water above the wing in the direction of the flow, and to the water below against it, so it speeds the flow above up and slows the flow below down. This first swirl is the starting vortex.

Side note

Circulation is how chapter 5 measured a vortex’s strength: the mean speed of the water along a loop times the loop’s length. In a fluid without viscosity, the circulation round a loop that moves with the water does not change with time. That is Kelvin’s circulation theorem. In the uniform flow at the start, every loop had zero circulation, so after the start, the circulation round the wing and the starting vortex’s circulation must still add up to zero. This experiment has viscosity, but the two circulations measured half a second after the start were the same size to within 0.1 %.

Now follow the white line. At the leading edge it splits above and below the wing. The water passing above reaches the trailing edge first. At that moment the water passing below is still only 83 % of the way along the chord. The water above reached the trailing edge 1.42 s after the experiment laid the line, the water below 1.55 s after: 0.13 s apart. The water below never catches up, even behind the wing: 0.3 past the trailing edge it is still 0.13 s behind. Water that sets off together does not arrive together.

The gap comes from the faster flow above the wing. At mid-chord, 0.1 from the surface, the speed was 1.16 times UU above and 1.03 times UU below. In the pressure layer, the pressure is low above the wing (sky blue) and high below it (orange). At the same places, p/ρp/\rho was −0.035 above and 0.007 below. The pressure is highest just under the leading edge, where the flow runs into the wing and stops (0.12), and lowest just above it (−0.14). Add up the push of the pressure over the whole surface of the wing, and it points up. That is lift. To compare lifts, divide by the oncoming flow’s 12ρU2\tfrac12\rho U^2 and the chord to get the lift coefficient.

CL=L12ρU2cC_L = \frac{L}{\tfrac{1}{2}\rho U^2 c}

LL is the lift per unit width of the wing. In this experiment at 8°, CLC_L was 0.49.

Side note

Thin aerofoil theory estimates a thin wing’s lift coefficient at a high Reynolds number as 2π2\pi times the angle of attack in radians: 0.88 at 8°. This experiment’s 0.49 is smaller. The thick slow layer on the surface at a Reynolds number of 1000 looks like one reason, but the nearby walls and the staircase surface play in too, and this experiment does not tell their shares apart. On a phone’s coarser grid it was 0.42, and the water above still reached the trailing edge 0.1 s first.

Faster flow together with lower pressure is chapter 3’s Bernoulli equation at work. The equation holds under four conditions: steady, inviscid, incompressible, along one streamline. A few seconds after the water started, the flow has settled (steady); outside the slow layers on the surface, viscosity hardly works (inviscid); and the water’s density is constant (incompressible). The fourth condition holds too. Every streamline above and below the wing comes from the uniform flow in front of it, where all the streamlines have the same pressure and speed. So the constant CC is the same on all of them, and you can compare above with below. Measured above, below and in front of the wing, p/ρ+12∣u∣2p/\rho + \tfrac12 |\mathbf{u}|^2 differed from its value in front of the wing by 2.1 % of 12U2\tfrac12 U^2 above and 2.4 % below (4.5 % between above and below). Each is small beside the 34 % of 12U2\tfrac12 U^2 by which the pressure itself differed between above and below.

Behind the wing the flow turns down. At 0.2 behind the trailing edge, from its height to 0.3 above it, the water flowed downwards at 0.035 times UU on average. The wing turned the water down, and by Newton’s third law (action and reaction) the water pushes the wing up. The pressure and the downturned flow are not two rival explanations of lift but two faces of one flow. The difference in pressure above and below the wing pushes the wing up and, at the same time, turns the water around it down. The water pushes on the wing through that pressure, and the wing pushes back on the water.

Switch the angle to 0°, and the wing sheds a swirl as the flow adjusts, and a few seconds later the flow settles again. A level symmetric wing is the same above and below, as if mirrored. The two parts of the line reach the trailing edge together (within one frame, 1/60 s), and the lift coefficient is below 0.001. There is no lift. Tilting the wing is what makes lift.

Tilt it to 16°, and the flow cannot follow the upper surface to the end: it comes away and sheds large swirls behind the wing. Take the water’s speed along the flow just above the rear half of the upper surface (within 0.06 of it) and look at the slowest place: at 4° it still ran forward, though only at 0.0035 times UU, so it was almost standing still; at 8° it ran backwards at 0.035 times UU, and at 16° at 0.42 times UU (the reading at 16° was taken 3.4 s after the tilt, and it changes with time while the flow sheds swirls). On a real wing, too large an angle makes the flow come away from the upper surface and the lift fall sharply. That is a stall.

One misconception

A popular explanation of lift goes like this: the upper side of a wing bulges more, so it is longer, and the water that splits at the leading edge must meet again at the trailing edge, so the water above must go faster. This explanation is wrong. Above all, in this experiment the water that set off together did not arrive together. The water above arrived first, and the water below never caught up. Nothing in the equations requires the two to meet again at the trailing edge. The difference in length the explanation relies on does not hold up either. This experiment’s NACA 0006 is a symmetric wing whose upper and lower surfaces are exactly the same length, so by that explanation it should make no lift at 8°. Its lift coefficient was 0.49. A paper aeroplane’s flat wings fly too, and aeroplanes can fly upside down. That the flow above the wing is faster and its pressure lower is right. The real speed is much higher than the explanation gives, so a speed worked out from arriving together gives far less lift than a real wing makes. What sets the speed is the whole flow round the wing: the way the wing turns the flow down.

Side note

NASA Glenn Research Center’s Beginner’s Guide to Aeronautics sets out the same objections on its page “Equal Transit Theory Interactive”: symmetric and flat wings make lift too, aeroplanes fly upside down, air molecules close together at the leading edge are not together at the trailing edge, and the lift this explanation predicts is much less than the lift observed.

One step of the engine

As chapter 0 said, the terms you built up are the lab’s engine. The engine moves on in short slices of time, one step at a time (in this chapter’s experiments a step is 1/60 s), and every step it runs five stages in a fixed order. Stages ① to ④ are in charge of the terms of the equation (④ of two terms), and ⑤ sets the velocity at the edges and on the obstacles. The picture below is that one step. Pick a stage, and the row along the top shows the fields the stage only reads (dashed outline) and the ones it writes (filled). A field it reads and then writes back in place is drawn filled too. Each stage’s real code follows the picture under the same number. The code is WGSL, WebGPU’s shader language, which runs on the graphics hardware, and what you see here are the very files the engine uses.

Side note

The engine keeps the velocity on the cells’ sides (faces) and the pressure and ink at the cells’ centres. The horizontal velocity sits on a cell’s left and right faces, the vertical velocity on its bottom and top faces. Staggered like this, the face velocities count directly what flows into and out of each cell (chapter 2’s box). The picture stretches one step to a 5 s loop. The boundary stage runs after each of the other stages, but the picture draws it once. The engine skips the stage of a switched-off term (stages ②, ③ and ④ take the boundary stage after them with them), and moving the ink always runs.

① Force

This stage adds the outside force f\frc{\mathbf{f}}. It adds the push on all the water (the body force), times one step’s time, to the velocity on every face, and it adds the stirring hand’s push and its ink in a bell shape centred on the hand’s path. It leaves alone the faces that touch an obstacle or the edges of the experiment window: stage ⑤, the boundaries, sets those.

f\frc{\mathbf{f}}forces.wgsl

src/lab/engine/shaders/forces.wgsl

// Stage 1 · external force f: body force plus Gaussian "splats" (stirring) of velocity and dye.
struct Splat {
  pos: vec2f,
  vel: vec2f,
  dye: vec4f,
  radius: f32,
  pad0: f32,
  pad1: f32,
  pad2: f32,
};

@group(0) @binding(1) var<storage, read_write> uF: array<f32>;
@group(0) @binding(2) var<storage, read_write> vF: array<f32>;
@group(0) @binding(3) var<storage, read_write> dyeF: array<vec4f>;
@group(0) @binding(4) var<storage, read> splats: array<Splat>;
@group(0) @binding(5) var<storage, read> solid: array<u32>;

fn splatWeight(p: vec2f, s: Splat) -> f32 {
  var d = p - s.pos;
  let L = domainSize();
  if (periodicX()) { d.x = d.x - L.x * round(d.x / L.x); }
  if (periodicY()) { d.y = d.y - L.y * round(d.y / L.y); }
  return exp(-dot(d, d) / (s.radius * s.radius));
}

@compute @workgroup_size(8, 8)
fn forces(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  if (i <= P.nx && j < P.ny && !uFixed(i, j)) {
    let p = vec2f(f32(i), f32(j) + 0.5) * P.h;
    var du = P.gx * P.dt;
    for (var s = 0u; s < P.splatCount; s++) { du += splats[s].vel.x * splatWeight(p, splats[s]); }
    uF[iU(i, j)] += du;
  }
  if (i < P.nx && j <= P.ny && !vFixed(i, j)) {
    let p = vec2f(f32(i) + 0.5, f32(j)) * P.h;
    var dv = P.gy * P.dt;
    for (var s = 0u; s < P.splatCount; s++) { dv += splats[s].vel.y * splatWeight(p, splats[s]); }
    vF[iV(i, j)] += dv;
  }
  if (i < P.nx && j < P.ny && solid[iC(i, j)] == 0u) {
    let p = (vec2f(f32(i), f32(j)) + 0.5) * P.h;
    var dd = vec4f(0.0);
    for (var s = 0u; s < P.splatCount; s++) { dd += splats[s].dye * splatWeight(p, splats[s]); }
    dyeF[iC(i, j)] = clamp(dyeF[iC(i, j)] + dd, vec4f(0.0), vec4f(1.0));
  }
}

② Advection

This stage solves (u⋅∇)u\adv{(\mathbf{u}\cdot\nabla)\mathbf{u}}, the flow carrying its own velocity along. For each face it traces the flow backwards to where the water that arrives there in this step came from, and brings the velocity from there (advect_vel.wgsl). The trace takes the velocity at the point half a step back and traces again with it. Running the same shader once more in reverse, to carry the result forwards along the flow again, does not quite give back the starting values; the second shader, correct_vel.wgsl, corrects by half that mismatch. Where the corrected value would leave the range of the values around it that it was brought from, the uncorrected value stays. The ink moves the same way, with a pair of its own (advect_dye.wgsl, correct_dye.wgsl), except that a corrected value outside that range is clipped to the range instead of being replaced by the uncorrected one. The ink moves at the very end of the step, with the flow after the pressure projection, and it always moves, even with the advection term off. That is why the ink kept flowing in the lab with advection off.

(u⋅∇)u\adv{(\mathbf{u}\cdot\nabla)\mathbf{u}}advect_vel.wgsl · correct_vel.wgsl · advect_dye.wgsl · correct_dye.wgsl

src/lab/engine/shaders/advect_vel.wgsl

// Stage 2 · inertia (u.grad)u, part 1: semi-Lagrangian transport of the velocity field.
// DT_SIGN = +1: source u^n traced backward -> hat.  DT_SIGN = -1: source hat traced forward -> back.
// Both passes trace through u^n (uIn / vIn), as MacCormack requires.
override DT_SIGN: f32 = 1.0;

@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> srcU: array<f32>;
@group(0) @binding(4) var<storage, read> srcV: array<f32>;
@group(0) @binding(5) var<storage, read_write> dstU: array<f32>;
@group(0) @binding(6) var<storage, read_write> dstV: array<f32>;
@group(0) @binding(7) var<storage, read> solid: array<u32>;

@compute @workgroup_size(8, 8)
fn advect_vel(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  let dts = DT_SIGN * P.dt;
  if (i <= P.nx && j < P.ny) {
    let k = iU(i, j);
    if (uFixed(i, j)) {
      dstU[k] = srcU[k];
    } else {
      dstU[k] = sampleU_src(trace(vec2f(f32(i), f32(j) + 0.5) * P.h, dts));
    }
  }
  if (i < P.nx && j <= P.ny) {
    let k = iV(i, j);
    if (vFixed(i, j)) {
      dstV[k] = srcV[k];
    } else {
      dstV[k] = sampleV_src(trace(vec2f(f32(i) + 0.5, f32(j)) * P.h, dts));
    }
  }
}

src/lab/engine/shaders/correct_vel.wgsl

// Stage 2 · inertia, part 2: MacCormack correction  u_new = hat + (u^n - back) / 2,
// falling back to hat where the result would leave the range of the values it was interpolated from.
@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read_write> hatU: array<f32>;
@group(0) @binding(4) var<storage, read_write> hatV: array<f32>;
@group(0) @binding(5) var<storage, read> backU: array<f32>;
@group(0) @binding(6) var<storage, read> backV: array<f32>;
@group(0) @binding(7) var<storage, read> solid: array<u32>;

@compute @workgroup_size(8, 8)
fn correct_vel(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  if (i <= P.nx && j < P.ny && !uFixed(i, j)) {
    let k = iU(i, j);
    let r = rangeU(trace(vec2f(f32(i), f32(j) + 0.5) * P.h, P.dt));
    let c = hatU[k] + 0.5 * (uIn[k] - backU[k]);
    hatU[k] = select(c, hatU[k], c < r.x || c > r.y);
  }
  if (i < P.nx && j <= P.ny && !vFixed(i, j)) {
    let k = iV(i, j);
    let r = rangeV(trace(vec2f(f32(i) + 0.5, f32(j)) * P.h, P.dt));
    let c = hatV[k] + 0.5 * (vIn[k] - backV[k]);
    hatV[k] = select(c, hatV[k], c < r.x || c > r.y);
  }
}

src/lab/engine/shaders/advect_dye.wgsl

// Dye transport (visualisation): the same semi-Lagrangian trace as the velocity, cell-centred RGBA.
override DT_SIGN: f32 = 1.0;

@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> dyeSrc: array<vec4f>;
@group(0) @binding(4) var<storage, read_write> dyeDst: array<vec4f>;

@compute @workgroup_size(8, 8)
fn advect_dye(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny) { return; }
  let p = (vec2f(f32(id.x), f32(id.y)) + 0.5) * P.h;
  dyeDst[iC(id.x, id.y)] = sampleDye(trace(p, DT_SIGN * P.dt));
}

src/lab/engine/shaders/correct_dye.wgsl

// Dye MacCormack correction, clamped per channel to the interpolation range (no new extrema).
@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> dyeSrc: array<vec4f>;
@group(0) @binding(4) var<storage, read_write> dyeHat: array<vec4f>;
@group(0) @binding(5) var<storage, read> dyeBack: array<vec4f>;

@compute @workgroup_size(8, 8)
fn correct_dye(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny) { return; }
  let k = iC(id.x, id.y);
  let p = (vec2f(f32(id.x), f32(id.y)) + 0.5) * P.h;
  let r = rangeDye(trace(p, P.dt));
  dyeHat[k] = clamp(dyeHat[k] + 0.5 * (dyeSrc[k] - dyeBack[k]), r.lo, r.hi);
}

③ Viscosity

This stage solves ν∇2u\vis{\nu\nabla^2\mathbf{u}}, the sharing of velocity between neighbours. The engine solves an equation that takes the differences from the neighbours with the new velocity, the one at the end of the step. The new velocity appears on both sides, so the equation does not solve in one go: the engine colours the cells like a chessboard and corrects one colour, then the other, over and over. The larger ν Δt/h2\nu\,\Delta t/h^2, the more rounds it needs. This chapter’s lab runs 4 to 25 of them depending on the viscosity (7 at the starting value).

ν∇2u\vis{\nu\nabla^{2}\mathbf{u}}diffuse.wgsl

src/lab/engine/shaders/diffuse.wgsl

// Stage 3 · viscosity nu*lap(u): one red-black Gauss-Seidel sweep of the implicit step
//   (1 + 4a) u - a * sum(u_nb) = u0,   a = nu * dt / h^2
// Neighbours beyond walls come from the ghost rules in fetchU_d / fetchV_d.
override COLOR: u32 = 0u;

@group(0) @binding(1) var<storage, read_write> uD: array<f32>;
@group(0) @binding(2) var<storage, read_write> vD: array<f32>;
@group(0) @binding(3) var<storage, read> u0: array<f32>;
@group(0) @binding(4) var<storage, read> v0: array<f32>;
@group(0) @binding(5) var<storage, read> solid: array<u32>;

@compute @workgroup_size(8, 8)
fn diffuse(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  if (((i + j) & 1u) != COLOR) { return; }
  let a = P.nu * P.dt / (P.h * P.h);
  let ii = i32(i);
  let jj = i32(j);
  if (i <= P.nx && j < P.ny && !uFixed(i, j)) {
    let s = fetchU_d(ii - 1, jj) + fetchU_d(ii + 1, jj) + fetchU_d(ii, jj - 1) + fetchU_d(ii, jj + 1);
    uD[iU(i, j)] = (u0[iU(i, j)] + a * s) / (1.0 + 4.0 * a);
  }
  if (i < P.nx && j <= P.ny && !vFixed(i, j)) {
    let s = fetchV_d(ii - 1, jj) + fetchV_d(ii + 1, jj) + fetchV_d(ii, jj - 1) + fetchV_d(ii, jj + 1);
    vD[iV(i, j)] = (v0[iV(i, j)] + a * s) / (1.0 + 4.0 * a);
  }
}

④ Projection

This stage solves −∇p/ρ\prs{-\nabla p/\rho} and ∇⋅u=0\cont{\nabla\cdot\mathbf{u} = 0} together. First it counts, for each cell, what flows in and out through its faces: the divergence (divergence.wgsl). Then it solves for the pressure that removes that divergence. When chapter 3 said that pressure adjusts everywhere at once, it meant that every cell has one equation tying its pressure to its neighbours’, and the engine has to solve the equations of all the cells together. By default the engine solves them on grids it coarsens by half again and again. It cuts the fine ripples of error on the fine grid, carries what is left to a coarser grid to cut the broad patterns of error cheaply, and brings the correction back to the fine grid (the multigrid method, the five files that start with mg_). When its setting is the one-grid method, or when halving the grid does not reach a coarsest grid small enough (4096 cells or fewer), pressure.wgsl solves on the one grid, round after round. Last, gradient.wgsl subtracts the pressure’s gradient from the velocity, which leaves a flow that neither crowds nor thins.

u←u−Δt ∇pρ\mathbf{u} \leftarrow \mathbf{u} - \Delta t\,\frac{\nabla p}{\rho}
−∇pρ,∇⋅u=0\prs{-\frac{\nabla p}{\rho}},\quad \cont{\nabla\cdot\mathbf{u} = 0}divergence.wgsl · mg_smooth.wgsl · mg_residual.wgsl · mg_restrict.wgsl · mg_coarsest.wgsl · mg_prolong.wgsl · pressure.wgsl · gradient.wgsl

src/lab/engine/shaders/divergence.wgsl

// Stage 4a · mass conservation check: net outflow of each fluid cell, (u_E - u_W + v_N - v_S) / h.
@group(0) @binding(1) var<storage, read> uIn: array<f32>;
@group(0) @binding(2) var<storage, read> vIn: array<f32>;
@group(0) @binding(3) var<storage, read> solid: array<u32>;
@group(0) @binding(4) var<storage, read_write> divOut: array<f32>;

@compute @workgroup_size(8, 8)
fn divergence(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny) { return; }
  let i = id.x;
  let j = id.y;
  let k = iC(i, j);
  if (solid[k] != 0u) {
    divOut[k] = 0.0;
    return;
  }
  divOut[k] = (uIn[iU(i + 1u, j)] - uIn[iU(i, j)] + vIn[iV(i, j + 1u)] - vIn[iV(i, j)]) / P.h;
}

src/lab/engine/shaders/mg_smooth.wgsl

// Multigrid smoother: one red-black Gauss-Seidel half-sweep of  sum(p_nb - p) / h^2 = f  on one level.
override COLOR: u32 = 0u;

@group(0) @binding(1) var<storage, read_write> p: array<f32>;
@group(0) @binding(2) var<storage, read> f: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;

@compute @workgroup_size(8, 8)
fn mg_smooth(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.nx || id.y >= L.ny || ((id.x + id.y) & 1u) != COLOR) { return; }
  let k = id.y * L.nx + id.x;
  if (mask[k] != 0u) { return; }
  let s = levelStencil(i32(id.x), i32(id.y));
  if (s.y == 0.0) { return; } // sealed off on all sides: nothing to balance
  p[k] = (s.x - L.h * L.h * L.fScale * f[k]) / s.y;
}

src/lab/engine/shaders/mg_residual.wgsl

// Multigrid residual  r = f - sum(p_nb - p) / h^2  on one level (0 in solid and sealed cells).
@group(0) @binding(1) var<storage, read> p: array<f32>;
@group(0) @binding(2) var<storage, read> f: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;
@group(0) @binding(4) var<storage, read_write> r: array<f32>;

@compute @workgroup_size(8, 8)
fn mg_residual(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.nx || id.y >= L.ny) { return; }
  let k = id.y * L.nx + id.x;
  let s = levelStencil(i32(id.x), i32(id.y));
  if (mask[k] != 0u || s.y == 0.0) {
    r[k] = 0.0;
    return;
  }
  r[k] = L.fScale * f[k] - (s.x - s.y * p[k]) / (L.h * L.h);
}

src/lab/engine/shaders/mg_restrict.wgsl

// Multigrid restriction: a coarse cell's right-hand side is the mean residual of its fine children inside
// the grid (one or two of them on the last row or column of an odd side); its correction starts at 0.
@group(0) @binding(1) var<storage, read> r: array<f32>;
@group(0) @binding(2) var<storage, read_write> fc: array<f32>;
@group(0) @binding(3) var<storage, read_write> pc: array<f32>;

@compute @workgroup_size(8, 8)
fn mg_restrict(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.cnx || id.y >= L.cny) { return; }
  var sum = 0.0;
  var count = 0.0;
  for (var b = 0u; b < 2u; b++) {
    for (var a = 0u; a < 2u; a++) {
      let i = 2u * id.x + a;
      let j = 2u * id.y + b;
      if (i < L.nx && j < L.ny) {
        sum += r[j * L.nx + i];
        count += 1.0;
      }
    }
  }
  let k = id.y * L.cnx + id.x;
  fc[k] = sum / count;
  pc[k] = 0.0;
}

src/lab/engine/shaders/mg_coarsest.wgsl

// Multigrid coarsest level: L.iterations red-black SOR sweeps by one workgroup of 256 invocations, with a
// storage barrier between colours so every cell sees its neighbours' new values.
@group(0) @binding(1) var<storage, read_write> p: array<f32>;
@group(0) @binding(2) var<storage, read> f: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;

@compute @workgroup_size(256)
fn mg_coarsest(@builtin(local_invocation_index) t: u32) {
  let n = L.nx * L.ny;
  for (var it = 0u; it < L.iterations; it++) {
    for (var color = 0u; color < 2u; color++) {
      for (var k = t; k < n; k += 256u) {
        let i = k % L.nx;
        let j = k / L.nx;
        if (((i + j) & 1u) == color && mask[k] == 0u) {
          let s = levelStencil(i32(i), i32(j));
          if (s.y > 0.0) {
            p[k] = mix(p[k], (s.x - L.h * L.h * L.fScale * f[k]) / s.y, L.omega);
          }
        }
      }
      storageBarrier();
    }
  }
}

src/lab/engine/shaders/mg_prolong.wgsl

// Multigrid prolongation: each fluid fine cell adds the bilinear interpolation of the coarse correction,
// its weights renormalised over the fluid coarse cells, so nothing leaks in from a solid one. The parent
// cell always carries a weight of at least 9/16, so the weights never sum to 0.
@group(0) @binding(1) var<storage, read_write> p: array<f32>;
@group(0) @binding(2) var<storage, read> pc: array<f32>;
@group(0) @binding(3) var<storage, read> mask: array<u32>;
@group(0) @binding(4) var<storage, read> maskC: array<u32>;

// Coarse cell index, or -1 beyond a non-periodic edge or inside a solid.
fn coarseCell(I: i32, J: i32) -> i32 {
  let cnx = i32(L.cnx);
  let cny = i32(L.cny);
  var ii = I;
  var jj = J;
  if (I < 0 || I >= cnx) {
    if (L.edgeL != LEVEL_PERIODIC) { return -1; }
    ii = ((I % cnx) + cnx) % cnx;
  }
  if (J < 0 || J >= cny) {
    if (L.edgeB != LEVEL_PERIODIC) { return -1; }
    jj = ((J % cny) + cny) % cny;
  }
  let k = jj * cnx + ii;
  if (maskC[u32(k)] != 0u) { return -1; }
  return k;
}

@compute @workgroup_size(8, 8)
fn mg_prolong(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= L.nx || id.y >= L.ny) { return; }
  let k = id.y * L.nx + id.x;
  if (mask[k] != 0u) { return; }
  let g = vec2f(f32(id.x), f32(id.y)) * 0.5 - 0.25; // fine cell centre in coarse index space
  let c0 = vec2i(floor(g));
  let t = g - vec2f(c0);
  var sum = 0.0;
  var weight = 0.0;
  for (var b = 0; b < 2; b++) {
    for (var a = 0; a < 2; a++) {
      let K = coarseCell(c0.x + a, c0.y + b);
      if (K >= 0) {
        let w = select(1.0 - t.x, t.x, a == 1) * select(1.0 - t.y, t.y, b == 1);
        sum += w * pc[u32(K)];
        weight += w;
      }
    }
  }
  if (weight > 0.0) { p[k] += sum / weight; }
}

src/lab/engine/shaders/pressure.wgsl

// Stage 4b · pressure: one red-black SOR sweep of  sum(p_nb - p) / h^2 = div(u*) / dt.
// Walls, inflow and obstacles are Neumann (no term); outflow edges are Dirichlet with ghost p = 0.
override COLOR: u32 = 0u;

@group(0) @binding(1) var<storage, read_write> pR: array<f32>;
@group(0) @binding(2) var<storage, read> divIn: array<f32>;
@group(0) @binding(3) var<storage, read> solid: array<u32>;

// Neighbour contribution: x = pressure, y = 1 when the neighbour enters the stencil.
fn neighbour(i: i32, j: i32, side: u32) -> vec2f {
  let outsideX = i < 0 || i >= i32(P.nx);
  let outsideY = j < 0 || j >= i32(P.ny);
  if ((outsideX && !periodicX()) || (outsideY && !periodicY())) {
    return vec2f(0.0, select(0.0, 1.0, side == EDGE_OUTFLOW));
  }
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); }
  let k = iC(u32(ii), u32(jj));
  if (solid[k] != 0u) { return vec2f(0.0, 0.0); }
  return vec2f(pR[k], 1.0);
}

@compute @workgroup_size(8, 8)
fn pressure(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= P.nx || id.y >= P.ny || ((id.x + id.y) & 1u) != COLOR) { return; }
  let k = iC(id.x, id.y);
  if (solid[k] != 0u) { return; }
  let i = i32(id.x);
  let j = i32(id.y);
  let s = neighbour(i - 1, j, P.edgeL) + neighbour(i + 1, j, P.edgeR)
        + neighbour(i, j - 1, P.edgeB) + neighbour(i, j + 1, P.edgeT);
  if (s.y == 0.0) { return; } // a cell sealed off on all sides has nothing to balance
  let gaussSeidel = (s.x - (P.h * P.h / P.dt) * divIn[k]) / s.y;
  pR[k] = mix(pR[k], gaussSeidel, P.omega);
}

src/lab/engine/shaders/gradient.wgsl

// Stage 4c · pressure push: u -= (dt/h) * grad p on every free face, making the field divergence-free.
@group(0) @binding(1) var<storage, read_write> uG: array<f32>;
@group(0) @binding(2) var<storage, read_write> vG: array<f32>;
@group(0) @binding(3) var<storage, read> pIn: array<f32>;
@group(0) @binding(4) var<storage, read> solid: array<u32>;

// Pressure of cell (i, j); beyond a non-periodic edge the ghost pressure is 0 (used by outflow faces only).
fn pAt(i: i32, j: i32) -> f32 {
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); } else if (i < 0 || i >= i32(P.nx)) { return 0.0; }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); } else if (j < 0 || j >= i32(P.ny)) { return 0.0; }
  return pIn[iC(u32(ii), u32(jj))];
}

fn uProjected(i: u32, j: u32) -> bool {
  if (!periodicX() && i == 0u) { return P.edgeL == EDGE_OUTFLOW && !isSolid(0, i32(j)); }
  if (!periodicX() && i == P.nx) { return P.edgeR == EDGE_OUTFLOW && !isSolid(i32(P.nx) - 1, i32(j)); }
  return !uFixed(i, j);
}

fn vProjected(i: u32, j: u32) -> bool {
  if (!periodicY() && j == 0u) { return P.edgeB == EDGE_OUTFLOW && !isSolid(i32(i), 0); }
  if (!periodicY() && j == P.ny) { return P.edgeT == EDGE_OUTFLOW && !isSolid(i32(i), i32(P.ny) - 1); }
  return !vFixed(i, j);
}

@compute @workgroup_size(8, 8)
fn gradient(@builtin(global_invocation_id) id: vec3u) {
  let i = id.x;
  let j = id.y;
  let s = P.dt / P.h;
  if (i <= P.nx && j < P.ny && uProjected(i, j)) {
    uG[iU(i, j)] -= s * (pAt(i32(i), i32(j)) - pAt(i32(i) - 1, i32(j)));
  }
  if (i < P.nx && j <= P.ny && vProjected(i, j)) {
    vG[iV(i, j)] -= s * (pAt(i32(i), i32(j)) - pAt(i32(i), i32(j) - 1));
  }
}

⑤ Boundaries

This stage sets the velocity on the faces at the edges of the experiment window and on the faces that touch obstacles. A face across a wall gets 0, so no water flows through the wall. An inlet face gets the inflow speed, an outlet face follows the face just inside it (except after stage ④, when it keeps the value the projection gave it), and a face touching an obstacle takes the obstacle’s velocity. This stage does not set the speed at which water slides along a wall: the other stages set it when they work out the values beyond the wall, according to the kind of wall. This stage runs after each of the other stages that run, so whatever another stage works out at the edges, the next stage always starts from a velocity whose boundaries are right.

boundary.wgsl

src/lab/engine/shaders/boundary.wgsl

// Stage 5 · boundaries: edge faces (walls, inflow, outflow, periodic copies) and faces touching obstacles.
// Tangential wall conditions are applied when sampling (ghost values in wgsl.ts), not stored here.
override SET_OUTFLOW: bool = true; // false after projection: outflow faces keep their projected values

@group(0) @binding(1) var<storage, read_write> uB: array<f32>;
@group(0) @binding(2) var<storage, read_write> vB: array<f32>;
@group(0) @binding(3) var<storage, read> solid: array<u32>;
@group(0) @binding(4) var<storage, read> solidVel: array<vec2f>;

fn solidVelAt(i: i32, j: i32) -> vec2f {
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); }
  return solidVel[iC(u32(ii), u32(jj))];
}

fn uTouches(i: i32, j: i32) -> bool { return isSolid(i - 1, j) || isSolid(i, j); }
fn vTouches(i: i32, j: i32) -> bool { return isSolid(i, j - 1) || isSolid(i, j); }

// Velocity of a face touched by an obstacle: mean of the adjacent solid cells.
fn uSolid(i: i32, j: i32) -> f32 {
  let a = isSolid(i - 1, j);
  let b = isSolid(i, j);
  if (a && b) { return 0.5 * (solidVelAt(i - 1, j).x + solidVelAt(i, j).x); }
  if (a) { return solidVelAt(i - 1, j).x; }
  return solidVelAt(i, j).x;
}

fn vSolid(i: i32, j: i32) -> f32 {
  let a = isSolid(i, j - 1);
  let b = isSolid(i, j);
  if (a && b) { return 0.5 * (solidVelAt(i, j - 1).y + solidVelAt(i, j).y); }
  if (a) { return solidVelAt(i, j - 1).y; }
  return solidVelAt(i, j).y;
}

// Final value of an interior face after this pass, computed instead of read when this pass rewrites it.
fn uFinal(i: i32, j: i32) -> f32 {
  if (uTouches(i, j)) { return uSolid(i, j); }
  return uB[iU(u32(i), u32(j))];
}

fn vFinal(i: i32, j: i32) -> f32 {
  if (vTouches(i, j)) { return vSolid(i, j); }
  return vB[iV(u32(i), u32(j))];
}

fn edgeValue(kind: u32, inflowValue: f32, inner: f32, current: f32) -> f32 {
  if (kind == EDGE_INFLOW) { return inflowValue; }
  if (kind == EDGE_OUTFLOW) { return select(current, inner, SET_OUTFLOW); }
  return 0.0; // no-slip and slip walls: no flow through the wall
}

@compute @workgroup_size(8, 8)
fn boundary(@builtin(global_invocation_id) id: vec3u) {
  let i = i32(id.x);
  let j = i32(id.y);
  let nx = i32(P.nx);
  let ny = i32(P.ny);
  if (i <= nx && j < ny) {
    let k = iU(id.x, id.y);
    if (periodicX() && i == nx) {
      uB[k] = uFinal(0, j);
    } else if (uTouches(i, j)) {
      uB[k] = uSolid(i, j); // obstacles win over edge rules, also on edge faces
    } else if (!periodicX() && i == 0) {
      uB[k] = edgeValue(P.edgeL, P.inflowL, uFinal(1, j), uB[k]);
    } else if (!periodicX() && i == nx) {
      uB[k] = edgeValue(P.edgeR, P.inflowR, uFinal(nx - 1, j), uB[k]);
    }
  }
  if (i < nx && j <= ny) {
    let k = iV(id.x, id.y);
    if (periodicY() && j == ny) {
      vB[k] = vFinal(i, 0);
    } else if (vTouches(i, j)) {
      vB[k] = vSolid(i, j);
    } else if (!periodicY() && j == 0) {
      vB[k] = edgeValue(P.edgeB, P.inflowB, vFinal(i, 1), vB[k]);
    } else if (!periodicY() && j == ny) {
      vB[k] = edgeValue(P.edgeT, P.inflowT, vFinal(i, ny - 1), vB[k]);
    }
  }
}

The code put in front of the shaders

The engine runs each shader file with shared code put in front of it. Most shaders start with COMMON: a structure holding values such as the grid’s size, one step’s time, the viscosity, the kinds of edges and the body force, and functions that work out the numbers of faces and cells. SOLID, which tells solid cells apart, and TRACE, the backward trace of stage ② advection, go in front of the shaders that need them. The multigrid shaders start with LEVEL_WGSL, the size of one level of the grid, instead of COMMON. Functions in wgsl.ts write in the rest to match the names of the values each shader uses: a function that estimates the velocity anywhere from the faces around it, one that estimates the ink, and the pressure equation on one level.

COMMON · SOLID · TRACE · LEVEL_WGSL

src/lab/engine/wgsl.ts · COMMON

struct Params {
  nx: u32, ny: u32, h: f32, dt: f32,
  nu: f32, omega: f32, edgeL: u32, edgeR: u32,
  edgeB: u32, edgeT: u32, inflowL: f32, inflowR: f32,
  inflowB: f32, inflowT: f32, gx: f32, gy: f32,
  splatCount: u32, pad0: u32, pad1: u32, pad2: u32,
};
@group(0) @binding(0) var<uniform> P: Params;

const EDGE_NOSLIP: u32 = 0u;
const EDGE_SLIP: u32 = 1u;
const EDGE_INFLOW: u32 = 2u;
const EDGE_OUTFLOW: u32 = 3u;
const EDGE_PERIODIC: u32 = 4u;

fn periodicX() -> bool { return P.edgeL == EDGE_PERIODIC; }
fn periodicY() -> bool { return P.edgeB == EDGE_PERIODIC; }
fn wrapi(i: i32, n: i32) -> i32 { return ((i % n) + n) % n; }

// u face (i, j): x = i*h, y = (j+0.5)*h | v face (i, j): x = (i+0.5)*h, y = j*h | cell (i, j): centre
fn iU(i: u32, j: u32) -> u32 { return j * (P.nx + 1u) + i; }
fn iV(i: u32, j: u32) -> u32 { return j * P.nx + i; }
fn iC(i: u32, j: u32) -> u32 { return j * P.nx + i; }

fn domainSize() -> vec2f { return vec2f(f32(P.nx), f32(P.ny)) * P.h; }

// Keeps a traced point inside the domain (wraps along periodic axes).
fn confine(p: vec2f) -> vec2f {
  let L = domainSize();
  var q = p;
  if (periodicX()) { q.x = q.x - floor(q.x / L.x) * L.x; } else { q.x = clamp(q.x, 0.0, L.x); }
  if (periodicY()) { q.y = q.y - floor(q.y / L.y) * L.y; } else { q.y = clamp(q.y, 0.0, L.y); }
  return q;
}

src/lab/engine/wgsl.ts · SOLID

fn isSolid(i: i32, j: i32) -> bool {
  var ii = i;
  var jj = j;
  if (periodicX()) { ii = wrapi(i, i32(P.nx)); }
  if (periodicY()) { jj = wrapi(j, i32(P.ny)); }
  if (ii < 0 || jj < 0 || ii >= i32(P.nx) || jj >= i32(P.ny)) { return false; }
  return solid[iC(u32(ii), u32(jj))] != 0u;
}

// A u face (i, j) lies between cells (i-1, j) and (i, j). Fixed faces are owned by the boundary kernel:
// edge faces, the duplicate face nx on a periodic axis, and every face touching an obstacle.
fn uFixed(i: u32, j: u32) -> bool {
  if (periodicX()) {
    if (i == P.nx) { return true; }
  } else if (i == 0u || i == P.nx) {
    return true;
  }
  return isSolid(i32(i) - 1, i32(j)) || isSolid(i32(i), i32(j));
}

fn vFixed(i: u32, j: u32) -> bool {
  if (periodicY()) {
    if (j == P.ny) { return true; }
  } else if (j == 0u || j == P.ny) {
    return true;
  }
  return isSolid(i32(i), i32(j) - 1) || isSolid(i32(i), i32(j));
}

src/lab/engine/wgsl.ts · TRACE

// dtSigned > 0 traces back in time (semi-Lagrangian departure point); dtSigned < 0 traces forward.
fn trace(p: vec2f, dtSigned: f32) -> vec2f {
  let mid = confine(p - 0.5 * dtSigned * velAt(p));
  return confine(p - dtSigned * velAt(mid));
}

src/lab/engine/wgsl.ts · LEVEL_WGSL

struct Level {
  nx: u32, ny: u32, cnx: u32, cny: u32,
  h: f32, fScale: f32, edgeL: u32, edgeR: u32,
  edgeB: u32, edgeT: u32, omega: f32, iterations: u32,
};
@group(0) @binding(0) var<uniform> L: Level;

const LEVEL_OUTFLOW: u32 = 3u;
const LEVEL_PERIODIC: u32 = 4u;

The whole equation again

Here are the equations this book built up, chapter by chapter.

DuDt=∂u∂t+(u⋅∇)u=−∇pρ+ν∇2u+f,∇⋅u=0\mat{\frac{D\mathbf{u}}{Dt}} = \frac{\partial \mathbf{u}}{\partial t} + \adv{(\mathbf{u}\cdot\nabla)\mathbf{u}} = \prs{-\frac{\nabla p}{\rho}} + \vis{\nu\nabla^2\mathbf{u}} + \frc{\mathbf{f}}, \qquad \cont{\nabla\cdot\mathbf{u} = 0}
  • u\vel{\mathbf{u}} is the flow’s velocity. You can follow a parcel of water or watch one point (chapter 1).
  • ∇⋅u=0\cont{\nabla\cdot\mathbf{u} = 0} is conservation of mass: what flows in must flow out. With the pressure off, the water crowds and thins (chapter 2).
  • −∇p/ρ\prs{-\nabla p/\rho} is the pressure, which adjusts everywhere at once so the water neither crowds nor thins, and pushes the water (chapter 3). It is also what holds up a wing.
  • ν∇2u\vis{\nu\nabla^2\mathbf{u}} is viscosity, which shares velocity evenly between neighbours. It makes a stirred swirl die away (chapter 4).
  • (u⋅∇)u\adv{(\mathbf{u}\cdot\nabla)\mathbf{u}} is the flow carrying its own velocity along. Switch it off, and a swirl stands where it is (chapter 5).
  • f\frc{\mathbf{f}} is the force from outside. Switch it off, and stirring moves no water (this chapter).

In one sentence

Each term of the Navier–Stokes equations does one job, the engine runs that work as stages in order every step, and even a flow that all the terms make together, such as the lift on a wing, reads as the work of those terms.

Something to think about

In the full lab, switching off only the pressure empties the space behind the cylinder and crowds the water in places. With the pressure term on, what prevents this?

Show answer

The pressure. Pressure adjusts everywhere at once so that the water neither crowds nor thins (chapter 3). In front of the cylinder the pressure rises and turns the oncoming water aside in advance; behind it the pressure drops and draws the water around into the gap. In the engine, stage ④, the projection, measures the divergence, solves for the pressure that removes it, and subtracts its gradient from the velocity. When the pressure went back on in this chapter’s lab, the divergence came back to within 0.000005 of zero.

Something to think about

This chapter’s wing is symmetric top to bottom. At an angle of attack of −8°, with the leading edge tipped down, which way does the lift point, and which part of a line that set off together reaches the trailing edge first?

Show answer

The flow past a symmetric wing at −8° is the flow at 8° turned upside down. So the flow below the wing is faster and its pressure lower, and the lift points down. The lift coefficient is as large as at 8°, 0.49, with the opposite sign. The water passing below reaches the trailing edge first, about 0.13 s ahead, and the flow behind the wing turns up. The wing’s shape is the same, and the lift points the other way. What decides lift is not the length of the surfaces but the way the wing turns the flow.