From 5dc6a29fe279131d4cd5059078e4340b1ebd5645 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 12:02:59 +1000 Subject: [PATCH] docs: correct why IndexSwarmVariable level sets keep inverse distance The recorded reason was wrong. It argued that a material indicator is piecewise constant, so linear exactness has nothing to gain away from an interface and cannot help at one. That conflates two fields: the indicator is the PARTICLE data, but the quantity estimated at a node is the local material FRACTION, which near an interface is a smooth ramp -- precisely what a constants-only scheme cannot reproduce. The supporting measurement was also unsound: it applied a gather-form variant, which models update_type=1, while the shipped default is update_type=0, a scatter (each particle accumulates 1/d into its nearest node, masked by is_nearest). There is no stencil there to re-weight, so the branches have to be measured separately. Its interface metric grouped unstructured nodes into rows by y, and its maximum was set by node spacing rather than by the weights. Re-measured against a known analytic fraction: particles take material 1 with probability p(x) = x, so the exact nodal value is p(x_node). scheme bias rms 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 conclusion stands but the mechanism is variance, not bias. A node estimates the fraction from a handful of INTEGER samples, so Var(estimate) = sum_j w_j^2 Var(f_j). Uniform weights minimise sum(w^2) at 1/nnn; inverse distance stays near that floor while linear-exact weights, being signed, sit 6-12x above it -- and, unlike inverse distance, do not come down as the stencil grows: dim nnn sum w^2 IDW 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 The bias linear exactness would remove is already ~1e-3 for every scheme, because a roughly symmetric stencil reproduces a linear ramp in expectation. Partition of unity holds for any weight sign (shared denominator, indicator flags sum to one per particle) and was never the objection. Actionable consequence: the lever for material-fraction accuracy is more samples per node, not better polynomial reproduction -- consistent with the scatter branch, which aggregates every particle nearest to a node, having the lowest rms. Script and output: ~/+Simulations/linear-rbf-proxy/step4_*. Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 9 +-- docs/developer/subsystems/interpolation.md | 83 ++++++++++++++-------- 2 files changed, 58 insertions(+), 34 deletions(-) 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`: