Reakcja-dyfuzja: wzory Turinga na GPU

Równania Graya–Scotta

Model Graya–Scotta śledzi dwie substancje chemiczne, U i V, na dwuwymiarowej siatce. U jest ciągle dostarczane do systemu; V jest zużywane. Ich interakcja to reakcja U + 2V → 3V:

Dwa parametry, f i k, kontrolują wszystko. Drobne zmiany ich wartości dają radykalnie różne typy wzorów: plamki, paski, labirynty, dziury, koral, odciski palców.

Dlaczego CPU to złe miejsce na to

Dla siatki 512×512 to 262 144 komórek. Każda komórka potrzebuje wartości ze swoich czterech sąsiadów (Laplasjan), potem dwóch mnożeń zmiennoprzecinkowych i dwóch dodawań. Uruchomienie tego w JavaScript przy 60 FPS wymaga ~31 milionów operacji na sekundę — osiągalne, ale ledwo, i bez żadnego zapasu na większą siatkę czy jakiekolwiek przetwarzanie końcowe.

Rozwiązanie GPU jest oczywiste: aktualizacja Graya–Scotta to klasyczna operacja wybitnie równoległa — każda komórka jest niezależna. Shader fragmentów, który czyta z jednej tekstury i zapisuje do drugiej, wykonuje tę samą aktualizację dla wszystkich 262 144 komórek jednocześnie.

Bufory ping-pong

Nie można czytać i zapisywać tej samej tekstury w jednym przebiegu renderowania (niezdefiniowane zachowanie w WebGL). Rozwiązanie to ping-pong: utrzymuj dwa cele renderowania i zamieniaj co klatkę, z którego czytasz, a do którego zapisujesz.

// Konfiguracja: dwa identyczne cele renderowania
const size = 512;
const opts = {
  type: THREE.FloatType,        // float32 dla precyzji
  format: THREE.RGBAFormat,     // R = stężenie U, G = stężenie V
  minFilter: THREE.LinearFilter,
  magFilter: THREE.LinearFilter
};
let rtA = new THREE.WebGLRenderTarget(size, size, opts);
let rtB = new THREE.WebGLRenderTarget(size, size, opts);

// Co klatkę: zamiana
function tick() {
  // Krok 1: renderuj aktualizację Graya-Scotta (czytaj A, zapisz B)
  uniforms.tState.value = rtA.texture;
  renderer.setRenderTarget(rtB);
  renderer.render(quadScene, orthoCamera);

  // Krok 2: wyświetl wynik na ekranie
  displayUniforms.tState.value = rtB.texture;
  renderer.setRenderTarget(null);
  renderer.render(displayScene, orthoCamera);

  // Zamiana
  [rtA, rtB] = [rtB, rtA];
  requestAnimationFrame(tick);
}

Shader fragmentów aktualizacji

Shader aktualizacji GLSL próbkuje bieżącą teksturę stanu i stosuje równania Graya–Scotta. Dyskretny Laplasjan wykorzystuje szablon 5-punktowy — klasyczne przybliżenie różnicami skończonymi.

uniform sampler2D tState;   // bieżąca tekstura RG (R=U, G=V)
uniform vec2 uTexelSize;    // 1/szerokość, 1/wysokość
uniform float uF;           // szybkość podawania
uniform float uK;           // szybkość usuwania
uniform float uDu;          // dyfuzja U
uniform float uDv;          // dyfuzja V
uniform float uDt;          // krok czasowy na klatkę

