WebGL Fluid Simulation: Pressure Solving

October 4, 2025•10 min read

Table of Contents

From Static to Dynamic
Mathematical Foundation
Pressure Equation
Jacobi Iteration Method
Implementing Core Physics Shaders
Adding Interaction with a Splat Shader
Implementing the Display Shader
Orchestrating the Complete Interactive Loop
Demo Implementation Notes
Summary

From Static to Dynamic

The previous article shown how to implement advection with a static velocity field. That animation missed what makes fluids feel real: the velocity field remained unchanged regardless of how much fluid moved through it. This article completes the physics simulation by implementing divergence correction and pressure solving, making the fluid incompressible and realistic. Additionally, it transforms the simulation from a passive animation into an interactive experience by implementing user input.

We will add velocity and density into the fluid simulation based on user mouse input and use the fluid's density to create a liquid painting effect.

The Mathematical Foundation

Making fluid incompressible and responsive to user input requires a grasp of the mathematical foundation behind these effects. The physics relies on key concepts from vector calculus.

Differential operators in vector calculus are the tools we use to describe how scalar and vector fields change in space. The most important ones for fluid simulation are gradient, divergence, curl, and the Laplacian.

Gradient (∇) The gradient operates on scalar fields and produces vector fields. If you have a scalar function f(x, y, z), the gradient ∇f points in the direction of steepest increase and has magnitude equal to the rate of change in that direction.

In Cartesian coordinates:

∇f = (∂f/∂x, ∂f/∂y, ∂f/∂z)

Think of it like this: if f represents temperature, the gradient tells you which way to move to get warmer fastest.

Divergence (∇·) Divergence operates on vector fields and produces scalar fields. For a vector field F = (Fx, Fy, Fz), the divergence measures how much the field is "spreading out" from a point:

∇·F = ∂Fx/∂x + ∂Fy/∂y + ∂Fz/∂z

Curl (∇×) Curl also operates on vector fields and produces another vector field. It measures the rotation or "twisting" of the field around a point:

∇×F = (∂Fz/∂y - ∂Fy/∂z, ∂Fx/∂z - ∂Fz/∂x, ∂Fy/∂x - ∂Fx/∂y)

If you imagine tiny paddle wheels placed in the field, curl tells you how fast they'd spin and in which direction.

The Laplacian (∇²) This is the divergence of the gradient: ∇²f = ∇·(∇f). It's crucial in physics, appearing in equations for heat flow, wave propagation, and quantum mechanics.

These operators are the computational tools that turn a simple advection animation into fluid that looks and behaves realistically. Each has a role in shaping how fluids move. Curl measures rotation in fluids and is important for vorticity effects, but this article focuses on incompressibility. To build up our intuition, let’s begin with divergence.

Understanding Divergence

Divergence is a measure of how much a vector field is expanding (positive divergence) or converging (negative divergence) at any given point. In fluid simulations, this vector field represents the velocity field of the fluid.

The connection to velocity comes through mass conservation:

If velocity has positive divergence, fluid flows outward from that point and density would decrease there
If velocity has negative divergence, fluid flows inward to that point and density would increase there

Real fluids are incompressible - water cannot be compressed or expanded under normal conditions. This means the velocity field must satisfy an important constraint: the divergence should be zero everywhere in regions where no sources or sinks should exist.

After each advection step, the velocity field develops non-zero divergence, breaking the incompressibility rule. This happens because the bilinear interpolation used in the semi-Lagrangian advection scheme to sample between grid points introduces interpolation errors. The pressure correction step fixes this by adjusting the velocities to restore the zero-divergence constraint, making the fluid behave like real, incompressible liquid.

The divergence of a 2D velocity field is calculated as:

∇ · v = (∂/∂x, ∂/∂y) · (v_x, v_y) = ∂v_x/∂x + ∂v_y/∂y

Here, ∇ = (∂/∂x, ∂/∂y) defines the nabla operator (or del operator) as a vector of partial derivative operators in 2D.

∇ · v means: take each component of the differential operator and apply it to the corresponding component of the vector field, then sum the results:

Apply ∂/∂x to v_x → get ∂v_x/∂x
Apply ∂/∂y to v_y → get ∂v_y/∂y
Sum them: ∂v_x/∂x + ∂v_y/∂y

Moving from continuous mathematics to GPU implementation requires discretizing these operators on a grid. On a discrete grid, we approximate these partial derivatives using finite differences. For a pixel at position (i, j), the divergence becomes:

div[i, j] = (v_x[i+1, j] - v_x[i-1, j]) / (2 × h) + (v_y[i, j+1] - v_y[i, j-1]) / (2 × h)

Where h is the grid spacing. In normalized texture coordinates [0, 1], this would be 1/texture_width. Since we work in texture space with unit spacing, where each grid cell has size 1 × 1, this simplifies to h = 1:

