Cahn-Hilliard Spinodal Decomposition

// Cahn-Hilliard spinodal decomposition on a periodic N x N grid.
//   du/dt = Laplacian( W'(u) - eps^2 * Laplacian(u) )
//   W(u)  = (1 - u^2)^2 / 4   =>   W'(u) = u^3 - u
// Explicit IMEX-ish solver: compute mu = W'(u) - eps^2 * lap(u),
// then advance u += dt * lap(mu).  The biharmonic term is the
// stability driver, so dt has to be small (we substep per frame).
// Mass is conserved exactly by construction (sum of Laplacian of
// anything periodic is zero up to roundoff).

const N = 128;
const EPS2 = 1.0;            // interface width^2 (in grid units)
// Explicit stability for the biharmonic requires dt < (dx)^4 / (16 * eps^2);
// with dx=1, eps^2=1 that's dt < 1/16 ≈ 0.0625. We pick 0.02 with a healthy
// margin (W''(u) at |u|=1 adds to the worst-case eigenvalue) and substep
// 40x per frame so coarsening still looks lively.
const DT_SIM = 0.02;         // sim timestep (was 0.05 — blew up to NaN)
const SUBSTEPS = 40;         // per render frame
const RESEED_SEC = 12;       // restart cycle for the loop
const U_CLAMP = 2.0;         // safety clamp on |u| after each step

let u, mu, lapU;
let img;                     // ImageData for the field buffer
let buf;                     // OffscreenCanvas at N x N
let bctx;
let W;
let H;
let elapsed;                 // sim-time since last reseed
let wallStart;               // wall clock at last reseed (for HUD)
let rngState;

// tiny deterministic PRNG (mulberry32) so reseeds are reproducible-ish
function rand() {
  rngState |= 0;
  rngState = (rngState + 0x6d2b79f5) | 0;
  let t = rngState;
  t = Math.imul(t ^ (t >>> 15), t | 1);
  t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
  return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
}

function reseed() {
  // Small zero-mean noise around u=0 — symmetric mixture.
  for (let i = 0; i < N * N; i++) {
    u[i] = (rand() - 0.5) * 0.2;
  }
  elapsed = 0;
}

// Periodic 5-point Laplacian into `out`.
function laplacian(field, out) {
  for (let j = 0; j < N; j++) {
    const jm = j === 0 ? N - 1 : j - 1;
    const jp = j === N - 1 ? 0 : j + 1;
    const row = j * N;
    const rowM = jm * N;
    const rowP = jp * N;
    for (let i = 0; i < N; i++) {
      const im = i === 0 ? N - 1 : i - 1;
      const ip = i === N - 1 ? 0 : i + 1;
      out[row + i] =
        field[row + ip] + field[row + im] +
        field[rowP + i] + field[rowM + i] -
        4 * field[row + i];
    }
  }
}

function step(dt) {
  // mu = u^3 - u - eps^2 * lap(u)
  laplacian(u, lapU);
  for (let k = 0; k < N * N; k++) {
    const v = u[k];
    mu[k] = v * v * v - v - EPS2 * lapU[k];
  }
  // u += dt * lap(mu), then clamp to keep the explicit scheme from
  // running away if the user pokes a bad value or numeric drift hits.
  laplacian(mu, lapU);
  for (let k = 0; k < N * N; k++) {
    let v = u[k] + dt * lapU[k];
    if (!(v > -U_CLAMP)) v = -U_CLAMP;       // catches NaN too
    else if (v > U_CLAMP) v = U_CLAMP;
    u[k] = v;
  }
}

// Divergent purple <-> off-white <-> orange colormap.
// u in roughly [-1, 1].  Map to t in [0, 1].
const C_PURPLE = [60, 30, 90];     // deep purple
const C_MID    = [245, 240, 230];  // warm off-white
const C_ORANGE = [225, 120, 40];   // burnt orange

function colormap(v, out, off) {
  // clamp to [-1.15, 1.15] and rescale
  let t = v;
  if (t < -1.15) t = -1.15;
  else if (t > 1.15) t = 1.15;
  t = (t + 1.15) / 2.3;
  let r, g, b;
  if (t < 0.5) {
    const s = t / 0.5;
    r = C_PURPLE[0] + (C_MID[0] - C_PURPLE[0]) * s;
    g = C_PURPLE[1] + (C_MID[1] - C_PURPLE[1]) * s;
    b = C_PURPLE[2] + (C_MID[2] - C_PURPLE[2]) * s;
  } else {
    const s = (t - 0.5) / 0.5;
    r = C_MID[0] + (C_ORANGE[0] - C_MID[0]) * s;
    g = C_MID[1] + (C_ORANGE[1] - C_MID[1]) * s;
    b = C_MID[2] + (C_ORANGE[2] - C_MID[2]) * s;
  }
  out[off] = r;
  out[off + 1] = g;
  out[off + 2] = b;
  out[off + 3] = 255;
}

function render() {
  const d = img.data;
  for (let k = 0; k < N * N; k++) {
    colormap(u[k], d, k * 4);
  }
  bctx.putImageData(img, 0, 0);
}

function init({ ctx, width, height }) {
  W = width;
  H = height;
  u = new Float32Array(N * N);
  mu = new Float32Array(N * N);
  lapU = new Float32Array(N * N);
  buf = new OffscreenCanvas(N, N);
  bctx = buf.getContext('2d');
  img = bctx.createImageData(N, N);
  rngState = 0x9e3779b9;
  reseed();
  // burn in a few steps so frame 1 isn't flat noise
  for (let i = 0; i < 40; i++) step(DT_SIM);
  wallStart = 0;
  // background
  ctx.fillStyle = '#1a0e22';
  ctx.fillRect(0, 0, W, H);
}

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

  // tap to re-noise (also resets the 12s loop timer)
  if (input && input.consumeClicks && input.consumeClicks() > 0) {
    reseed();
    wallStart = time;
    for (let i = 0; i < 40; i++) step(DT_SIM);
  }

  // advance sim
  for (let i = 0; i < SUBSTEPS; i++) step(DT_SIM);
  elapsed += dt;

  // reseed loop
  if (elapsed > RESEED_SEC) {
    reseed();
    wallStart = time;
    for (let i = 0; i < 40; i++) step(DT_SIM);
  }

  render();

  // letterbox background
  ctx.fillStyle = '#1a0e22';
  ctx.fillRect(0, 0, W, H);

  // fit the square field into the canvas, centered
  const side = Math.min(W, H);
  const ox = (W - side) * 0.5;
  const oy = (H - side) * 0.5;
  ctx.imageSmoothingEnabled = true;
  ctx.drawImage(buf, ox, oy, side, side);

  // subtle elapsed-time tick in the corner (Peyré-style minimal HUD)
  const t = elapsed.toFixed(1);
  ctx.font = '11px ui-monospace, monospace';
  ctx.fillStyle = 'rgba(245, 240, 230, 0.55)';
  ctx.textBaseline = 'top';
  ctx.fillText(`t = ${t}`, ox + 8, oy + 6);
}