docs: correct why IndexSwarmVariable level sets keep inverse distance

docs/developer/CHANGELOG.md

@@ -40,10 +40,11 @@ first order with refinement and **not at all** with stencil size.
 Two deliberate exclusions, both measured rather than assumed:
 `MeshVariable.rbf_interpolate` keeps inverse distance because it is the
 fallback rung of the point-location ladder, whose documented contract is that
-it is bounded; and `IndexSwarmVariable` material level sets keep it because a
-material indicator is piecewise constant — there is nothing for linear
-exactness to gain at a discontinuity, and signed weights push level sets
-outside `[0, 1]`.
+it is bounded; and `IndexSwarmVariable` material level sets keep it because
+they estimate a fraction from a handful of *integer* samples, where the error
+is dominated by variance rather than bias. Signed weights amplify that variance
+by roughly an order of magnitude and push level sets outside `[0, 1]`, while
+the bias they would remove is already negligible.
 
 Related: swarm proxy refresh no longer fails under an active units model
 (#426, #434); the units the proxy advertises are tracked separately (#439).

docs/developer/subsystems/interpolation.md

@@ -234,48 +234,71 @@ Note what this does and does not promise: it bounds *new oscillation relative
 to the local trend*, not absolute range. A quantity that must stay inside hard
 physical bounds (a fraction in $[0,1]$) needs its own clip on top.
 
-### Material level sets stay on `order=0` — measured, not assumed
+### Material level sets stay on `order=0` — variance, not bias
 
 `IndexSwarmVariable` builds one level-set MeshVariable per material index and
 keeps its own inverse-distance weighting. It was tested against `order=1` and
-**deliberately not changed**.
+**deliberately not changed** — but not for the reason one might expect.
 
-The reason is structural. A material indicator is **piecewise constant**, not
-smooth. Away from an interface both schemes reproduce it exactly, because both
-reproduce constants — linear exactness has nothing to add. At the interface the
-field is *discontinuous*, so no polynomial-reproducing scheme is exact either;
-signed weights simply add overshoot where the data has a jump.
+The tempting argument, *a material indicator is piecewise constant so linear
+exactness has nothing to gain*, is **wrong**. The indicator is the *particle*
+data; the quantity estimated at a node is the local material **fraction**,
+which near an interface is a smooth ramp. Reproducing a ramp is exactly what a
+constants-only scheme cannot do.
 
-Measured on a straight interface at `x = 0.5` (exactly representable, so any
-displacement of the recovered 0.5 contour is scheme error):
+The real reason is that a level-set node estimates that fraction from a handful
+of **integer** samples, so its error has a variance term as well as a bias term:
 
-| scheme | interface error, median | level-set range |
-|---|---|---|
-| inverse distance, `nnn=5` | 5.2e-3 – 1.1e-2 | `[0, 1]` exactly |
-| `order=1`, `nnn=6` | 6.9e-3 – 7.8e-3 | `[0, 1]` exactly |
-| `order=1`, `nnn=8` | 5.2e-3 – 7.9e-3 | **`[-0.038, 1.038]`** |
+$$ \mathrm{Var}(\text{estimate}) = \sum_j w_j^2 \,\mathrm{Var}(f_j) $$
+
+For weights summing to one, uniform weighting minimises $\sum_j w_j^2$ at
+$1/nnn$. Inverse-distance weights are positive and stay near that floor;
+linear-exact weights are signed and sit far above it:
+
+| dim | `nnn` | $\sum w^2$ inverse distance | $\sum w^2$ `order=1` | $1/nnn$ | amplification |
+|---|---|---|---|---|---|
+| 2 | 6 | 0.223 | 2.314 | 0.167 | **10.4x** |
+| 2 | 20 | 0.079 | 0.922 | 0.050 | **11.7x** |
+| 3 | 6 | 0.188 | 1.290 | 0.167 | **6.9x** |
+| 3 | 20 | 0.061 | 0.657 | 0.050 | **10.9x** |
+
+Note the second row of each pair: inverse distance averages the noise *down* as
+the stencil grows, while `order=1` barely moves. So the noise cannot be bought
+off with more neighbours.
+
+Measured end to end, with particles assigned material 1 with probability
+$p(x)=x$ so that the exact nodal fraction is a known linear function
+(three seeds, `cellSize` 1/16 and 1/32):
+
+| scheme | bias | rms | level-set range |
+|---|---|---|---|
+| shipped `update_type=0` (scatter) | ~1e-3 | **0.09 – 0.11** | `[0, 1]` |
+| shipped `update_type=1` (gather) | ~1e-3 | 0.18 | `[0, 1]` |
+| gather, `order=1` | ~1e-3 | 0.18 – 0.20 | **`[-0.35, 1.32]`** |
 
-The accuracy result is a wash — `order=1` is better at the coarse resolution
-and equal or worse at the fine one, with the ordering flipping between cases —
-while `nnn=8` violates the `[0, 1]` bound by ~3.8%.
+The bias that linear exactness would remove is already ~1e-3 for every scheme,
+because a roughly symmetric stencil reproduces a linear ramp in expectation
+anyway. What is left is variance — and `order=1` amplifies it and breaks the
+`[0, 1]` bound that `constitutive_models.py` relies on.
 
-Partition of unity survives either way (all indices share one weight set, and
-the indicator flags sum to one per particle, so the level sets sum to
-`Σ w_j = 1` regardless of sign). But a *negative* material fraction is still
-physically wrong, and `constitutive_models.py` consumes these directly.
+Partition of unity survives regardless of weight sign (all indices share one
+denominator, and the per-particle indicator flags sum to one, so the level sets
+sum to $\sum_j w_j = 1$). It was never the objection.
 
 ```{note}
-The interface metric groups nodes into rows by `y` and interpolates the 0.5
-crossing, which is crude on an unstructured simplex mesh — the *maximum* error
-is identical across all schemes because it is set by node spacing, not by the
-weights. Only the median is informative, and it is the median that shows no
-consistent gain.
+The two branches are different algorithms, and only one has a stencil to
+re-weight. `update_type=0`, the **default**, is a *scatter*: it queries `nnn`
+nodes but masks with `is_nearest`, so each particle accumulates $1/d$ into its
+nearest node only. `update_type=1` is the gather form. That the scatter has the
+lowest rms is consistent with the variance argument — a node aggregates every
+particle that is nearest to it, typically many more than `nnn`.
 ```
 
-So the swarm story is deliberately split: the plain `SwarmVariable` proxy takes
-`order=1` because its fields are smooth and the gain is two orders of magnitude;
-`IndexSwarmVariable` keeps inverse distance because its field is a jump, where
-there is nothing to gain and a bound to lose.
+So the lever for material-fraction accuracy is **more samples per node**, not
+better polynomial reproduction. The swarm story splits accordingly: the plain
+`SwarmVariable` proxy takes `order=1` because it carries exact real values and
+only bias matters; `IndexSwarmVariable` keeps inverse distance because it
+carries quantised values and variance dominates.
 
 Consumers that depend on absolute boundedness, and are therefore deliberately
 left on `order=0`: