Vorticity Confinement in GLSL Fluid Solvers

Vorticity confinement is the standard algorithmic technique used in real-time GLSL fluid simulations to counteract numerical dissipation and regenerate turbulent micro-details. In grid-based Eulerian solvers, semi-Lagrangian advection frequently smooths out small eddies, causing smoke, water, and fire to appear overly viscous. By identifying the local centers of rotation and injecting artificial momentum back into the velocity field, vorticity confinement restores high-frequency rotational energy, ensuring fluids retain their chaotic, natural appearance without requiring denser simulation grids.

The Algorithm: Steinhoff and Underhill's Vorticity Confinement

Originally developed by John Steinhoff and David Underhill for aerodynamics and later introduced to computer graphics by Ronald Fedkiw, Jos Stam, and Henrik Wann Jensen, the vorticity confinement algorithm calculates a non-physical body force that pushes energy back into swirling regions.

The algorithm operates in four distinct mathematical steps across the discrete simulation grid:

1. Compute the Vorticity Field

The algorithm first evaluates the curl of the current velocity field \(\mathbf{u}\), which yields the vorticity vector \(\omega\):

\[\omega = \nabla \times \mathbf{u}\]

In a 2D GLSL simulation, the curl simplifies to a scalar field representing rotation along the perpendicular z-axis:

\[\omega = \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y}\]

In GLSL fragment shaders, this is evaluated using central finite differences by sampling neighboring texels (left, right, top, bottom) of the velocity texture.

2. Determine the Vorticity Gradients

Next, the solver computes the scalar magnitude of the vorticity, \(|\omega|\). It then computes the gradient vector \(\mathbf{\eta}\) of this magnitude:

\[\mathbf{\eta} = \nabla |\omega|\]

This gradient points directly toward the local centers of the vortices, where rotation is most intense.

3. Normalize the Direction Vectors

To isolate direction from magnitude, \(\mathbf{\eta}\) is normalized to produce unit vector \(\mathbf{N}\):

\[\mathbf{N} = \frac{\mathbf{\eta}}{\|\mathbf{\eta}\|}\]

A small epsilon value is added to the denominator in shader code to prevent division-by-zero errors in areas with zero or uniform vorticity.

4. Apply the Confinement Force

The confinement force \(\mathbf{f}_{vc}\) is calculated as the cross product of the unit vector \(\mathbf{N}\) and the original vorticity vector \(\omega\), scaled by a user-defined strength parameter \(\epsilon\) and grid spacing \(h\):

\[\mathbf{f}_{vc} = \epsilon h (\mathbf{N} \times \omega)\]

This generated force vector acts perpendicular to both the vortex direction and the gradient pointing to the vortex center. When added to the fluid's velocity field across time step \(\Delta t\), it accelerates the rotational flow without expanding or compressing the fluid volume, successfully counteracting numerical damping while preserving the divergence-free nature required by incompressible Navier-Stokes solvers.