void main() {
  vec2 uv = vUv;
  vec2 ts = uTexelSize;

  // Próbkuj bieżącą komórkę i 4 sąsiadów
  vec2 curr  = texture2D(tState, uv).rg;
  vec2 top   = texture2D(tState, uv + vec2( 0, ts.y)).rg;
  vec2 bot   = texture2D(tState, uv + vec2( 0,-ts.y)).rg;
  vec2 left  = texture2D(tState, uv + vec2(-ts.x, 0)).rg;
  vec2 right = texture2D(tState, uv + vec2( ts.x, 0)).rg;

  // Dyskretny Laplasjan 5-punktowy
  vec2 lap = top + bot + left + right - 4.0 * curr;

  float u = curr.r;
  float v = curr.g;
  float uvv = u * v * v;

  // Równania reakcji-dyfuzji Graya-Scotta
  float du = uDu * lap.r - uvv + uF * (1.0 - u);
  float dv = uDv * lap.g + uvv - (uF + uK) * v;

  gl_FragColor = vec4(
    clamp(u + du * uDt, 0.0, 1.0),
    clamp(v + dv * uDt, 0.0, 1.0),
    0.0, 1.0
  );
}

Zoo parametrów

Najbardziej fascynującym aspektem Graya–Scotta jest diagram fazowy. Trzymając D_u = 0,21 i D_v = 0,105 na stałe, oto wzory, które powstają przy zmianie samego F i k:

Zasiewanie: skąd biorą się wzory

Warunek początkowy ma ogromne znaczenie. Start z jednorodnego stanu nie daje żadnego wzoru — to stabilna równowaga. Aby wprawić system w tryb tworzenia wzorów, zasiewam małe „krople" wysokiego stężenia V w polu wysokiego U:

// Zasiej teksturę początkową (na CPU, wysyłana raz)
const data = new Float32Array(size * size * 4);
for (let i = 0; i < size * size; i++) {
  data[i * 4 + 0] = 1.0; // U = 1 wszędzie
  data[i * 4 + 1] = 0.0; // V = 0 wszędzie
}
// Wrzuć 20 losowych ziaren V
for (let s = 0; s < 20; s++) {
  const cx = Math.floor(Math.random() * size);
  const cy = Math.floor(Math.random() * size);
  for (let dx = -4; dx <= 4; dx++) {
    for (let dy = -4; dy <= 4; dy++) {
      const idx = ((cy + dy) * size + (cx + dx)) * 4;
      data[idx + 0] = 0.5; // zmniejsz U
      data[idx + 1] = 0.25; // wstrzyknij V
    }
  }
}
const initTexture = new THREE.DataTexture(data, size, size, THREE.RGBAFormat, THREE.FloatType);

Interaktywne malowanie

Najbardziej satysfakcjonująca funkcja: kliknij i przeciągnij, aby wstrzyknąć V bezpośrednio na teksturę w czasie działania. Uchwyt zdarzenia myszy konwertuje współrzędne ekranowe na przestrzeń UV, a następnie renderuje mały „pędzel" wysokiego V do bieżącego bufora ramki. Namalowane wzory rosną i ewoluują organicznie, łącząc się z istniejącą strukturą.

Wymagało to uwagi przy transformacjach współrzędnych — kwadrat jest rysowany w znormalizowanych współrzędnych urządzenia, ale pozycja myszy jest w pikselach CSS. Jedno wywołanie getBoundingClientRect() i konwersja z viewportu na UV rozwiązały to elegancko.

Co zaskoczyło mnie najbardziej

Dwie rzeczy. Po pierwsze, jak powolne są wzory w skali biologicznej — symulacja wykonuje tysiące kroków Eulera na sekundę, podczas gdy prawdziwe wzorowanie organizmu zachodzi w ciągu godzin do dni rozwoju embrionalnego. Oszałamiająco złożona struktura, którą widzisz krystalizującą się w sekundach w symulacji, wymaga od organizmu pełnego cyklu rozwojowego.

Po drugie, że te same równania opisują zjawiska w siedmiu rzędach wielkości skali: od milimetrowych łatek sygnalizacji komórkowej po pola wydm o rozmiarze kilometrów. Kombinacja Laplasjanu i nieliniowości reakcji wydaje się być uniwersalną maszyną formującą wzory, którą natura odkryła wcześnie i nigdy nie przestała używać.