# Cross the Line

*Task 2 of 5 · [The Heat Equation & Stability](https://gpu.rocks/learn/heat-and-stability-514063bb.md) · GPU.js Learn*

Nothing about the kernel you just wrote is wrong. Hand it a step size that is
slightly too big and it will still compute exactly what you asked for — and the field
will be at 10²⁴ in eighty steps.

Look again at that centre weight, `1 − 4α`. At `α = 0.1` it is
`+0.6`, the new value really is an average of five old ones, and an average
can never leave the range of the things it averaged. At `α = 0.4` it is
`−0.6`, and "average" has become a lie: a cell that sits above its
neighbours is now pushed *further* above them every step. In
`dims` dimensions the centre weight is `1 − 2·dims·α`, so the
tipping point is `α = 1/(2·dims)` — which, written in the quantities you
actually control, is the limit every explicit solver lives under:

```js
dt  ≤  dx² / (2 · D · dims)
       dims = 2 on a 2D grid
```

Past it, the fastest pattern the grid can hold — a checkerboard of alternating hot
and cold cells — is multiplied by `|1 − 8α|` every step instead of damped.
At `α = 0.4` that is **2.2× per step**, and there is always
some checkerboard in there, if only from rounding. Ten steps: ×2,700. Eighty steps:
the run below.

## Figures

- **the fastest pattern the grid can hold, flipping and growing, every step**

## Goal

**Goal:** work out `dtMax`, then run the same eighty-step
simulation twice — once safely inside the limit, once past it — and watch the console.
The stepping loop is written for you: N-Body already made a meal of *how* to
advance a simulation, and the only thing that matters here is how far each step goes.

## Requirements

- Compute `dtMax = dx * dx / (2 * D * DIMS)` — it is `0.5` for this grid
- Call `run` once at `0.4 * dtMax` and once at `1.6 * dtMax`
- Leave `STEPS` at 80 — the trace prints the hottest cell every 20 steps

## Hint 1 — where the numbers come from

Each of the `2 · dims` neighbours takes an `α`-sized bite
out of the centre, so the centre keeps `1 − 2·dims·α`. Set that to zero,
substitute `α = D·dt/dx²`, and solve for `dt`.

## Hint 2 — the two runs

```js
const dtMax = dx * dx / (2 * D * DIMS);
await run('SAFE', 0.4 * dtMax);
await run('PAST THE LINE', 1.6 * dtMax);
```

## Same idea elsewhere

Every production explicit solver computes this number and refuses to exceed it:
CFL conditions in fluid codes, the diffusion-number check in a thermal simulation, the
substepping loop in a cloth or fluid solver on the GPU. It is also why GPU simulations
so often become *launch-bound* — halving `dx` to sharpen a picture
quarters the legal `dt`, so the same second of simulated time costs four
times as many kernel launches.

## Starter code

```js
// Two runs of the same simulation. Only the step size differs.
const gpu = new GPU({ mode });

const D = 8;      // diffusivity
const dx = 4;     // cell spacing
const DIMS = 2;   // a 2D grid: four neighbours
const STEPS = 80;

// TODO: the explicit stability limit — dt ≤ dx² / (2 · D · dims)
const dtMax = 0;

// max |u| over the field, written so a NaN cannot hide: every comparison with
// NaN is false, so `if (a > m)` would skip it. `!(a <= m)` is true for NaN.
function hottest(u) {
  let m = 0;
  for (let y = 0; y < u.length; y++) {
    for (let x = 0; x < u[y].length; x++) {
      const a = Math.abs(u[y][x]);
      if (!(a <= m)) m = a;
    }
  }
  return m;
}

async function run(label, dt) {
  const alpha = D * dt / (dx * dx);
  const step = gpu.createKernel(function (u) {
    const x = this.thread.x;
    const y = this.thread.y;
    let xl = x - 1; if (xl < 0) xl = this.constants.size - 1;
    let xr = x + 1; if (xr > this.constants.size - 1) xr = 0;
    let yd = y - 1; if (yd < 0) yd = this.constants.size - 1;
    let yu = y + 1; if (yu > this.constants.size - 1) yu = 0;
    const c = u[y][x];
    const lap = u[y][xl] + u[y][xr] + u[yd][x] + u[yu][x] - 4 * c;
    return c + this.constants.alpha * lap;
  }, { output: [48, 48], constants: { size: 48, alpha } });

  console.log(label, '— dt =', dt, ' alpha =', alpha, ' centre weight =', 1 - 4 * alpha);
  let u = seed;
  for (let i = 1; i <= STEPS; i++) {
    u = await step(u);
    if (i % 20 === 0) console.log('   step', i, '→ hottest |u| =', hottest(u));
  }
  return u;
}

console.log('stability limit: dt <=', dtMax);

// TODO: run twice — at 0.4 * dtMax, then at 1.6 * dtMax.
//       run() is async now (it awaits the kernel), so await each call.
```

---

Interactive version: https://gpu.rocks/learn/heat-and-stability-514063bb/2

[Previous task](https://gpu.rocks/learn/heat-and-stability-514063bb/1.md) · [Next task](https://gpu.rocks/learn/heat-and-stability-514063bb/3.md)