GLSL Jacobi Iteration for Fluid Pressure Poisson

This article details the discrete mathematical formula and GLSL implementation of the Jacobi iteration method for solving the pressure Poisson equation in grid-based fluid simulations. You will learn the mathematical derivation of the finite difference stencil, the standard iteration formula, and how to structure the fragment shader code for GPU execution.

The Pressure Poisson Equation

In incompressible Navier-Stokes fluid simulations, enforcing incompressibility requires projecting an intermediate, divergent velocity field \(\mathbf{u}^*\) into a divergence-free velocity field \(\mathbf{u}\). This projection yields the pressure Poisson equation:

\[\nabla^2 p = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^*\]

Assuming normalized density (\(\rho = 1\)) and grid spacing (\(\Delta x\)), the right-hand side is simplified to the divergence field scaled by the time step, often represented simply as \(b\):

\[\nabla^2 p = b\]

Finite Difference Discretization

Using a standard central difference approximation on a 2D Cartesian grid with cell size \(\Delta x\), the continuous Laplacian operator \(\nabla^2 p\) expands to:

\[\nabla^2 p \approx \frac{p_{i+1, j} + p_{i-1, j} + p_{i, j+1} + p_{i, j-1} - 4p_{i, j}}{(\Delta x)^2}\]

Setting this equal to the divergence term \(b_{i, j}\) yields:

\[\frac{p_{i+1, j} + p_{i-1, j} + p_{i, j+1} + p_{i, j-1} - 4p_{i, j}}{(\Delta x)^2} = b_{i, j}\]

The Jacobi Iteration Formula

The Jacobi method iteratively isolates the center cell pressure \(p_{i, j}\) at step \(k+1\) using the neighboring pressure values from the previous step \(k\):

\[p_{i, j}^{(k+1)} = \frac{p_{i-1, j}^{(k)} + p_{i+1, j}^{(k)} + p_{i, j-1}^{(k)} + p_{i, j+1}^{(k)} - (\Delta x)^2 b_{i, j}}{4}\]

In graphics programming literature (such as GPU Gems), this is generalized to solve the linear system \(Ax = b\) using the recurrence relation:

\[x_{i, j}^{(k+1)} = \frac{x_{i-1, j}^{(k)} + x_{i+1, j}^{(k)} + x_{i, j-1}^{(k)} + x_{i, j+1}^{(k)} + \alpha b_{i, j}}{\beta}\]

For the 2D pressure Poisson equation:

• \(\alpha = -(\Delta x)^2\)
• \(\beta = 4\)

For a 3D grid with 6 neighbors:

• \(\alpha = -(\Delta x)^2\)
• \(\beta = 6\)

GLSL Fragment Shader Implementation

Because GLSL cannot read from and write to the same texture simultaneously, the Jacobi solver executes across multiple render passes using two ping-pong framebuffer textures.

#version 300 es
precision highp float;

out vec4 fragColor;

uniform sampler2D u_pressure;   // Texture containing p^(k)
uniform sampler2D u_divergence; // Texture containing b (div(u*))

uniform vec2 u_texelSize;       // 1.0 / resolution
uniform float u_alpha;          // -(dx * dx)
uniform float u_beta;           // 4.0

in vec2 v_uv;

void main() {
    // Sample 4-point neighboring pressure values from step k
    float pL = texture(u_pressure, v_uv - vec2(u_texelSize.x, 0.0)).r;
    float pR = texture(u_pressure, v_uv + vec2(u_texelSize.x, 0.0)).r;
    float pB = texture(u_pressure, v_uv - vec2(0.0, u_texelSize.y)).r;
    float pT = texture(u_pressure, v_uv + vec2(0.0, u_texelSize.y)).r;

    // Sample divergence at current cell
    float bC = texture(u_divergence, v_uv).r;

    // Compute updated pressure p^(k+1)
    float pNew = (pL + pR + pB + pT + u_alpha * bC) / u_beta;

    fragColor = vec4(pNew, 0.0, 0.0, 1.0);
}

Practical Considerations

• Iteration Count: Jacobi iterations converge slowly. In real-time applications, 20 to 50 iterations per frame provide a balance between visual stability and performance.
• Boundary Conditions: Neumann boundary conditions (\(\frac{\partial p}{\partial n} = 0\)) are applied at walls. If a neighboring cell lies inside a solid or outside the boundary, its pressure value must be substituted with the center cell's pressure (\(p_{\text{neighbor}} = p_{i, j}\)) to ensure zero pressure gradient across the boundary.