Deriving 3D Curl Noise from Vector Fields in GLSL

This article explains how to derive 3D curl noise from analytical potential vector fields within OpenGL Shading Language (GLSL). It covers the underlying mathematical principles of the curl operator, explains why the resulting vector field is naturally divergence-free, and demonstrates how to implement this technique efficiently in a fragment or compute shader using both analytical approaches and finite differences.

The Mathematical Principle

Curl noise generates divergence-free velocity fields, which are critical for simulating fluid-like motion, smoke, and turbulence without unnatural compression or expansion artifacts.

The process begins with a 3D vector potential field, denoted as \(\vec{\Psi}(x, y, z) = (\psi_x, \psi_y, \psi_z)\). The velocity field \(\vec{v}\) is obtained by taking the mathematical curl of \(\vec{\Psi}\):

\[\vec{v} = \nabla \times \vec{\Psi}\]

According to vector calculus identities, the divergence of the curl of any twice continuously differentiable vector field is identically zero:

\[\nabla \cdot (\nabla \times \vec{\Psi}) = 0\]

Because the divergence is zero, the resulting field is incompressible.

Expanded into Cartesian coordinates, the individual velocity components are computed as:

\[v_x = \frac{\partial \psi_z}{\partial y} - \frac{\partial \psi_y}{\partial z}\]

\[v_y = \frac{\partial \psi_x}{\partial z} - \frac{\partial \psi_z}{\partial x}\]

\[v_z = \frac{\partial \psi_y}{\partial x} - \frac{\partial \psi_x}{\partial y}\]

Defining the Analytical Potential Field

To generate varied and organic patterns, the vector potential field \(\vec{\Psi}(p)\) can be constructed by combining three independent scalar potential fields. Common choices include analytic trigonometric functions, Perlin noise, or Simplex noise:

vec3 potentialField(vec3 p) {
    // Example using three offset scalar noise calls or analytical functions
    float psiX = snoise(p + vec3(0.0, 0.0, 0.0));
    float psiY = snoise(p + vec3(43.12, 17.56, 11.23));
    float psiZ = snoise(p + vec3(123.4, 78.9, 45.6));
    return vec3(psiX, psiY, psiZ);
}

If the scalar fields are completely analytical (e.g., sums of sine and cosine waves), you can calculate the exact partial derivatives analytically using the chain rule. However, procedural noise algorithms are typically differentiated numerically.

Computing Derivatives in GLSL

When analytical derivatives of the underlying noise functions are complex or unavailable, the partial derivatives can be computed using central finite differences. Central differences provide second-order accuracy (\(O(\epsilon^2)\)) and avoid directional bias.

Given a small sampling step \(\epsilon\):

\[\frac{\partial \psi}{\partial x} \approx \frac{\psi(x + \epsilon, y, z) - \psi(x - \epsilon, y, z)}{2\epsilon}\]

GLSL Implementation

Below is a complete GLSL function that samples a potential field and computes the 3D curl:

vec3 computeCurl(vec3 p, float epsilon) {
    vec2 e = vec2(epsilon, 0.0);

    // Partial derivatives with respect to X
    vec3 dPsi_dx = (potentialField(p + e.xyy) - potentialField(p - e.xyy)) / (2.0 * epsilon);

    // Partial derivatives with respect to Y
    vec3 dPsi_dy = (potentialField(p + e.yxy) - potentialField(p - e.yxy)) / (2.0 * epsilon);

    // Partial derivatives with respect to Z
    vec3 dPsi_dz = (potentialField(p + e.yyx) - potentialField(p - e.yyx)) / (2.0 * epsilon);

    // Compute curl components:
    // v.x = dPsi_z / dy - dPsi_y / dz
    // v.y = dPsi_x / dz - dPsi_z / dx
    // v.z = dPsi_y / dx - dPsi_x / dy
    vec3 curl;
    curl.x = dPsi_dy.z - dPsi_dz.y;
    curl.y = dPsi_dz.x - dPsi_dx.z;
    curl.z = dPsi_dx.y - dPsi_dy.x;

    return curl;
}

Performance Considerations

1. Sampling Overhead: Using central differences requires 6 evaluations of the vector potential field (or 18 individual scalar noise evaluations). If performance is constrained, forward differences require only 3 vector samples (plus the evaluation at \(p\)), though this sacrifices numerical symmetry.
2. Step Size Selection: The value of \(\epsilon\) dictates the fidelity of the derivatives. If \(\epsilon\) is too small, floating-point precision issues arise; if it is too large, fine turbulent details are smoothed out. A value between 1e-4 and 1e-2 relative to the noise frequency typically yields the best balance.