# Colour the Lattice

*Task 4 of 6 · [The Ising Model: Colour to Break the Race](https://gpu.rocks/learn/ising-model-1f12d841.md) · GPU.js Learn*

The cure is a chessboard, and if you have been through **Iterative Linear
Solvers** you have already met it: red-black Gauss-Seidel colours a grid in exactly this
way, for exactly this reason. One idea in two costumes. There it rescues a solver that is
sequential by construction; here it repairs a Monte Carlo update that is silently wrong when
it is parallel. Both times the argument is one sentence — *the stencil reads only the four
direct neighbours, and on a chessboard every direct neighbour of a red square is black*.

So call a cell **red** when `(x + y)` is even and
**black** when it is odd, and update all 8,192 reds at once. No red cell reads
another red cell, so nothing a red thread looked at can move while it is deciding: the
`ΔE` it computed is exact, not an estimate, and the flip it accepts is the flip
Metropolis meant. Then do the blacks, reading the reds that were just written. Two data-parallel
half-sweeps, and the race is gone — not mitigated, gone.

Black cells are not "skipped" in the red half. Every thread still writes its own cell;
a black thread writes back the value it already had, ready for the half-sweep that is about
to need it exactly as it is.

## Goal

**Goal:** write the red half-sweep — cells with
`(x + y) % 2 === 0` take the Metropolis decision, every other cell comes through
untouched.

## Requirements

- Keep the parity in a *number*: `const parity = (x + y) % 2;` — a boolean in a kernel variable does not compile on WebGL
- Cells with parity `1` (black) return their spin unchanged
- Red cells with `dE ≤ 0` flip unconditionally
- Red cells with `dE > 0` flip when `u < Math.exp(-dE / temperature)`

## Hint 1 — the trap this task is built around

The natural spelling is a boolean, and it is the one thing gpu.js cannot do:

```js
const isRed = (x + y) % 2 === 0;   // throws on WebGL
```

The GL backend has no way to store a `bool` in a kernel variable, so it fails at
shader-compile time with *cannot convert from 'bool' to 'lowp float'* — while the CPU
backend runs it happily, which is how this reaches production. Keep the number:

```js
const parity = (x + y) % 2;        // 0 or 1
if (parity !== 0) return spin;
```

## Hint 2 — the Metropolis decision

Two exits, in this order:

```js
if (dE <= 0) return -spin;
if (u < Math.exp(-dE / temperature)) return -spin;
return spin;
```

The first line is technically implied by the second — `u` is always below 1 and
`exp(−dE/T)` is at least 1 whenever `dE ≤ 0`, so the test would pass
anyway. Write it out regardless: it is the rule, and it says that a flip which lowers the
energy is never a gamble.

## Hint 3 — the shape of the body

The parity guard comes first, so a black thread never computes a neighbour sum it is
not going to use:

```js
const spin = s[y][x];
const parity = (x + y) % 2;
if (parity !== 0) return spin;
const sum = s[y][(x + 1) % n] + s[y][(x + n - 1) % n]
          + s[(y + 1) % n][x] + s[(y + n - 1) % n][x];
const dE = 2 * spin * sum;
```

Then the hash (already written for you) and the two exits above.

## Same idea elsewhere

Red-black ordering, and its multi-colour generalisation, is the standard way to put
a Gauss-Seidel smoother, a lattice Monte Carlo or a physics solver's constraint pass on a GPU:
one dispatch per colour, and an unstructured mesh gets its colours from a graph-colouring pass
first. The idea is bigger than any of them — a colour is simply a set of updates guaranteed
not to depend on each other, which is the same permission slip a wavefront's anti-diagonal or
a task graph's level hands out.

## Starter code

```js
// Half the lattice at a time. Reds first — and no red reads a red.
const gpu = new GPU({ mode });

const red = gpu.createKernel(function (s, temperature, seed) {
  const x = this.thread.x;
  const y = this.thread.y;
  const n = this.constants.size;
  const spin = s[y][x];
  // TODO 1: the parity, as a NUMBER — const parity = (x + y) % 2;
  //         (a boolean in a kernel variable will not compile on WebGL)
  // TODO 2: parity 1 is black — return spin unchanged.
  // TODO 3: sum the four wrapped neighbours and form dE = 2 * spin * sum.

  // a random number for THIS thread, this sweep — no shared state, no ordering
  let h = (x * 1103 + y * 2749 + seed * 3571) % 65536;
  let q = h % 2048;
  h = (h + q * q) % 65536;
  h = (h % 256) * 256 + Math.floor(h / 256);
  h = (h * 253 + 30011) % 65536;
  q = h % 2048;
  h = (h + q * q) % 65536;
  const u = h / 65536;

  // TODO 4: dE <= 0 flips unconditionally; otherwise flip when u < exp(-dE / temperature).
  return spin;
}, { output: [128, 128], constants: { size: 128 } });

const after = await red(lattice, 1.5, 0);

let movedRed = 0;
let movedBlack = 0;
for (let y = 0; y < 128; y++) {
  for (let x = 0; x < 128; x++) {
    if (after[y][x] === lattice[y][x]) continue;
    if ((x + y) % 2 === 0) movedRed++;
    else movedBlack++;
  }
}
console.log('red cells that flipped:', movedRed, 'of 8192');
console.log('black cells that flipped:', movedBlack, '(must be 0)');
```

---

Interactive version: https://gpu.rocks/learn/ising-model-1f12d841/4

[Previous task](https://gpu.rocks/learn/ising-model-1f12d841/3.md) · [Next task](https://gpu.rocks/learn/ising-model-1f12d841/5.md)