Allen-Cahn Phase Separation

// Allen-Cahn equation: phase separation under the L^2 gradient flow of the
// Ginzburg-Landau energy. Solved on a periodic grid by semi-implicit Euler.
//
// Inspired by Gabriel Peyre's gradient-flow visualizations.

const GW = 200;
const GH = 150;
const N = GW * GH;

// Numerical parameters.
const EPS = 0.05;          // interface width (in units of grid spacing L=1)
const EPS2 = EPS * EPS;
const DT = 0.08;           // time step (semi-implicit so generous)
const SUBSTEPS = 2;        // sim steps per render frame
const RESEED_SECONDS = 10; // re-noise + restart every ~10 s

let u;        // current field
let uNext;    // scratch
let buf;      // pixel buffer (GW x GH)
let bufCtx;
let bufImg;

let W, H;
let simTime;
let secondsSinceReseed;

let prevMouseDown;
let mouseStrokeSign;

function reseed() {
  for (let i = 0; i < N; i++) {
    u[i] = (Math.random() * 2 - 1) * 0.4;
  }
  simTime = 0;
  secondsSinceReseed = 0;
}

// W'(u) = u^3 - u
function step(dtSim) {
  // semi-implicit Euler:
  //   (u^{n+1} - u^n)/dt = eps^2 * Lap(u^{n+1}) - ((u^n)^3 - u^n)
  // Approximate with a single Jacobi sweep on the implicit Laplacian:
  //   u^{n+1} ~ ( u^n + dt*u^n - dt*(u^n)^3 + dt*eps^2 * lap5(u^n) ) / 1
  // This is effectively explicit but the small eps^2 and modest dt keep it stable.
  // (Stability: dt < 1 for the reaction part; dt*eps^2*4 < 1 for diffusion.
  //  With eps=0.05 and dt=0.08, dt*eps^2*4 = 8e-4 << 1.)
  for (let y = 0; y < GH; y++) {
    const ym = (y === 0 ? GH - 1 : y - 1) * GW;
    const yp = (y === GH - 1 ? 0 : y + 1) * GW;
    const y0 = y * GW;
    for (let x = 0; x < GW; x++) {
      const xm = x === 0 ? GW - 1 : x - 1;
      const xp = x === GW - 1 ? 0 : x + 1;
      const c = u[y0 + x];
      const lap = u[y0 + xm] + u[y0 + xp] + u[ym + x] + u[yp + x] - 4 * c;
      // reaction: -W'(u) = u - u^3
      const react = c - c * c * c;
      uNext[y0 + x] = c + dtSim * (EPS2 * lap + react);
    }
  }
  const tmp = u;
  u = uNext;
  uNext = tmp;
}

// Diverging red/blue colormap. v in [-1,1] -> rgb.
// v = -1: cool blue. v = +1: warm red. v near 0: near-white.
function colormap(v, out, off) {
  // clamp
  if (v > 1) v = 1;
  else if (v < -1) v = -1;
  // map magnitude to a saturation curve so the interface is a thin white band
  const a = Math.abs(v);
  // ease-out so values near 0 stay near white
  const s = a; // linear is fine; thin interface comes from the dynamics
  let r, g, b;
  if (v >= 0) {
    // white -> red
    r = 255;
    g = Math.round(255 * (1 - 0.78 * s));
    b = Math.round(255 * (1 - 0.85 * s));
  } else {
    // white -> blue
    r = Math.round(255 * (1 - 0.85 * s));
    g = Math.round(255 * (1 - 0.55 * s));
    b = 255;
  }
  out[off] = r;
  out[off + 1] = g;
  out[off + 2] = b;
  out[off + 3] = 255;
}

function paintField() {
  const data = bufImg.data;
  let off = 0;
  for (let i = 0; i < N; i++) {
    colormap(u[i], data, off);
    off += 4;
  }
  bufCtx.putImageData(bufImg, 0, 0);
}

function paintBrush(canvasX, canvasY, sign) {
  // map canvas pixel to grid cell
  const gx = Math.floor((canvasX / W) * GW);
  const gy = Math.floor((canvasY / H) * GH);
  if (gx < 0 || gx >= GW || gy < 0 || gy >= GH) return;
  const R = 6; // grid cells
  const R2 = R * R;
  for (let dy = -R; dy <= R; dy++) {
    const yy = gy + dy;
    if (yy < 0 || yy >= GH) continue;
    for (let dx = -R; dx <= R; dx++) {
      const xx = gx + dx;
      if (xx < 0 || xx >= GW) continue;
      const d2 = dx * dx + dy * dy;
      if (d2 > R2) continue;
      const w = 1 - d2 / R2;
      const idx = yy * GW + xx;
      // push toward sign
      u[idx] = u[idx] * (1 - w) + sign * w;
    }
  }
}

function init({ ctx, width, height }) {
  W = width;
  H = height;

  u = new Float32Array(N);
  uNext = new Float32Array(N);

  // Offscreen pixel buffer at simulation resolution.
  buf = new OffscreenCanvas(GW, GH);
  bufCtx = buf.getContext('2d');
  bufImg = bufCtx.createImageData(GW, GH);

  reseed();

  prevMouseDown = false;
  mouseStrokeSign = 1;

  ctx.fillStyle = '#ffffff';
  ctx.fillRect(0, 0, W, H);
}

function tick({ ctx, dt, width, height, input }) {
  if (width !== W || height !== H) {
    W = width;
    H = height;
  }

  // Painting: on mouse press, pick a sign opposite to the local field so the
  // stroke actually creates a perturbation. While held, keep the same sign.
  if (input.mouseDown) {
    if (!prevMouseDown) {
      const gx = Math.floor((input.mouseX / W) * GW);
      const gy = Math.floor((input.mouseY / H) * GH);
      if (gx >= 0 && gx < GW && gy >= 0 && gy < GH) {
        mouseStrokeSign = u[gy * GW + gx] >= 0 ? -1 : 1;
      } else {
        mouseStrokeSign = 1;
      }
    }
    paintBrush(input.mouseX, input.mouseY, mouseStrokeSign);
  }
  prevMouseDown = input.mouseDown;

  // Advance the PDE.
  for (let s = 0; s < SUBSTEPS; s++) {
    step(DT);
    simTime += DT;
  }

  secondsSinceReseed += dt;
  if (secondsSinceReseed >= RESEED_SECONDS) {
    reseed();
  }

  // Render: paint the grid to an offscreen buffer, then upscale to the canvas
  // with nearest-neighbor disabled for a soft look; the interface stays thin
  // because the dynamics make it thin, not the rendering.
  paintField();

  ctx.imageSmoothingEnabled = true;
  ctx.imageSmoothingQuality = 'high';
  ctx.drawImage(buf, 0, 0, W, H);

  // Tiny time readout in the corner.
  ctx.font = '12px ui-sans-serif, system-ui, sans-serif';
  ctx.textBaseline = 'top';
  ctx.fillStyle = 'rgba(30, 30, 30, 0.55)';
  ctx.fillText('t = ' + simTime.toFixed(1), 10, 8);
}