`IBMStepper`: the IBM loop is not a fixed point iteration — `ibm_tolerance` cannot fire, and `ibm_relaxation`/`ibm_max_iterations` act as a forcing gain

## Summary

In `xlb/operator/stepper/ibm_stepper.py`, the loop in `warp_implementation` never updates `f_1`. The only kernel that writes `f_1` — `correct_population_ibm` — is launched once *after* the loop exits. As a result every sweep recomputes the same `u` from the same populations, so the loop does not converge in the sense the docstring describes.

This makes the documented behaviour

> 3) Residual-based stopping instead of a fixed iteration count. The loop monitors the maximum incremental change in Lagrangian forces and stops early when it falls below `ibm_tolerance` [...]

unreachable in practice.

## Mechanism

Inside the loop:

- `compute_velocity_and_correct` reads `f_1` to get `u` and writes `eul_velocities`, but does not modify `f_1`.
- `interpolate_velocity_and_update_force` interpolates `eul_velocities` back to the markers and does `lag_forces += v_solid - u_interp`.

Since `f_1` is unchanged, `u` is identical on every sweep, hence `delta_F = v_solid - u_interp` is identical on every sweep. Therefore:

- `lag_forces` after sweep `k` is exactly `k * delta_F` (a linear ramp)
- the residual `|lag_forces - lag_forces_prev|` equals `|delta_F|` — **a constant**

A constant residual can never cross a threshold it did not already satisfy on the first sweep.

## Reproduction

A faithful 1-D transcription of the kernel sequence (Peskin weights, Voronoi areas, same spread/interpolate operators), `omega=1.0`, 6 sweeps:

```
residual : 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
|lag|max : 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000
```

The residual is bit-for-bit constant; the accumulated force grows linearly.

We also instrumented the real Warp path on GPU and counted iterations: **100% of steps hit the `ibm_max_iterations` ceiling with the convergence flag still set**, across two different resolutions and two different `ibm_relaxation` values (60 steps each). That matches the analysis.

## Consequences

**1. `ibm_tolerance` is effectively dead.** The early exit never triggers; the loop always runs `ibm_max_iterations` times.

**2. `ibm_relaxation` and `ibm_max_iterations` are a boundary forcing gain, not convergence controls.** Working through the algebra, the correction applied after the loop is

```
eul_forces = omega * ( (n-1) * spread(delta_F)/weights - u )
```

so the gain on the target term is `ibm_relaxation * (ibm_max_iterations - 1)`. Note the `- 1`: the final sweep spreads the forces accumulated by the *previous* sweep.

**3. `ibm_max_iterations = 1` silently drops the wall velocity entirely.** With `n = 1` the target gain is zero and only `-omega * u` is applied. For a static body that still drives `u -> 0`, which is the correct no-slip target, so it looks fine. For a **moving** body the prescribed surface velocity `v_solid` never enters the equations at all. This one is easy to miss because the static case masks it.

**4. Stability.** Because the scheme is an explicit forcing with gain rather than a contraction, the near-surface amplification factor flips sign and exceeds 1 as the gain rises. In our reference implementation, for a static body: `omega=1.0, n=4` gives **-2.61**, `omega=0.5` gives `-0.805`, `omega=0.25` gives `+0.098`. We see the corresponding ordering in practice — at fine grid resolutions the run diverges at the higher gains and survives longer as the gain is lowered.

(The absolute numbers overstate the real amplification, since the reference implementation omits the collide/stream between corrections. The robust findings are the ordering and the sign flip.)

## Suggested fix

Apply each sweep's correction to `f_1` inside the loop, so the next sweep reads an updated `u`. Two details matter:

- spread the **per-sweep increment**, not the accumulated `lag_forces` — the accumulated total is already baked into `f_1` by the previous sweeps' corrections, so spreading it again double-counts;
- drop the `- u` term in the normalization for the same reason.

With that change the residual actually decreases, `ibm_tolerance` and the early exit become meaningful, and the surface correction is monotonically damping for `0 < ibm_relaxation <= 1` (no sign flip at any of the gains above).

We have implemented this in our fork behind an opt-in flag, keeping the current behaviour as the default so existing results stay reproducible. Happy to open a PR against `main` if that shape is useful to you — either the flag approach or a straight fix.

## Environment

- `xlb` `main` (the code path above is unchanged as of today)
- Warp backend, D3Q27, `warp-lang` 1.12.1

## Alternative reading

If the current behaviour is intentional — i.e. this is meant to be a single explicit forcing step with a tunable gain rather than a multi-direct forcing iteration — then the fix is documentation rather than code: the `ibm_tolerance` parameter should be removed or marked non-functional, and the docstring's points (2) and (3) reworded, since `ibm_relaxation` and `ibm_max_iterations` are gain knobs. Either way the current docs describe behaviour the code does not have.

Start in xlb/operator/stepper/ibm_stepper.py at warp_implementation and trace compute_velocity_and_correct, interpolate_velocity_and_update_force, and correct_population_ibm to confirm how f_1 and lag_forces change across sweeps. Reproduce the reported residual and iteration behavior, then determine whether the intended result is an iterative correction or documented explicit forcing; done means the chosen behavior and its tolerance, relaxation, and max-iteration semantics are covered by validation.

Written by the indexing model from the issue text.