Gray-Scott Reaction-Diffusion

Math essence

dU/dt = Du*lap(U) - U*V*V + F*(1 - U); dV/dt = Dv*lap(V) + U*V*V - (F + k)*V. The isotropic 5-point + 4-diagonal 9-point Laplacian is computed with numpy.roll for periodic boundaries; integration is explicit Euler with clipping to [0,1].

Parameters

Reference

First 60 lines of reference.py. Full source: library_external/reas_processing/grayScottRD/reference.py

from __future__ import annotations

from typing import ClassVar

import numpy as np


class GrayScottRD:
    """Gray-Scott Reaction-Diffusion.

    Two scalar species U and V diffuse on a square grid with rates Du and Dv
    and react with feed rate F and kill rate k:

        dU/dt = Du*lap(U) - U*V*V + F*(1 - U)
        dV/dt = Dv*lap(V) + U*V*V - (F + k)*V

    The system is the canonical Gray-Scott model from Pearson 1993, popular in
    generative art since the early 2000s (it appears in the Processing
    community's Reaction-Diffusion tutorials). Different (F, k) pairs produce
    mitosis spots, stripes, coral, worms, solitons, and other Turing patterns.
    """

    family: ClassVar[str] = "reas_processing"
    sketch: ClassVar[str] = "grayScottRD"

    def __init__(
        self,
        dim: int = 256,
        Du: float = 1.0,
        Dv: float = 0.5,
        F: float = 0.0367,
        k: float = 0.0649,
        dt: float = 1.0,
        steps_per_render: int = 16,
        palette_mode: str = "viridis",
        seed: int | None = None,
        initial_seeds: int = 8,
    ) -> None:
        self.dim = dim
        self.Du = Du
        self.Dv = Dv
        self.F = F
        self.k = k
        self.dt = dt
        self.steps_per_render = steps_per_render
        self._palette_mode = palette_mode
        self.initial_seeds = initial_seeds
        self._seed = seed
        self.reset()

    def reset(self, seed: int | None = None) -> None:
        if seed is not None:
            self._seed = seed
        self.rng = np.random.default_rng(self._seed)
        self.U = np.ones((self.dim, self.dim), dtype=np.float32)
        self.V = np.zeros((self.dim, self.dim), dtype=np.float32)
        n = self.initial_seeds
        cx = self.rng.integers(self.dim // 4, 3 * self.dim // 4, size=n)
        cy = self.rng.integers(self.dim // 4, 3 * self.dim // 4, size=n)
        rad = np.arange(-3, 4)