div[i, j] = (v_x[i+1, j] - v_x[i-1, j] + v_y[i, j+1] - v_y[i, j-1]) / 2

The goal is to find a pressure field p such that when we subtract its gradient from the velocity field, the result has zero divergence everywhere.

Understanding Gradient

The gradient of a scalar field indicates the direction of steepest increase and the rate of change at any given point. In fluid simulations, the gradient of the pressure field represents the forces that drive fluid motion: fluid naturally flows from high-pressure regions toward low-pressure regions.

For a 2D pressure field p, the gradient is a vector field:

∇p = (∂p/∂x, ∂p/∂y)

At each point, this vector points toward higher pressure regions. Its magnitude shows how rapidly pressure increases in that direction. On a discrete grid, the components are approximated using finite differences:

∇p[i, j].x = (p[i+1, j] - p[i-1, j]) / 2 ∇p[i, j].y = (p[i, j+1] - p[i, j-1]) / 2

The pressure correction step subtracts this gradient from the velocity field:

v_corrected = v_original - ∇p

This pushes fluid away from high-pressure regions toward low-pressure regions. The challenge is finding the pressure field whose gradient effectively cancels the divergence errors.

Pressure Equation

In fluid simulation, pressure isn’t the physical pressure in a pipe. It’s a mathematical field representing the "correction force" needed to make the fluid incompressible. Think of it as a map showing how strongly to push or pull the velocity at each point.

The relationship comes from a fundamental principle: subtracting the gradient of pressure from the velocity field produces a result with zero divergence everywhere. In other words, the gradient of the pressure field determines how to adjust the velocity so the fluid behaves incompressibly.

To find such a pressure field, we require a condition: the corrected velocity must have zero divergence:

∇·(v - ∇p) = 0

Expanding this using the linearity of the divergence operator:

∇·v - ∇·∇p = 0

Writing this mathematically leads directly to Poisson’s equation, where ∇·∇ is the Laplacian (∇²):

∇²p = ∇·v

Poisson's equation is a partial differential equation that appears throughout physics when one field (pressure) controls the sources and sinks of another field (velocity). This equation means: "find the pressure field that, when its gradient is subtracted from velocity, eliminates all the divergence."

The left side, ∇²p, is the Laplacian of pressure. The Laplacian measures how much a value differs from the average of its neighbors - it captures the "curvature" or "bending" of the field:

∇²p = ∂²p/∂x² + ∂²p/∂y²

On a discrete grid, the Laplacian becomes a simple average of the four neighboring values minus four times the center value:

p[i+1, j] + p[i-1, j] + p[i, j+1] + p[i, j-1] - 4×p[i, j] = div[i, j]

This discrete form says: "the pressure at this point should equal the average of its neighbors, adjusted by the divergence error that needs fixing".

Solving this equation gives us the pressure field needed to correct the velocities.

Jacobi Iteration Method

Because a fragment shader can only write to its assigned pixel and cannot modify neighboring pixels during the same rendering pass, it's impossible to solve the Poisson equation for the entire grid in a single step. The Poisson equation creates a system where each pixel's correct pressure value depends on the correct pressure values of its neighbors - but those neighbors also depend on their neighbors, creating circular dependencies throughout the grid.

We use an iterative solver called the Jacobi iteration method, where each pixel updates its pressure based on the neighbor values from the previous iteration. By repeatedly applying this update across the entire grid, the solution gradually converges:

p[i, j] = (p[i + 1, j] + p[i - 1, j] + p[i, j + 1] + p[i, j - 1] - div[i, j]) * 0.25

Where:

The factor 0.25 computes the average of the four neighboring pressure values.

The Jacobi method requires its own ping-pong setup since we need to read from the pressure field and simultaneously write to it.

Implementing Core Physics Shaders

The physics calculations are handled by three main shaders that implement the concepts we've discussed.

Divergence Shader

The divergence shader calculates the divergence of the velocity field and stores it in a texture for use by the pressure solver.

glsl
1// GLSL - Divergence Shader
2#version 300 es
3precision highp float;
4// ... (in, out, uniforms) ...
5uniform sampler2D uVelocity;
6uniform vec2 uTexelSize;
7
8void main() {
9    vec2 vL = texture(uVelocity, vUV - vec2(uTexelSize.x, 0.0)).xy;
10    vec2 vR = texture(uVelocity, vUV + vec2(uTexelSize.x, 0.0)).xy;
11    vec2 vB = texture(uVelocity, vUV - vec2(0.0, uTexelSize.y)).xy;
12    vec2 vT = texture(uVelocity, vUV + vec2(0.0, uTexelSize.y)).xy;
13
14    float divergence = 0.5 * ((vR.x - vL.x) + (vT.y - vB.y));
15
16    outColor = vec4(divergence, 0.0, 0.0, 1.0);
17}

