Gray-Scott Reaction-Diffusion in GLSL on 2D Grids

This article details the mathematical formulation of the Gray-Scott reaction-diffusion system and demonstrates how to implement its numerical evaluation on a two-dimensional discrete grid using OpenGL Shading Language (GLSL). It covers the governing differential equations, spatial discretization via a discrete Laplacian operator, and the final fragment shader code used for ping-pong texture computation.

Mathematical Formulation

The Gray-Scott model simulates the interaction of two virtual chemical substances, typically denoted as \(U\) and \(V\). Chemical \(U\) is continuously added at a feed rate \(F\), while \(V\) is converted and removed at a combined kill rate \(F + k\). Two units of \(V\) react with one unit of \(U\) to produce more \(V\) (\(U + 2V \rightarrow 3V\)).

The continuous partial differential equations governing their concentrations over time \(t\) are:

\[\frac{\partial U}{\partial t} = D_u \nabla^2 U - UV^2 + F(1 - U)\]

\[\frac{\partial V}{\partial t} = D_v \nabla^2 V + UV^2 - (F + k)V\]

Where:

\(U\) and \(V\) are the concentrations of the chemicals (\(0 \le U, V \le 1\)).
\(D_u\) and \(D_v\) are the diffusion rates of \(U\) and \(V\), respectively (usually \(D_u = 1.0\) and \(D_v = 0.5\)).
\(\nabla^2\) is the Laplacian operator, representing spatial diffusion across the 2D surface.
\(F\) is the feed rate parameter.
\(k\) is the kill rate parameter.

Grid Discretization

To evaluate this system on a 2D digital grid (such as a texture in a GPU pipeline), space and time must be discretized.

Using Euler's forward method, time stepping is computed as:

\[U_{t+\Delta t} = U_t + \left( D_u \nabla^2 U - UV^2 + F(1 - U) \right) \Delta t\]

\[V_{t+\Delta t} = V_t + \left( D_v \nabla^2 V + UV^2 - (F + k)V \right) \Delta t\]

The spatial Laplacian \(\nabla^2\) calculates the second spatial derivative by comparing the value of a target cell to its neighbors. A standard 9-point convolution stencil provides smooth, isotropic diffusion across a 2D texture grid:

Center cell: weight \(-1.0\)
Orthogonal neighbors (up, down, left, right): weight \(0.20\) each
Diagonal neighbors: weight \(0.05\) each

GLSL Fragment Shader Implementation

In GLSL, the simulation is executed as a fragment shader rendering to an off-screen framebuffer. A ping-pong buffer technique is required: one texture stores the current state \((U, V)\) in its red and green channels, while the shader writes the updated state to an alternate texture.

#version 330 core

out vec4 FragColor;

uniform sampler2D u_texture; // Previous state (R = U, G = V)
uniform vec2 u_resolution;   // Texture dimensions in pixels
uniform float u_dt;          // Time step (e.g., 1.0)
uniform float u_feed;        // Feed rate (F)
uniform float u_kill;        // Kill rate (k)
uniform float u_diffU;       // Diffusion rate Du (e.g., 1.0)
uniform float u_diffV;       // Diffusion rate Dv (e.g., 0.5)

void main()
{
    vec2 texel = 1.0 / u_resolution;
    vec2 uv = gl_FragCoord.xy * texel;

    // Sample current cell
    vec2 state = texture(u_texture, uv).rg;
    float u = state.r;
    float v = state.g;

    // 9-point Laplacian stencil
    vec2 laplacian = vec2(0.0);

    // Orthogonal neighbors (weight: 0.20)
    laplacian += texture(u_texture, uv + vec2( texel.x, 0.0)).rg * 0.20;
    laplacian += texture(u_texture, uv + vec2(-texel.x, 0.0)).rg * 0.20;
    laplacian += texture(u_texture, uv + vec2(0.0,  texel.y)).rg * 0.20;
    laplacian += texture(u_texture, uv + vec2(0.0, -texel.y)).rg * 0.20;

    // Diagonal neighbors (weight: 0.05)
    laplacian += texture(u_texture, uv + vec2( texel.x,  texel.y)).rg * 0.05;
    laplacian += texture(u_texture, uv + vec2(-texel.x,  texel.y)).rg * 0.05;
    laplacian += texture(u_texture, uv + vec2( texel.x, -texel.y)).rg * 0.05;
    laplacian += texture(u_texture, uv + vec2(-texel.x, -texel.y)).rg * 0.05;

    // Center cell (weight: -1.0)
    laplacian -= state * 1.0;

    // Reaction term UV^2
    float uvv = u * v * v;

    // Differential equations
    float du = (u_diffU * laplacian.r) - uvv + (u_feed * (1.0 - u));
    float dv = (u_diffV * laplacian.g) + uvv - ((u_feed + u_kill) * v);

    // Euler integration step
    float newU = clamp(u + du * u_dt, 0.0, 1.0);
    float newV = clamp(v + dv * u_dt, 0.0, 1.0);

    FragColor = vec4(newU, newV, 0.0, 1.0);
}

Execution Considerations

To achieve numerical stability and avoid artifacts:

Use 32-bit floating-point textures (GL_RG32F or GL_RGBA32F) to prevent rounding errors that cause the reaction to collapse prematurely.
Set texture wrapping to GL_REPEAT for seamless periodic boundary conditions, or GL_CLAMP_TO_EDGE for closed walls.
Keep the integration time step \(u\_dt \le 1.0\) when using the standard diffusion rates (\(D_u = 1.0, D_v = 0.5\)) to ensure the explicit numerical integration does not diverge.