Boundary Conditions in GLSL Eulerian Fluid Simulations

In Eulerian fluid simulations, fluid motion is simulated over a static grid where velocity and pressure fields are computed iteratively using GPU fragment shaders. Enforcing boundary conditions along solid obstacles is critical to preventing fluid from penetrating surfaces and ensuring physically plausible behavior like drag and friction. This article explains how solid obstacles are represented on the grid, how Dirichlet and Neumann boundary conditions are applied to velocity and pressure fields, and how to implement these constraints directly within GLSL shaders.

1. Obstacle Representation on the Grid

In a GLSL-based Eulerian simulation, the computational domain is mapped to 2D or 3D textures. Obstacles are typically represented using an obstacle mask texture (often encoded as a single-channel R8 or R16F texture) or an analytical signed distance field (SDF).

Discrete Mask: A value of 1.0 denotes a solid cell, while 0.0 denotes fluid.
Signed Distance Field (SDF): A continuous distance representation where values \(\le 0.0\) define the solid interior, allowing for sub-grid boundary accuracy.

During rendering passes, fragment shaders sample this obstacle texture to identify whether the current cell or its orthogonal neighbors are solid.

2. Velocity Boundary Conditions

Velocity at a solid boundary must prevent fluid from penetrating the obstacle while accounting for fluid-surface friction. Two primary boundary conditions are used:

Free-Slip Boundaries (Inviscid Fluids)

Fluid can slide freely along the surface, but cannot pass through it: \[\vec{u} \cdot \vec{n} = \vec{u}_{\text{solid}} \cdot \vec{n}\] The normal velocity component relative to the boundary is zero, while the tangential component is unaffected.

No-Slip Boundaries (Viscous Fluids)

Fluid adheres to the solid surface due to viscosity: \[\vec{u} = \vec{u}_{\text{solid}}\] For stationary obstacles, the velocity is set to zero (\(\vec{u} = 0\)).

GLSL Implementation for Velocity

In the advection step, if a cell is solid, its velocity is immediately forced to the obstacle's velocity. When advecting fluid near boundaries, back-tracing coordinates that land inside an obstacle should be clamped or reflected back into the fluid domain.

// Enforcing velocity boundary conditions in a fragment shader
uniform sampler2D u_velocityTexture;
uniform sampler2D u_obstacleTexture;
uniform vec2 u_texelSize;

void main() {
    vec2 uv = gl_FragCoord.xy * u_texelSize;
    float isSolid = texture(u_obstacleTexture, uv).r;

    if (isSolid > 0.5) {
        // Stationary obstacle: no velocity
        gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
        return;
    }

    vec2 vel = texture(u_velocityTexture, uv).xy;
    
    // Check neighbor obstacles (Left, Right, Bottom, Top)
    float obsL = texture(u_obstacleTexture, uv - vec2(u_texelSize.x, 0.0)).r;
    float obsR = texture(u_obstacleTexture, uv + vec2(u_texelSize.x, 0.0)).r;
    float obsB = texture(u_obstacleTexture, uv - vec2(0.0, u_texelSize.y)).r;
    float obsT = texture(u_obstacleTexture, uv + vec2(0.0, u_texelSize.y)).r;

    // Zero out velocity components pointing into solid cells
    if (obsL > 0.5 && vel.x < 0.0) vel.x = 0.0;
    if (obsR > 0.5 && vel.x > 0.0) vel.x = 0.0;
    if (obsB > 0.5 && vel.y < 0.0) vel.y = 0.0;
    if (obsT > 0.5 && vel.y > 0.0) vel.y = 0.0;

    gl_FragColor = vec4(vel, 0.0, 1.0);
}

3. Pressure Boundary Conditions (Poisson Equation)

To enforce incompressibility (\(\nabla \cdot \vec{u} = 0\)), the pressure Poisson equation \(\nabla^2 p = \nabla \cdot \vec{u}^*\) is solved iteratively, typically using Jacobi iterations.

Solid boundaries require a homogeneous Neumann boundary condition for pressure: \[\frac{\partial p}{\partial \vec{n}} = 0\]

This condition dictates that the pressure gradient normal to the boundary is zero, meaning fluid acceleration perpendicular to the wall is zero.

Modified Jacobi Stencil

In a standard Jacobi solver without obstacles, the center cell pressure is updated using the average of its four neighbors minus the divergence: \[p_{i,j}^{(k+1)} = \frac{p_{i+1,j} + p_{i-1,j} + p_{i,j+1} + p_{i,j-1} - \Delta x^2 \cdot (\nabla \cdot \vec{u})}{4}\]

When a neighbor cell is solid, the Neumann condition (\(\frac{\partial p}{\partial \vec{n}} = 0\)) means the pressure at the solid cell is treated as equal to the pressure of the fluid cell itself (\(p_{\text{solid}} = p_{\text{center}}\)).

// Modified Jacobi pressure iteration accounting for boundaries
uniform sampler2D u_pressureTexture;
uniform sampler2D u_divergenceTexture;
uniform sampler2D u_obstacleTexture;
uniform vec2 u_texelSize;

void main() {
    vec2 uv = gl_FragCoord.xy * u_texelSize;
    
    // Return early if the current cell is inside an obstacle
    if (texture(u_obstacleTexture, uv).r > 0.5) {
        gl_FragColor = vec4(0.0);
        return;
    }

    float pC = texture(u_pressureTexture, uv).r;
    float div = texture(u_divergenceTexture, uv).r;

    // Sample neighbors
    vec2 uvL = uv - vec2(u_texelSize.x, 0.0);
    vec2 uvR = uv + vec2(u_texelSize.x, 0.0);
    vec2 uvB = uv - vec2(0.0, u_texelSize.y);
    vec2 uvT = uv + vec2(0.0, u_texelSize.y);

    float obsL = texture(u_obstacleTexture, uvL).r;
    float obsR = texture(u_obstacleTexture, uvR).r;
    float obsB = texture(u_obstacleTexture, uvB).r;
    float obsT = texture(u_obstacleTexture, uvT).r;

    // If neighbor is solid, copy center pressure (dp/dn = 0)
    float pL = (obsL > 0.5) ? pC : texture(u_pressureTexture, uvL).r;
    float pR = (obsR > 0.5) ? pC : texture(u_pressureTexture, uvR).r;
    float pB = (obsB > 0.5) ? pC : texture(u_pressureTexture, uvB).r;
    float pT = (obsT > 0.5) ? pC : texture(u_pressureTexture, uvT).r;

    float newPressure = (pL + pR + pB + pT - div) * 0.25;
    gl_FragColor = vec4(newPressure, 0.0, 0.0, 1.0);
}

4. Projection Step

In the final step of the projection method, velocity is updated by subtracting the pressure gradient: \[\vec{u} = \vec{u}^* - \nabla p\]

When calculating \(\nabla p\) near solid boundaries, the boundary conditions must remain consistent:

If the neighbor along an axis is solid, the pressure gradient along that axis should evaluate to zero, preventing false accelerations into or out of the solid wall.
After the gradient subtraction is applied across all fluid cells, the velocity texture inside all solid cells is explicitly overwritten to zero (or the obstacle's velocity vector).