Pressure Shader

The pressure shader performs one iteration of the Jacobi method to solve Poisson's equation.

glsl
1// GLSL - Pressure Jacobi Iteration Shader
2#version 300 es
3precision highp float;
4// ... (in, out, uniforms) ...
5uniform sampler2D uPressure;
6uniform sampler2D uDivergence;
7uniform vec2 uTexelSize;
8
9void main() {
10    float pL = texture(uPressure, vUV - vec2(uTexelSize.x, 0.0)).r;
11    float pR = texture(uPressure, vUV + vec2(uTexelSize.x, 0.0)).r;
12    float pB = texture(uPressure, vUV - vec2(0.0, uTexelSize.y)).r;
13    float pT = texture(uPressure, vUV + vec2(0.0, uTexelSize.y)).r;
14
15    float divergence = texture(uDivergence, vUV).r;
16
17    float newPressure = (pL + pR + pB + pT - divergence) * 0.25;
18
19    outColor = vec4(newPressure, 0.0, 0.0, 1.0);
20}

This shader needs to run multiple times per frame (typically 20-50 iterations) with its own ping-pong buffer pair to converge to the correct pressure solution.

Gradient Subtraction Shader

Once we have the pressure field, we subtract its gradient from the velocity field to create a divergence-free result.

glsl
1// GLSL - Gradient Subtraction Shader
2#version 300 es
3precision highp float;
4
5in vec2 vUV;
6out vec4 outColor;
7
8uniform sampler2D uVelocity;
9uniform sampler2D uPressure;
10uniform vec2 uTexelSize;
11
12void main() {
13    // Current velocity
14    vec2 velocity = texture(uVelocity, vUV).xy;
15
16    // Sample neighboring pressure values
17    float pL = texture(uPressure, vUV - vec2(uTexelSize.x, 0.0)).r;
18    float pR = texture(uPressure, vUV + vec2(uTexelSize.x, 0.0)).r;
19    float pB = texture(uPressure, vUV - vec2(0.0, uTexelSize.y)).r;
20    float pT = texture(uPressure, vUV + vec2(0.0, uTexelSize.y)).r;
21
22    // Calculate pressure gradient
23    vec2 gradient = 0.5 * vec2(pR - pL, pT - pB);
24
25    // Subtract gradient from velocity
26    velocity -= gradient;
27
28    outColor = vec4(velocity, 0.0, 1.0);
29}

This final step ensures the velocity field has zero divergence, satisfying the incompressibility constraint.

Adding Interaction with a Splat Shader

To make the simulation interactive, we need a way to add velocity and dye in response to user input. We achieve this with a "splat" shader. This shader takes a center point (the cursor position), a color (for dye) or vector (for velocity), and a radius, and draws a smooth, Gaussian-like splat into a target texture.

glsl
1// GLSL - Splat Shader
2#version 300 es
3precision highp float;
4
5in vec2 vUV;
6out vec4 outColor;
7
8uniform sampler2D uTarget;
9uniform float uAspectRatio;
10uniform vec2 uPoint;
11uniform vec3 uColor;
12uniform float uRadius;
13
14void main() {
15    vec2 p = vUV - uPoint;
16    p.x *= uAspectRatio;
17
18    float sqDist = dot(p, p);
19    float radiusSq = uRadius * uRadius;
20
21    // Use a Gaussian falloff for a smooth splat
22    vec3 splat = exp(-sqDist / radiusSq) * uColor;
23
24    vec3 base = texture(uTarget, vUV).rgb;
25
26    // Add the new splat to the existing texture values
27    outColor = vec4(base + splat, 1.0);
28}

When the user drags their cursor, we call this shader twice: once to add a colored splat to the dye texture, and once to add a directional velocity splat to the velocity texture.

Implementing Display Shader

After computing the velocity field corrections and advecting the dye through the fluid, the final step is rendering the results to the screen. The display shader samples the dye texture - which now contains colors that have been transported and mixed by the realistic fluid motion - and determines how to visualize this data. For a dynamic painting effect, the shader renders the dye colors directly, creating fluid brushstrokes that respond to user interaction:

glsl
1// GLSL - Display Shader for Dynamic Painting
2#version 300 es
3precision highp float;
4
5in vec2 vUV;
6out vec4 outColor;
7
8uniform sampler2D uDisplayTexture; // The dye texture
9
10void main() {
11    vec3 color = texture(uDisplayTexture, vUV).rgb;
12
13    // Apply slight tone mapping for better visuals
14    color = pow(color, vec3(0.9));
15
16    outColor = vec4(color, 1.0);
17}

This approach treats the fluid as liquid paint, where the dye colors flow and swirl according to the physics simulation, creating organic painting effects as the user interacts with the canvas.

