Gray-Scott Reaction Diffusion on the GPU

1. The Gray-Scott Equations

The Gray-Scott model (1983, by Peter Gray and Stephen Scott) describes a continuously stirred tank reaction between two chemical species U and V:

• U is the "food" chemical — replenished from outside at feed rate F
• V is the "activator" — produced by the reaction U + 2V → 3V, killed at rate (F+k)

The three terms in the U equation are:

• D_u·∇²u — U diffuses across the domain (spreads, smooths)
• −u·v² — consumed by autocatalytic reaction with V
• F·(1−u) — continuously fed in; acts as a saturating source (faster when u is low)

2. Discrete Laplacian on a Grid

On a 2D grid with spacing h = 1 pixel, the standard 5-point stencil Laplacian is:

The 9-point stencil produces more isotropic patterns — spots and stripes without visible grid-axis bias. Use this in the shader.

3. Pattern Parameter Space

The F-k parameter space is a zoo of qualitatively different patterns. Classic named regions:

4. GPU Architecture: Ping-Pong Framebuffers

The problem: each step of the simulation reads the current state and writes the next state. You can't read from and write to the same texture simultaneously on a GPU.

The solution: ping-pong — allocate two framebuffers (A and B). Each step:

• Read from texture A, write to framebuffer B (the simulation shader)
• Next step: read from B, write to A
• Swap references and repeat

// Ping-pong state
let frontBuffer = texA;  // GPU reads from this
let backFBO     = fboB;  // GPU writes to this

function swap() {
  [frontBuffer, backFBO] = [backBuffer.texture, frontFBO];
  // After swap: front ← what was just written; back ← old buffer to overwrite next step
}

The render pass is separate: a full-screen quad that reads from frontBuffer and maps concentration values to colours for display. It does NOT write to the simulation textures.

5. WebGL 2 Setup

// Create two pairs of (texture + framebuffer) for ping-pong
function createPingPongBuffers(gl, width, height) {
  function createTex() {
    const tex = gl.createTexture();
    gl.bindTexture(gl.TEXTURE_2D, tex);
    gl.texImage2D(gl.TEXTURE_2D, 0, gl.RG32F, width, height, 0,
                  gl.RG, gl.FLOAT, null);  // 2-channel float: [u, v]
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.REPEAT);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.REPEAT);
    return tex;
  }
  function createFBO(tex) {
    const fbo = gl.createFramebuffer();
    gl.bindFramebuffer(gl.FRAMEBUFFER, fbo);
    gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0,
                            gl.TEXTURE_2D, tex, 0);
    gl.bindFramebuffer(gl.FRAMEBUFFER, null);
    return fbo;
  }
  const texA = createTex(), texB = createTex();
  return {
    texA, texB,
    fboA: createFBO(texA), fboB: createFBO(texB),
    front: texA, back: { fbo: createFBO(texB), tex: texB },
  };
}

6. Simulation Fragment Shader

#version 300 es
precision highp float;

uniform sampler2D uState;   // RG = [u, v]
uniform vec2  uResolution;
uniform float uF;           // feed rate
uniform float uK;           // kill rate
uniform float uDu;          // diffusion U
uniform float uDv;          // diffusion V
uniform float uDt;          // timestep

out vec4 outColor;

// 9-point isotropic Laplacian
vec2 laplacian(vec2 pos) {
  vec2 d = 1.0 / uResolution;
  vec2 c  = texture(uState, pos).rg;
  vec2 N  = texture(uState, pos + vec2( 0, d.y)).rg;
  vec2 S  = texture(uState, pos + vec2( 0,-d.y)).rg;
  vec2 E  = texture(uState, pos + vec2( d.x, 0)).rg;
  vec2 W  = texture(uState, pos + vec2(-d.x, 0)).rg;
  vec2 NE = texture(uState, pos + vec2( d.x, d.y)).rg;
  vec2 NW = texture(uState, pos + vec2(-d.x, d.y)).rg;
  vec2 SE = texture(uState, pos + vec2( d.x,-d.y)).rg;
  vec2 SW = texture(uState, pos + vec2(-d.x,-d.y)).rg;
  return 0.2*(N+S+E+W) + 0.05*(NE+NW+SE+SW) - 1.0*c;
}

void main() {
  vec2 uv = gl_FragCoord.xy / uResolution;
  vec2 state = texture(uState, uv).rg;
  float u = state.r;
  float v = state.g;

  vec2 lap = laplacian(uv);
  float reaction = u * v * v;

  float du = uDu * lap.r - reaction + uF * (1.0 - u);
  float dv = uDv * lap.g + reaction - (uF + uK) * v;

  float newU = clamp(u + du * uDt, 0.0, 1.0);
  float newV = clamp(v + dv * uDt, 0.0, 1.0);

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

7. Colour Mapping Shader

#version 300 es
precision highp float;
uniform sampler2D uState;
out vec4 fragColor;

// Inferno-inspired colour ramp for V concentration
vec3 inferno(float t) {
  const vec3 c0 = vec3(0.0002, 0.0016, 0.0139);
  const vec3 c1 = vec3(0.1140, 0.0631, 0.4149);
  const vec3 c2 = vec3(0.5419, 0.1773, 0.3773);
  const vec3 c3 = vec3(0.8785, 0.3835, 0.1551);
  const vec3 c4 = vec3(0.9891, 0.7945, 0.2824);
  t = clamp(t, 0.0, 1.0);
  if (t < 0.25) return mix(c0, c1, t * 4.0);
  if (t < 0.50) return mix(c1, c2, (t - 0.25) * 4.0);
  if (t < 0.75) return mix(c2, c3, (t - 0.50) * 4.0);
  return mix(c3, c4, (t - 0.75) * 4.0);
}

void main() {
  vec2 uv = gl_FragCoord.xy / vec2(textureSize(uState, 0));
  float v = texture(uState, uv).g;
  fragColor = vec4(inferno(v), 1.0);
}

8. Full JavaScript Driver

const canvas = document.getElementById('canvas');
const gl = canvas.getContext('webgl2');
gl.getExtension('EXT_color_buffer_float');

const W = 512, H = 512;
canvas.width  = W;
canvas.height = H;

// --- Compile shaders (helpers omitted for brevity) ---
const simProgram    = compileProgram(gl, simVertSrc, simFragSrc);
const renderProgram = compileProgram(gl, renderVertSrc, renderFragSrc);

// --- Full-screen quad VAO ---
const quadVAO = makeQuadVAO(gl);

// --- Ping-pong buffers ---
let { texA, texB, fboA, fboB } = createPingPongBuffers(gl, W, H);
let front = texA, frontFBO = fboA;
let back  = texB, backFBO  = fboB;

function swap() {
  [front, back]       = [back,  front];
  [frontFBO, backFBO] = [backFBO, frontFBO];
}

// --- Seed initial state ---
const seed = new Float32Array(W * H * 2).fill(0);
for (let i = 0; i < W * H; i++) seed[i * 2] = 1.0; // u=1 everywhere
for (let cy = 220; cy < 260; cy++) {          // seed patch
  for (let cx = 220; cx < 260; cx++) {
    const i = cy * W + cx;
    seed[i * 2]     = 0.50;
    seed[i * 2 + 1] = 0.25 + 0.01 * (Math.random() - 0.5);
  }
}
uploadTexture(gl, front, W, H, seed);

// --- Simulation parameters ---
const params = { F: 0.035, k: 0.065, Du: 0.16, Dv: 0.08, dt: 1.0 };
const STEPS_PER_FRAME = 8;

function simStep() {
  gl.useProgram(simProgram);
  gl.viewport(0, 0, W, H);
  gl.bindFramebuffer(gl.FRAMEBUFFER, backFBO);
  gl.activeTexture(gl.TEXTURE0);
  gl.bindTexture(gl.TEXTURE_2D, front);
  gl.uniform1i(gl.getUniformLocation(simProgram, 'uState'), 0);
  gl.uniform2f(gl.getUniformLocation(simProgram, 'uResolution'), W, H);
  gl.uniform1f(gl.getUniformLocation(simProgram, 'uF'), params.F);
  gl.uniform1f(gl.getUniformLocation(simProgram, 'uK'), params.k);
  gl.uniform1f(gl.getUniformLocation(simProgram, 'uDu'), params.Du);
  gl.uniform1f(gl.getUniformLocation(simProgram, 'uDv'), params.Dv);
  gl.uniform1f(gl.getUniformLocation(simProgram, 'uDt'), params.dt);
  drawQuad(gl, quadVAO);
  swap();
}

function render() {
  gl.useProgram(renderProgram);
  gl.viewport(0, 0, canvas.width, canvas.height);
  gl.bindFramebuffer(gl.FRAMEBUFFER, null);
  gl.activeTexture(gl.TEXTURE0);
  gl.bindTexture(gl.TEXTURE_2D, front);
  gl.uniform1i(gl.getUniformLocation(renderProgram, 'uState'), 0);
  drawQuad(gl, quadVAO);
}

function loop() {
  for (let i = 0; i < STEPS_PER_FRAME; i++) simStep();
  render();
  requestAnimationFrame(loop);
}
loop();

Next Steps

• Add interactive seeding on click/drag — draw V=(0.25) patches where the mouse is
• Vary F and k spatially — create a gradient map to get different pattern types coexisting
• Try 3D → slice output at various Z planes for volumetric patterns
• Port to an 8-bit RGBA texture if RG32F is not available — pack U in R, V in G, quantise to 8-bit