Реакція-дифузія Грея-Скотта на GPU

1. Рівняння Грея-Скотта

Модель Грея-Скотта (1983, авторства Пітера Грея та Стівена Скотта) описує реакцію у безперервно перемішуваному резервуарі між двома хімічними речовинами U і V:

• U — це «їжа» (поживна речовина), що поповнюється ззовні зі швидкістю подачі F
• V — це «активатор», що утворюється реакцією U + 2V → 3V і знищується зі швидкістю (F+k)

Три члени в рівнянні U такі:

• D_u·∇²u — U дифундує по всій області (поширюється, згладжується)
• −u·v² — споживається автокаталітичною реакцією з V
• F·(1−u) — безперервно подається; діє як насичувальне джерело (швидше, коли u низьке)

2. Дискретний лапласіан на сітці

На 2D-сітці з кроком h = 1 піксель стандартний 5-точковий шаблонний лапласіан має вигляд:

9-точковий шаблон дає ізотропніші візерунки — плями й смуги без помітного зміщення вздовж осей сітки. Використовуйте його в шейдері.

3. Простір параметрів візерунків

Простір параметрів F-k — це справжній «зоопарк» якісно різних візерунків. Класичні іменовані області:

4. Архітектура GPU: ping-pong-фреймбуфери

Проблема: кожен крок симуляції читає поточний стан і записує наступний. На GPU не можна одночасно читати з тієї самої текстури й записувати в неї.

Рішення: ping-pong — виділіть два фреймбуфери (A і B). На кожному кроці:

• Читаємо з текстури A, записуємо у фреймбуфер B (шейдер симуляції)
• Наступний крок: читаємо з B, записуємо в A
• Міняємо посилання місцями та повторюємо

// Стан ping-pong
let frontBuffer = texA;  // GPU читає звідси
let backFBO     = fboB;  // GPU записує сюди

function swap() {
  [frontBuffer, backFBO] = [backBuffer.texture, frontFBO];
  // Після обміну: front ← те, що щойно записали; back ← старий буфер для перезапису на наступному кроці
}

Прохід рендерингу окремий: повноекранний квад, який читає з frontBuffer і відображає значення концентрацій у кольори для виводу. Він НЕ записує в текстури симуляції.

5. Налаштування WebGL 2

// Створюємо дві пари (текстура + фреймбуфер) для 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-канальний 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. Фрагментний шейдер симуляції

#version 300 es
precision highp float;

uniform sampler2D uState;   // RG = [u, v]
uniform vec2  uResolution;
uniform float uF;           // швидкість подачі
uniform float uK;           // швидкість знищення
uniform float uDu;          // дифузія U
uniform float uDv;          // дифузія V
uniform float uDt;          // крок за часом

out vec4 outColor;

// 9-точковий ізотропний лапласіан
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. Шейдер відображення кольорів

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

// Кольоровий градієнт у стилі Inferno для концентрації V
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. Повний драйвер на JavaScript

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;

// --- Компілюємо шейдери (допоміжні функції опущено задля стислості) ---
const simProgram    = compileProgram(gl, simVertSrc, simFragSrc);
const renderProgram = compileProgram(gl, renderVertSrc, renderFragSrc);

// --- VAO повноекранного квада ---
const quadVAO = makeQuadVAO(gl);

// --- Буфери ping-pong ---
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];
}

// --- Засіваємо початковий стан ---
const seed = new Float32Array(W * H * 2).fill(0);
for (let i = 0; i < W * H; i++) seed[i * 2] = 1.0; // u=1 усюди
for (let cy = 220; cy < 260; cy++) {          // засіяна ділянка
  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);

// --- Параметри симуляції ---
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();

Наступні кроки

• Додайте інтерактивний засів при кліку/перетягуванні — малюйте ділянки V=(0.25) там, де миша
• Змінюйте F і k у просторі — створіть карту градієнтів, щоб різні типи візерунків співіснували
• Спробуйте 3D → нарізайте вивід на різних площинах Z для об'ємних візерунків
• Перенесіть на 8-бітну текстуру RGBA, якщо RG32F недоступна — запакуйте U в R, V у G, квантуйте до 8 біт