Variations: Liquid Reveal Mask

Instead of simply drawing the colored dye, you can use it to create a visual effect. This shader turns the fluid into a mask that "reveals" content behind the canvas. It renders a black overlay, but uses the brightness of the dye to control the alpha channel. Where there is no dye, the overlay is opaque black. Where the dye is bright, the overlay becomes transparent.

glsl
1void main() {
2    vec3 dyeColor = texture(uDisplayTexture, vUV).rgb;
3    float brightness = length(dyeColor);
4    float alpha = 1.0 - clamp(brightness * 1.5, 0.0, 1.0);
5    outColor = vec4(0.0, 0.0, 0.0, alpha); // Black overlay with variable transparency
6}

This version requires creating the WebGL context with { alpha: true } and enabling alpha blending: gl.enable(gl.BLEND); gl.blendFunc(gl.SRC_ALPHA, gl.ONE_MINUS_SRC_ALPHA);

Orchestrating the Complete Interactive Loop

Now that we have all the shader components, we can orchestrate them in the render loop. The full render loop here is a multi-step process that must be executed in the correct order every frame.

We first advect the velocity through itself because momentum must also move with the flow. This step transports existing motion - like swirls and jets - to their new positions. Only after the velocity field has moved do we apply the pressure correction to ensure the result is incompressible.

javascript
1function render() {
2    // Check for and apply user input splats
3    if (pointer.isDown) {
4        applySplat(pointer.x, pointer.y, pointer.dx, pointer.dy, ...);
5    }
6
7    // --- Step 1: Advect the velocity field by itself ---
8    // This makes forces like vortices and jets move through the grid
9    gl.bindFramebuffer(gl.FRAMEBUFFER, velocityFboPair.write.framebuffer);
10    gl.useProgram(advectionProgram);
11    // ... bind velocity texture to both quantity and velocity uniforms
12    gl.drawElements(...);
13    velocityFboPair.swap();
14
15    // --- Step 2: Calculate divergence of the advected velocity ---
16    gl.bindFramebuffer(gl.FRAMEBUFFER, divergenceFbo.framebuffer);
17    gl.useProgram(divergenceProgram);
18    // ... bind velocity
19    gl.drawElements(...);
20
21    // --- Step 3: Solve for pressure using Jacobi iterations ---
22    gl.useProgram(pressureProgram);
23    for (let i = 0; i < pressureIterations; i++) {
24        gl.bindFramebuffer(gl.FRAMEBUFFER, pressureFboPair.write.framebuffer);
25        // ... bind pressure and divergence textures
26        gl.drawElements(...);
27        pressureFboPair.swap();
28    }
29
30    // --- Step 4: Subtract pressure gradient to correct velocity ---
31    gl.bindFramebuffer(gl.FRAMEBUFFER, velocityFboPair.write.framebuffer);
32    gl.useProgram(gradientProgram);
33    // ... bind pressure and velocity textures
34    gl.drawElements(...);
35    velocityFboPair.swap();
36
37    // --- Step 5: Advect the dye using the corrected velocity field ---
38    gl.bindFramebuffer(gl.FRAMEBUFFER, dyeFboPair.write.framebuffer);
39    gl.useProgram(advectionProgram);
40    // ... bind velocity and dye textures
41    gl.drawElements(...);
42    dyeFboPair.swap();
43
44    // --- Step 6: Render the dye texture to the screen ---
45    gl.enable(gl.BLEND);
46    gl.bindFramebuffer(gl.FRAMEBUFFER, null);
47    gl.useProgram(displayProgram);
48    // ... bind final dye texture
49    gl.drawElements(...);
50
51    requestAnimationFrame(render);
52}

This extended render loop creates a proper fluid simulation where the velocity field reacts to the fluid movement, creating realistic swirling and conservation of mass.

Demo Implementation Notes

Demo is available here. Demo code is available here

Performance: The pressure solver is the most expensive part of the simulation. Start with 20 iterations and adjust based on performance needs. More iterations create more accurate pressure fields but take longer to compute.
Resolution: The demo uses a lower resolution for the physics simulation (velocity, pressure) and a higher resolution for the visible dye. This is a common optimization, as the dye contains high-frequency details and the velocity is smoothed through the pressure projection.
Boundary Conditions: The shaders use gl.CLAMP_TO_EDGE texture wrapping, which creates "no-slip" boundaries where fluid appears to stick to the walls.

Summary

This article completed the core mathematical foundation of fluid simulation by implementing divergence calculation and pressure correction. The simulation now supports the fundamental physics of fluid flow and creates realistic flow patterns.

From here, you can explore more advanced visual effects like vorticity confinement, add obstacles, or incorporate more complex forces to create even richer experiences.