diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 1cdbc5dd6..39caa892e 100644 --- a/docs/developer/CHANGELOG.md +++ b/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). diff --git a/docs/developer/subsystems/interpolation.md b/docs/developer/subsystems/interpolation.md index a187be062..261fb9622 100644 --- a/docs/developer/subsystems/interpolation.md +++ b/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`: