# Softening the Singularity

*Task 3 of 5 · [N-Body Gravity](https://gpu.rocks/learn/n-body-gravity-5de47751.md) · GPU.js Learn*

Two of this task's bodies sit `0.001` apart. Plug that into
`1 / r²` and their mutual pull is about a *million* — one tick of the
clock later they're flung out of the galaxy. That's not physics; it's what happens when a
point-mass model meets a finite time step.

The standard fix is **Plummer softening**: replace `r²` with
`r² + ε²`. Far away, `ε` changes nothing; up close, the force
flattens out instead of diverging. Bonus: the `j !== i` self-check becomes dead
weight — your own term has `dx = dy = 0`, so it contributes exactly zero. Drop
the branch; GPUs run happiest when every thread takes the same path.

## Figures

- **close encounters flatten out instead of blowing up**

## Goal

**Goal:** soften the kernel — use `r² + soft²`, drop the
self-check, and return the full `[ax, ay]` pair.

## Requirements

- Squared distance becomes `dx·dx + dy·dy + soft·soft`
- Remove the `j !== this.thread.x` guard — the self term is now zero
- Accumulate *both* components and return `[ax, ay]`

## Hint 1 — why the guard can go

For `j === i`: `dx` and `dy` are 0, so the
contribution is `0 · something`. With `soft² > 0` the
denominator is never zero, so that something is a plain finite number.

## Hint 2 — share the weight

Compute `const w = mass[j] / (r2 * Math.sqrt(r2));` once, then
`ax += dx * w; ay += dy * w;` — one denominator, two components.

## Same idea elsewhere

Softening appears verbatim in production astrophysics codes (GADGET, Bonsai) on
CUDA and ROCm clusters. It's also a lesson in GPU numerics generally: shader float math
never throws — a divide-by-zero silently mints `Infinity` and then
`NaN`s spread through every sum they touch, on Metal and WebGPU alike.

## Starter code

```js
// Bodies 0 and 1 sit 0.001 apart. Unsoftened, their mutual pull
// is ~a million — one bad pair and the whole simulation explodes.
const gpu = new GPU({ mode });

const accel = gpu.createKernel(function (posX, posY, mass, soft) {
  const myX = posX[this.thread.x];
  const myY = posY[this.thread.x];
  let ax = 0;
  let ay = 0;
  for (let j = 0; j < this.constants.n; j++) {
    if (j !== this.thread.x) {
      const dx = posX[j] - myX;
      const dy = posY[j] - myY;
      // TODO: soften — add soft·soft to r² so close encounters stay
      // finite. Then the j !== i guard above is dead weight: delete it.
      const r2 = dx * dx + dy * dy;
      const w = mass[j] / (r2 * Math.sqrt(r2));
      ax += dx * w;
      ay += dy * w;
    }
  }
  return [ax, ay];
}, { output: [64], constants: { n: 64 } });

const acc = await accel(posX, posY, mass, 0.1);
console.log('acceleration on body 0:', acc[0][0], acc[0][1]);
```

---

Interactive version: https://gpu.rocks/learn/n-body-gravity-5de47751/3

[Previous task](https://gpu.rocks/learn/n-body-gravity-5de47751/2.md) · [Next task](https://gpu.rocks/learn/n-body-gravity-5de47751/4.md)