diff --git a/docs/developer/index.md b/docs/developer/index.md index 3e49e380..977e09ac 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -35,6 +35,7 @@ same topic are reference or historical material subordinate to the governing doc |-------|--------------------| | Coding style & API conventions | [UW3 Style Charter](UW3_STYLE_CHARTER.md) (detailed reference: [Style Guide](UW3_Style_and_Patterns_Guide.md)) | | Data access | [subsystems/data-access.md](subsystems/data-access.md) (internals reference: [NDArray System](UW3_Developers_NDArrays.md)) | +| Local scattered-point interpolation | [subsystems/interpolation.md](subsystems/interpolation.md) | | Units | [design/UNITS_SIMPLIFIED_DESIGN_2025-11.md](design/UNITS_SIMPLIFIED_DESIGN_2025-11.md) | | Testing tiers | [TESTING-RELIABILITY-SYSTEM.md](TESTING-RELIABILITY-SYSTEM.md) | | Branching & releases | [guides/branching-strategy.md](guides/branching-strategy.md) | @@ -189,6 +190,7 @@ subsystems/constitutive-models-theory subsystems/constitutive-models-anisotropy subsystems/swarm-system subsystems/data-access +subsystems/interpolation subsystems/expressions-functions subsystems/containers subsystems/checkpointing-system diff --git a/docs/developer/subsystems/data-access.md b/docs/developer/subsystems/data-access.md index cd1ef76a..764d5731 100644 --- a/docs/developer/subsystems/data-access.md +++ b/docs/developer/subsystems/data-access.md @@ -103,6 +103,10 @@ symbolic_expr = swarm_var.sym # Triggers RBF interpolation if stale This pattern prevents circular dependencies when callbacks would trigger nested PETSc field access. +The interpolation itself — which fields the proxy transfer reproduces exactly, +how to choose `nnn`, and when boundedness matters more than accuracy — is +covered in [Local Scattered-Point Interpolation](interpolation.md). + ### Vector Availability Variables set `_available=True` by default, ensuring solvers can access vectors without modification. Lazy initialization creates vectors on first access: diff --git a/docs/developer/subsystems/interpolation.md b/docs/developer/subsystems/interpolation.md new file mode 100644 index 00000000..6dc77eae --- /dev/null +++ b/docs/developer/subsystems/interpolation.md @@ -0,0 +1,320 @@ +--- +title: "Local Scattered-Point Interpolation" +--- + +## Purpose + +Underworld needs to move field values between point sets that have no shared +mesh topology: particles to proxy mesh nodes, nodes to arbitrary evaluation +points, one mesh to another after adaptation. The tool for that is a *local* +interpolant — for each target point, take its `nnn` nearest source points from +a kd-tree and form a weighted average. + +This is general-purpose numerical machinery, and the goal is straightforward: +**the interpolation should be accurate**. The single question that determines +that is **which fields the weights reproduce exactly**. A scheme that cannot +reproduce a uniform gradient will smear every field that has one, everywhere it +is used, and no amount of tuning elsewhere recovers what it threw away. + +That is what `order` selects on +{py:meth}`underworld3.ckdtree.KDTree.rbf_interpolator_local`. + +## The two schemes + +### `order=0` — inverse distance (Shepard) + +$$ w_j = \frac{d_j^{-p}}{\sum_k d_k^{-p}} $$ + +Weights are positive and sum to one. Consequences, both of them important: + +- A **constant** field is reproduced exactly. +- The result is a **convex combination** of the neighbouring values, so it can + never overshoot them. The interpolant is intrinsically bounded. +- A field with a **gradient is smeared**, and the error does *not* vanish as + the stencil tightens. Adding neighbours does not help; refining the point + spacing helps only at first order. + +### `order=1` — polyharmonic RBF with an affine tail + +For a target $x^*$ with neighbours $x_j$, solve the small saddle-point system + +$$ +\begin{bmatrix} A & P \\ P^{T} & 0 \end{bmatrix} +\begin{bmatrix} w \\ \lambda \end{bmatrix} += +\begin{bmatrix} \varphi(|x^* - x_j|) \\ 1,\; x^* \end{bmatrix} +$$ + +with $A_{ij} = \varphi(|x_i - x_j|)$, the polyharmonic (thin-plate) kernel +$\varphi(r) = r^2 \log r$, and the affine tail $P = [1,\; x_j]$. + +The lower block *is* the reproduction constraint $P^{T} w = [1,\; x^*]$, so + +$$ \sum_j w_j = 1, \qquad \sum_j w_j x_j = x^* $$ + +hold to round-off by construction: **constants and linear fields are exact**. + +Weights are signed, so the interpolant can overshoot the neighbouring values — +the price of giving up the convex combination. + +Cost is one dense $(nnn + d + 1)^3$ solve per target point, and the result is +still sparse: `nnn` non-zeros per row. The weights depend only on geometry, so +one solve serves every data component. + +## Measured + +Relative error of a swarm proxy variable against the analytic field, at the +proxy's own nodal coordinates. A **linear** field lies exactly inside both the +P1 and P2 proxy spaces, so the discretisation contributes nothing and every bit +of the measured error is particle-to-node transfer error. + +2D simplex box, `cellSize=1/8`, `fill_param=4`, proxy degree 1: + +| field | `order=0`, `nnn=3` | `order=1`, `nnn=6` | +|---|---|---| +| linear, max | 4.6e-3 | **5.1e-16** | +| quadratic, max | 9.2e-3 | 1.0e-4 | +| quadratic, rms | 2.1e-3 | 2.1e-5 | + +3D simplex box, `cellSize=1/8`, `fill_param=4`, proxy degree 1: + +| field | `order=0`, `nnn=4` | `order=1`, `nnn=8` | +|---|---|---| +| linear, max | 4.8e-3 | **5.5e-16** | +| quadratic, max | 8.5e-3 | 2.4e-4 | +| quadratic, rms | 2.0e-3 | 2.3e-5 | + +Two things to read off. First, the linear-field column is the clean statement: +one scheme is exact, the other is not. Second — and this is the part that makes +it a general accuracy improvement rather than a special case — the gain carries +over to a field with curvature, by roughly two orders of magnitude in rms. + +Widening the `order=0` stencil does **not** close the gap. At `nnn=12` the +inverse-distance error is unchanged to two significant figures, because the +error is set by the asymmetry of the neighbourhood, not by how many points are +in it. + +## Getting the operator, not just the values + +When the same transfer is applied more than once — or when the operator *is* +the product, as for a multigrid prolongation — build it once: + +```python +T = kdt.interpolation_matrix(target_coords, order=1) # scipy CSR +values = T @ data # identical to the value API +``` + +The weights depend only on geometry, so one build serves every field and every +component; the value API re-solves each call. `T @ data` is asserted equal to +`rbf_interpolator_local(...)` in the test suite so the two cannot drift. + +```{warning} +`rbf_stencil.linear_exact_weights` is the raw kernel underneath, and it returns +**zeros** for degenerate rows plus a mask. It has no kd-tree, so it can neither +widen a stencil nor fall back. A caller that ignores the mask silently +interpolates to zero. Use `interpolation_matrix` unless you specifically need +the bare weights — its rows are never empty. +``` + +### Row-wise construction carries no column guarantee + +Every row has `nnn` non-zeros. **Nothing guarantees every column is non-empty.** +A source point that is not among the `nnn` nearest neighbours of any target +produces an empty column, and a consumer forming a Galerkin coarse operator +$P^{T} A P$ then gets a singular matrix — the failure mode of issue #424. + +This is a property of *any* row-wise kNN construction, not of this kernel, and +it is not fixed by the scheme being linear-exact. Consumers that need full +column rank must check for and repair empty columns themselves. + +## Choosing `nnn` + +`order=1` requires **`nnn >= dim + 2`** and raises `ValueError` below that. + +The reason is worth knowing. At exactly `nnn == dim + 1` the affine block $P$ +is square, so the constraint $P^{T} w = [1, x^*]$ alone determines $w$: the RBF +block is inert and the scheme collapses to bare **barycentric interpolation on +the neighbour simplex**. That is singular whenever those `dim + 1` points are +collinear (2D) or coplanar (3D) — common near boundaries and on graded meshes. + +The swarm proxy default is `nnn = 2 * (dim + 1)` (6 in 2D, 8 in 3D), which +leaves enough slack that degenerate neighbourhoods are rare. That is also what +`nnn=None` resolves to at `order=1`; `order=0` keeps its historical default of +4. The default has to depend on `order` *and* dimension — a fixed default of 4 +is below the enforced minimum in 3D and would simply raise. + +## Things that fail silently if you get them wrong + +- **Build kd-trees from `.coords_nd`, never `.coords`.** `MeshVariable.coords` + dimensionalises once the model has reference quantities set, while particle + coordinates and `MeshVariable._get_kdtree` are non-dimensional. Mixing them + raises from `_convert_coords_to_tree_units` (issue #426). +- **The output carries no units.** Coordinates are converted into the tree's + frame, but `data` passes through untouched and the return is a plain array. + Re-attach units at the boundary if the caller needs them. +- **`order=1` costs about 5x `order=0`** — a dense $(nnn+d+1)^3$ solve per + target point. Irrelevant if you build an operator once; think twice in a hot + loop. +- **Do not validate against `uw.function.evaluate`.** It returns wrong values + for P1 fields at points lying exactly on cell edges in 3D (issue #432), and + proxy nodes sit on cell boundaries routinely. Compare against an analytic + field, or against an independent implementation, instead. + +## Determinism + +The weights are a pure function of the geometry: no global state, no random +numbers, no accumulation across calls. Repeated calls, and separate `KDTree` +instances built from equal point sets, return bit-identical results. + +Where several source points are exactly equidistant from a target, *which* of +them the kNN search selects is unspecified — but the choice is deterministic, +so results reproduce run to run. Both properties are pinned by tests. + +## Degenerate stencils + +If the neighbours cannot support an affine fit, the saddle system is singular. +This is not hypothetical: it happens at boundaries and in graded meshes. + +The implementation (`underworld3.utilities.rbf_stencil.linear_exact_weights`): + +1. **Pre-screens** the affine block by SVD — degenerate if + $\sigma_{\min} < 10^{-8}\,\sigma_{\max}$. +2. **Solves** only the healthy stencils, batched and chunked over targets. +3. **Validates the answer** rather than guessing a condition number: it checks + the two reproduction identities and finiteness on the computed weights, and + demotes anything that fails. +4. **Falls back** to the inverse-distance weights for those points — less + accurate there, but finite and bounded. It never returns NaN. +5. **Warns once** with the count and fraction that fell back. A silent + geometric fallback is the failure mode behind issue #424; a per-point + warning would be unusable, a single counted one is not. + +## Limiting: `monotone` + +`order=1` weights are signed, so the interpolant can oscillate where the data +is rough. `monotone=True` (or `"clamp"`) bounds that. + +The important part is *what* it bounds. The obvious limiter — clip the result +to the min/max of the stencil's source values — is **wrong for a linear-exact +scheme**, and wrong in a way that is easy to miss: + +```{warning} +A target that lies outside the convex hull of its own `nnn` neighbours has a +value outside their range **even for an exactly linear field**. Proxy nodes are +routinely in that position, and not only at domain boundaries. A raw min/max +clip therefore cannot distinguish legitimate extrapolation from ringing, and +fires on the linear part — destroying the one guarantee the scheme exists to +provide. +``` + +So the limiter here follows the slope-limiter discipline instead: **never limit +the linear reconstruction, limit only the correction on top of it.** + +1. Fit an affine function to the stencil data by least squares. +2. Keep that trend at the target untouched. +3. Bound the remaining RBF correction to the range of the non-affine residual + the stencil actually exhibits. + +Consequences, both measured: + +- On any field the scheme already reproduces exactly, the limiter is a **no-op** + (it moves the answer by ~1e-15). Linear reproduction survives it. +- On a rough field it does bite — and, like any limiter, it trades accuracy for + boundedness where the correction genuinely exceeds the observed residual. On + a `sin(6x) + |x|²` test field it moved the answer by 2e-3 (2D) and 5e-2 (3D); + in 3D that made the max error slightly worse (9.0e-2 → 1.3e-1), still well + below inverse distance (2.7e-1). + +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 + +`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**. + +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. + +Measured on a straight interface at `x = 0.5` (exactly representable, so any +displacement of the recovered 0.5 contour is scheme error): + +| 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]`** | + +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%. + +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. + +```{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. +``` + +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. + +Consumers that depend on absolute boundedness, and are therefore deliberately +left on `order=0`: + +- The **RBF rung of the point-location fallback ladder** in + `uw.function.evaluate` — `docs/developer/design/point-location-capability.md` + states its contract as "bounded, topology-free, honest". + `MeshVariable.rbf_interpolate` therefore keeps `order=0` as its default. +- `IndexSwarmVariable` level sets, which carry material fractions in $[0,1]$ and + are divided through by their sum in `constitutive_models.py`. That path + hand-rolls its own inverse-distance weighting and is untouched. + +## Where each is used + +| Caller | Default | Why | +|---|---|---| +| `SwarmVariable.rbf_interpolate` / `_rbf_to_meshVar` (proxy refresh) | `order=1`, `nnn=2*(dim+1)` | Accuracy; the proxy feeds integration, derivatives and the Lagrangian time derivative | +| `MeshVariable.rbf_interpolate` | `order=0` | Documented bounded contract (see above) | +| `SwarmVariable.read_timestep`, `MeshVariable.read_timestep` | `nnn=1` | Exact checkpoint round-trip, no averaging wanted | +| `IndexSwarmVariable._update_proxy_variables` | its own IDW | Level-set fractions must stay in $[0,1]$ | + +A rank holding fewer particles than the affine tail needs drops to `order=0` +automatically rather than failing the refresh. + +## Parallel behaviour + +Proxy refresh is **rank-local**: the kd-tree indexes only this rank's +particles and there is no halo exchange (SWARM-15, see +`docs/developer/design/SWARM_MODERNIZATION_DESIGN_2026-07.md` §4). A proxy node +near a partition seam therefore gathers from a one-sided neighbourhood. + +Linear exactness improves this but does not fix it. A linear-exact stencil +reproduces a linear field exactly from *any* neighbourhood, one-sided or not, +so for linear fields the seam error is zero and the proxy is np-independent +(pinned by `tests/parallel/test_0776_linear_rbf_proxy_parallel.py`). For a +field with curvature a one-sided stencil still differs from a centred one, so +np-dependence remains. The halo exchange in SWARM-15 is still the real fix. + +## Related + +- `docs/developer/subsystems/data-access.md` — lazy proxy updates +- `docs/developer/design/SWARM_MODERNIZATION_DESIGN_2026-07.md` §4 — rank-local seams +- `docs/developer/design/point-location-capability.md` — the evaluate fallback ladder +- `utilities/custom_mg.py` — a *global, dense* polyharmonic prolongation for + multigrid transfers. Same kernel and tail, but every coarse point enters every + row, so `PᵀAP` is dense. The local scheme documented here is the sparse + equivalent. diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index d8eb90b9..c52f6be5 100644 --- a/src/underworld3/ckdtree.pyx +++ b/src/underworld3/ckdtree.pyx @@ -34,6 +34,34 @@ def total_constructed(): return _total_constructed +def _normalise_monotone(monotone): + """Resolve the ``monotone`` argument to a bool, using UW3's one vocabulary. + + The canonical spelling lives with the evaluator + (``function.functions_unit_system._normalize_monotone``) and is reused here + so there is a single definition of what the word accepts. Only the + ``"clamp"`` mode has meaning for a local stencil: ``"pick"`` re-evaluates + out-of-bounds points through the FE path, which a kd-tree knows nothing + about. + """ + # Early out before the import: the overwhelmingly common call has no + # limiter, and ckdtree is a low-level module that should not take a + # dependency on the evaluator just to be told "no". + if monotone is False or monotone is None: + return False + + from underworld3.function.functions_unit_system import _normalize_monotone + + mode = _normalize_monotone(monotone) + if mode == "pick": + raise ValueError( + "monotone='pick' has no meaning for a local kd-tree stencil — it " + "re-evaluates out-of-bounds points through the finite-element path. " + "Use monotone='clamp' here, or uw.function.evaluate(..., monotone='pick')." + ) + return mode == "clamp" + + cdef class KDTree: """ Unit-aware KD-Tree for spatial indexing and queries. @@ -405,17 +433,32 @@ cdef class KDTree: def rbf_interpolator_local(self, coords, data, - nnn = 4, + nnn = None, p=2, verbose = False, + order = 0, + monotone = False, ): """ - Interpolate data to target coordinates using inverse distance weighting. + Interpolate data from the KD-tree points to arbitrary target coordinates. + + Two local schemes are available, selected by ``order`` — the highest + degree of polynomial the weights reproduce exactly: + + ``order=0`` (default) + Inverse-distance (Shepard) weighting. Weights are positive and sum + to one, so a **constant** field is reproduced exactly and the result + is a convex combination of the neighbouring values — bounded, but a + field with a gradient is smeared and the error does not vanish as + the stencil tightens. - This is a convenience wrapper around :meth:`rbf_interpolator_local_from_kdtree`. - It performs radial basis function (RBF) interpolation using inverse distance - weighting to map known data values from the KD-tree points to arbitrary - target coordinates. + ``order=1`` + Polyharmonic RBF with an affine tail, so **constant and linear** + fields are reproduced exactly (:math:`\\sum_j w_j = 1` and + :math:`\\sum_j w_j x_j = x^*`). Weights may be negative, so the + result can overshoot the neighbouring values; see ``monotone``. + + Both are local: ``nnn`` non-zero weights per target point. Parameters ---------- @@ -426,12 +469,27 @@ cdef class KDTree: Known data values at KD-tree points. Shape should be ``(n_points,)`` or ``(n_points, n_components)``. nnn : int, optional - Number of nearest neighbours to use for interpolation (default 4). - If 1, returns raw nearest-neighbour values without distance weighting. + Number of nearest neighbours to use. Defaults to 4 for ``order=0`` + and ``2 * (dim + 1)`` for ``order=1`` — the default has to depend + on both, because ``order=1`` requires ``nnn >= dim + 2`` and a + fixed default of 4 would raise in 3D. If 1, returns raw + nearest-neighbour values without distance weighting. p : int, optional - Power index for distance weighting: ``weight = 1/distance^p`` (default 2). + Power index for distance weighting: ``weight = 1/distance^p`` + (default 2). Used by ``order=0`` only. verbose : bool, optional Print progress messages (default False). + order : int, optional + Polynomial reproduction order, 0 (default) or 1. + monotone : bool or str, optional + ``False`` (default) or ``True`` / ``"clamp"``. Limits the + **non-affine part** of the interpolant: the local least-squares + affine trend is preserved exactly, and only the RBF correction on + top of it is bounded by the correction actually present in the + stencil. This is the slope-limiter discipline — the linear + reconstruction is never clipped — so linear reproduction survives + the limiter and it is a no-op on any field the scheme already + reproduces exactly. It bounds new oscillation, not absolute range. Returns ------- @@ -444,7 +502,7 @@ cdef class KDTree: query : Find nearest neighbours without interpolation. """ return self.rbf_interpolator_local_from_kdtree( - coords, data, nnn, p, verbose, + coords, data, nnn, p, verbose, order, monotone, ) def old_rbf_interpolator_local_from_kdtree(self, @@ -576,9 +634,235 @@ cdef class KDTree: return Values - def rbf_interpolator_local_from_kdtree(self, coords, data, nnn, p, verbose): + def _resolve_nnn(self, nnn, order): + """Default stencil size. ``order=1`` needs at least ``dim + 2``, so a + fixed default cannot serve both schemes (it would raise in 3D).""" + if nnn is not None: + return nnn + return 2 * (self.ndim + 1) if order == 1 else 4 + + def _retry_degenerate_stencils(self, coords_converted, tree_points, + degenerate, nnn): + """Re-solve the failed points on a wider stencil. + + A neighbourhood that cannot support an affine fit is usually a local + accident of point placement, so reaching further normally escapes it. + Only the failed points are re-queried, so the cost is set by how many + actually failed (measured: a handful in tens of thousands). + + Returns ``(rows, indices, weights, remaining)`` — the recovered rows and + their wider stencils, plus the mask of points that failed even after + widening and so keep their inverse-distance weights. """ - Performs an inverse distance (squared) mapping of data to the target `coords`. + from underworld3.utilities.rbf_stencil import linear_exact_weights + + wide = min(4 * nnn, self.n) + if wide <= nnn: + return None, None, None, degenerate + + failed = np.flatnonzero(degenerate) + targets = np.ascontiguousarray(coords_converted[failed], dtype=np.float64) + + # The raw query: coords_converted is already in tree units, so the + # unit-aware `query` would reject it when the tree carries units. + indices, _ = self.find_closest_n_points(wide, targets) + if np.any(indices >= self.n): + return None, None, None, degenerate + + weights, still_degenerate = linear_exact_weights( + targets, tree_points[indices] + ) + + remaining = degenerate.copy() + remaining[failed] = still_degenerate + + recovered = ~still_degenerate + if not recovered.any(): + return None, None, None, remaining + + return (failed[recovered], indices[recovered], weights[recovered], + remaining) + + def _local_stencil(self, coords_converted, nnn, p, order): + """Neighbour indices and weights for every target point. + + The single source of the interpolation weights: both + :meth:`rbf_interpolator_local` and :meth:`interpolation_matrix` are + built on this, so the operator and the values it produces cannot drift + apart. + + Returns + ------- + indices, weights : ndarray + Shape ``(n_targets, nnn)``. The primary stencil. + wide : tuple or None + ``(rows, indices, weights)`` for targets re-solved on a wider + stencil after their primary one proved degenerate. Their primary + row is superseded and must be discarded by the caller. + degenerate : ndarray + Bool, shape ``(n_targets,)``. Targets still degenerate after + widening; these carry inverse-distance weights, which are bounded + but not linear-exact. + """ + # find_closest_n_points returns (indices, dist_sqr) -- the reverse of + # query()'s (dist, indices). It is used here rather than query() + # because coords_converted is already in tree units, and query() would + # convert a second time. + closest_n, distance_n = self.find_closest_n_points( + nnn, np.ascontiguousarray(coords_converted, dtype=np.float64) + ) + + # valid indices are 0..n-1; the empty-tree sentinel (0 with n=0) + # must trip this guard, so the comparison is >= (issue #399). + if np.any(closest_n >= self.n): + raise RuntimeError( + "Error in rbf_interpolator_local_from_kdtree - a nearest neighbour wasn't found" + ) + + # np.bool_, not bool: this module cimports the C++ `bool` from libcpp, + # which shadows the Python builtin and will not compile here. + degenerate = np.zeros(coords_converted.shape[0], dtype=np.bool_) + + if nnn == 1: + return closest_n, np.ones(closest_n.shape), None, degenerate + + # can decompose weighting vecotrs as IDW is a linear relationship + # build normalise weight vectors and multiply that with known data + # TODO(BUG): issue #427 — `distance_n` holds SQUARED distances, so the + # decay is r^(-2p), not the documented r^(-p), and `epsilon` floors r + # at ~1e-6 rather than 1e-12. + epsilon = 1e-12 + weights = 1 / np.power(epsilon + distance_n[:], p) + n_weights = (weights.T / np.sum(weights, axis=1)).T + + if order == 0: + return closest_n, n_weights, None, degenerate + + from underworld3.utilities.rbf_stencil import linear_exact_weights + + tree_points = np.asarray(self.points) + linear_weights, degenerate = linear_exact_weights( + coords_converted, tree_points[closest_n[:]] + ) + # A stencil that cannot support an affine fit (collinear in 2D, + # coplanar in 3D) keeps the inverse-distance weights: finite and + # bounded, but not linear-exact. + linear_weights[degenerate] = n_weights[degenerate] + + wide = None + if degenerate.any(): + # Widen before surrendering. A single non-exact node sitting among + # exact ones is an isolated SPIKE, far more damaging to derivatives + # and integrals than a smooth error of the same size. + rows, wide_indices, wide_weights, degenerate = ( + self._retry_degenerate_stencils( + coords_converted, tree_points, degenerate, nnn + ) + ) + if rows is not None: + wide = (rows, wide_indices, wide_weights) + + n_degenerate = int(degenerate.sum()) + if n_degenerate: + import warnings + + warnings.warn( + f"rbf_interpolator_local(order=1): {n_degenerate} of " + f"{degenerate.size} stencils could not support an affine fit " + "(collinear/coplanar neighbours) even after widening, and fell " + "back to inverse-distance weighting.", + stacklevel=3, + ) + + return closest_n, linear_weights, wide, degenerate + + def interpolation_matrix(self, coords, nnn=None, p=2, order=0): + """Sparse operator mapping values on the KD-tree points to ``coords``. + + ``T @ data`` is exactly what :meth:`rbf_interpolator_local` returns for + the same arguments — both are built from the same weights. + + Use this instead of the value API when the same transfer is applied + more than once, or when the operator itself is the product (a multigrid + prolongation, for example). The weights depend only on geometry, so one + build serves every field and every component; the value API re-solves + each time. + + Parameters + ---------- + coords : array-like + Target coordinates, shape ``(n_coords, dim)``. + nnn : int, optional + Neighbours per target. Defaults to 4 for ``order=0`` and + ``2 * (dim + 1)`` for ``order=1``. + p : int, optional + Inverse-distance power; used by ``order=0`` only. + order : int, optional + Polynomial reproduction order, 0 or 1. See + :meth:`rbf_interpolator_local`. + + Returns + ------- + scipy.sparse.csr_matrix + Shape ``(n_coords, self.n)``. Rows carry ``nnn`` non-zeros, except + those re-solved on a wider stencil. + + Notes + ----- + Degenerate stencils are handled exactly as in the value path — widened, + then fall back to inverse-distance weights with a warning. Rows are + never empty. + + **Row-wise construction gives no column guarantee.** Every row has + ``nnn`` non-zeros, but a source point that is not among the neighbours + of any target produces an empty column. Consumers that need full column + rank -- a Galerkin coarse operator :math:`P^{T} A P`, for instance -- + must check for and repair empty columns themselves (see issue #424). + """ + import scipy.sparse as sp + + if order not in (0, 1): + raise ValueError( + f"order must be 0 (inverse distance) or 1 (linear-exact), got {order!r}." + ) + coords_converted = self._convert_coords_to_tree_units(coords) + if coords_converted.shape[1] != self.ndim: + raise RuntimeError( + f"Interpolation coordinates dimensionality " + f"({coords_converted.shape[1]}) is different to kD-tree " + f"dimensionality ({self.ndim})." + ) + nnn = self._resolve_nnn(nnn, order) + + indices, weights, wide, _ = self._local_stencil( + coords_converted, nnn, p, order + ) + + n_targets = coords_converted.shape[0] + keep = np.ones(n_targets, dtype=np.bool_) + if wide is not None: + keep[wide[0]] = False + + rows = np.repeat(np.flatnonzero(keep), nnn) + cols = indices[keep].ravel() + vals = weights[keep].ravel() + + if wide is not None: + wide_rows, wide_indices, wide_weights = wide + rows = np.concatenate( + [rows, np.repeat(wide_rows, wide_indices.shape[1])] + ) + cols = np.concatenate([cols, wide_indices.ravel()]) + vals = np.concatenate([vals, wide_weights.ravel()]) + + return sp.csr_matrix( + (vals, (rows, cols.astype(np.int64))), shape=(n_targets, self.n) + ) + + def rbf_interpolator_local_from_kdtree(self, coords, data, nnn, p, verbose, + order=0, monotone=False): + """ + Map data held on the KD-tree points onto the target `coords`. This method is unit-aware: if the KD-tree was built with unit-aware coordinates, it will automatically convert query coordinates to match before interpolation. @@ -596,15 +880,29 @@ cdef class KDTree: nnn : int The number of neighbour points to sample from. If `1`, no distance averaging is done. p : int - The power index to calculate weights, i.e., pow(distance, -p) + The power index to calculate weights, i.e., pow(distance, -p). + Used by ``order=0`` only. verbose : bool Print when mapping occurs + order : int, optional + Polynomial reproduction order: 0 for inverse-distance (constants + exact), 1 for polyharmonic + affine tail (constants and linears + exact). See :meth:`rbf_interpolator_local`. + monotone : bool or str, optional + ``False``, ``True`` or ``"clamp"`` — bound each result to the + min/max of its own stencil's source values. Returns ------- ndarray Interpolated data values at target coordinates """ + if order not in (0, 1): + raise ValueError( + f"order must be 0 (inverse distance) or 1 (linear-exact), got {order!r}." + ) + monotone_clamp = _normalise_monotone(monotone) + # Convert coordinates to match tree's coordinate system coords_converted = self._convert_coords_to_tree_units(coords) @@ -617,39 +915,75 @@ cdef class KDTree: f"Data does not match kd-tree size array ({data.shape[0]} v ({self.n}))" ) - coords_contiguous = np.ascontiguousarray(coords_converted) - # query nnn points to the coords - # distance_n is a list of distance to the nearest neighbours for all coords_contiguous - # closest_n is the index of the neighbours from ncoords for all coords_contiguous - # Note: query() returns sqr_dists=True by default, and we use the converted coords - distance_n, closest_n = self.query(coords, k=nnn) - - # valid indices are 0..n-1; the empty-tree sentinel (0 with n=0) - # must trip this guard, so the comparison is >= (issue #399). - if np.any(closest_n >= self.n): - raise RuntimeError( - "Error in rbf_interpolator_local_from_kdtree - a nearest neighbour wasn't found" - ) + nnn = self._resolve_nnn(nnn, order) if verbose and uw.mpi.rank == 0: - # For Debugging - # print(f"kd-tree diagnostics: d.shape - {distance_n.shape}, c.shape - {closest_n.shape}") print(f"Mapping values with nnn - {nnn} & p {p} ... start", flush=True) if nnn == 1: # only use nearest neighbour raw data + if order == 1: + raise ValueError( + "order=1 needs at least dim + 2 neighbours to determine the " + f"affine tail; nnn=1 selects the raw nearest-neighbour path." + ) + closest_n, _ = self.find_closest_n_points( + 1, np.ascontiguousarray(coords_converted, dtype=np.float64) + ) + # (n, 1) -> (n,): the nearest-neighbour path returns data rows + # directly, so the stencil axis must not survive into the result. + # query(k=1) used to do this reshape for us. + closest_n = closest_n.reshape(-1) + if np.any(closest_n >= self.n): + raise RuntimeError( + "Error in rbf_interpolator_local_from_kdtree - a nearest neighbour wasn't found" + ) return data[closest_n] - # can decompose weighting vecotrs as IDW is a linear relationship - # build normalise weight vectors and multiply that with known data - epsilon = 1e-12 - weights = 1 / np.power(epsilon + distance_n[:], p) - n_weights = (weights.T / np.sum(weights, axis=1)).T + closest_n, n_weights, wide, degenerate = self._local_stencil( + coords_converted, nnn, p, order + ) kdata = data[closest_n[:]] # magic with einstein summation power vals = np.einsum("sdc,sd->sc", kdata, n_weights) - # print(valz) + + if wide is not None: + wide_rows, wide_indices, wide_weights = wide + vals[wide_rows] = np.einsum( + "sdc,sd->sc", data[wide_indices], wide_weights + ) + + if monotone_clamp: + if order == 0: + # Already a convex combination; the clip is exact-arithmetic + # redundant and only guards round-off. + vals = np.clip(vals, kdata.min(axis=1), kdata.max(axis=1)) + else: + # Limit the CORRECTION, never the linear part -- the same + # discipline as a slope limiter in a second-order FV scheme. + # + # Clipping the total against the stencil's raw min/max would + # be wrong here: a target outside the convex hull of its own + # neighbours has a value outside their range even for an + # exactly linear field, so the naive clip cannot distinguish + # legitimate extrapolation from ringing and destroys the + # reproduction guarantee. + from underworld3.utilities.rbf_stencil import affine_trend + + # The bound comes from the primary stencil even for the few + # rows re-solved on a wider one: it is still a legitimate local + # bound, and the affine trend -- the part being preserved -- is + # what matters here. + trend_at_target, trend_at_stencil = affine_trend( + coords_converted, np.asarray(self.points)[closest_n[:]], kdata + ) + residual = kdata - trend_at_stencil + vals = trend_at_target + np.clip( + vals - trend_at_target, + residual.min(axis=1), + residual.max(axis=1), + ) if verbose and uw.mpi.rank == 0: print(f"Mapping values ... finished", flush=True) diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 1a71d2c4..3e6bdabc 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -886,24 +886,41 @@ def _get_kdtree(self): return self._kdtree - def rbf_interpolate(self, new_coords, meth=0, p=2, verbose=False, nnn=None, rubbish=None): + def rbf_interpolate(self, new_coords, meth=0, p=2, verbose=False, nnn=None, + rubbish=None, order=0, monotone=False): """Interpolate variable data to new coordinates using RBF. Uses inverse distance weighting with k-nearest neighbors to interpolate values from mesh nodes to arbitrary coordinates. + The default ``order=0`` is deliberate and differs from the swarm proxy + path: this method is the RBF rung of the point-location fallback ladder + in :func:`underworld3.function.evaluate`, whose documented contract is + that it is *bounded*. Inverse-distance weights are a convex combination + and so cannot overshoot; ``order=1`` weights can. Pass ``order=1`` + explicitly where linear exactness matters more than boundedness. + Parameters ---------- new_coords : numpy.ndarray Target coordinates of shape ``(n_points, dim)``. meth : int, optional Interpolation method (reserved, currently unused). + TODO(BUG): issue #428 — ``meth`` and ``rubbish`` are dead + parameters, and ``tests/test_0505_rbf_swarm_mesh.py`` passes its + ``nnn`` into ``meth`` positionally, so that test silently + discards it. p : float, optional Power parameter for inverse distance weighting (default: 2). verbose : bool, optional Print progress information. nnn : int, optional Number of nearest neighbors (default: 4 for 3D, 3 for 2D). + order : int, optional + Polynomial reproduction order, 0 (default, bounded) or 1 + (constants and linears exact; requires ``nnn >= dim + 2``). + monotone : bool or str, optional + Bound each value to the min/max of its own stencil. Returns ------- @@ -928,7 +945,9 @@ def rbf_interpolate(self, new_coords, meth=0, p=2, verbose=False, nnn=None, rubb # Use cached KDTree for interpolation kdt = self._get_kdtree() - values = kdt.rbf_interpolator_local(new_coords, D, nnn, p=p, verbose=verbose) + values = kdt.rbf_interpolator_local( + new_coords, D, nnn, p=p, verbose=verbose, order=order, monotone=monotone + ) return values diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 2c22be7d..8ef5fab8 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -1140,9 +1140,15 @@ def _update_proxy_if_stale(self): return # Maybe rbf_interpolate for this one and meshVar is a special case - def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False): + def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False, order=1, + monotone=False): """ - Here is how it works: for each particle, create a distance-weighted average on the node data + Refresh a proxy mesh variable from the particles. + + Each proxy node gathers from its ``nnn`` nearest particles. The + default weights reproduce linear fields exactly (``order=1``), so a + field with a uniform gradient transfers without smearing; ``nnn`` and + ``order`` are resolved in :meth:`rbf_interpolate`. Todo: caching the k-d trees etc for the proxy-mesh-variable nodal points Todo: some form of global fall-back for when there are no particles on a processor @@ -1151,9 +1157,6 @@ def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False): # Mapping to the coordinates of the variable from the # particle coords - if nnn is None: - nnn = self.swarm.mesh.dim + 1 - if meshVar.mesh != self.swarm.mesh: # If this is our own proxy variable and mesh has changed, recreate it if hasattr(self, "_meshVar") and meshVar is self._meshVar: @@ -1194,7 +1197,9 @@ def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False): ) Values = current_values else: - Values = self.rbf_interpolate(new_coords, verbose=verbose, nnn=nnn) + Values = self.rbf_interpolate( + new_coords, verbose=verbose, nnn=nnn, order=order, monotone=monotone + ) meshVar.data[...] = Values[...] @@ -1483,12 +1488,16 @@ def _object_viewer(self): display(self.data), return - def rbf_interpolate(self, new_coords, verbose=False, nnn=None): + def rbf_interpolate(self, new_coords, verbose=False, nnn=None, order=1, + monotone=False): """ Radial basis function interpolation of particle data to arbitrary points. - Uses inverse-distance weighting to interpolate particle values - to new coordinate locations. + By default this reproduces constant *and linear* fields exactly + (``order=1``): a polyharmonic kernel with an affine tail over the + ``nnn`` nearest particles. Inverse-distance weighting (``order=0``) + reproduces only constants, so any field with a gradient is smeared by + an error that does not vanish as the particles crowd together. Parameters ---------- @@ -1497,7 +1506,18 @@ def rbf_interpolate(self, new_coords, verbose=False, nnn=None): verbose : bool, default=False Print diagnostic information during interpolation. nnn : int, optional - Number of nearest neighbors to use. Defaults to ``mesh.dim + 1``. + Number of nearest neighbours to use. Defaults to + ``2 * (mesh.dim + 1)`` — comfortably above the ``dim + 2`` that the + affine tail needs, so that near-degenerate particle neighbourhoods + do not have to fall back. + order : int, default=1 + Polynomial reproduction order: 1 (constants and linears exact) or + 0 (constants only, inverse distance). Drops to 0 automatically on a + rank holding too few particles to determine the affine tail. + monotone : bool or str, default=False + Limit the non-affine part of the interpolant to the non-affine + variation present in the particle stencil. The local linear trend + is preserved, so this does not cost linear exactness. Returns ------- @@ -1532,15 +1552,23 @@ def rbf_interpolate(self, new_coords, verbose=False, nnn=None): return np.zeros((new_coords.shape[0], data_size[1])) if nnn is None: - nnn = self.swarm.mesh.dim + 1 + nnn = 2 * (self.swarm.mesh.dim + 1) if nnn > data_size[0]: nnn = data_size[0] + # A rank holding fewer particles than the affine tail needs cannot + # support a linear fit at all. Inverse distance still gives a sensible + # answer there, so degrade rather than fail the whole refresh. + if order == 1 and nnn < self.swarm.mesh.dim + 2: + order = 0 + # Use direct PETSc access to avoid callback circular dependency D = raw_data.copy() kdt = self.swarm._get_kdtree() - values = kdt.rbf_interpolator_local(new_coords, D, nnn, 2, verbose) + values = kdt.rbf_interpolator_local( + new_coords, D, nnn, 2, verbose, order=order, monotone=monotone + ) return values diff --git a/src/underworld3/utilities/rbf_stencil.py b/src/underworld3/utilities/rbf_stencil.py new file mode 100644 index 00000000..20f4acf4 --- /dev/null +++ b/src/underworld3/utilities/rbf_stencil.py @@ -0,0 +1,260 @@ +"""Local RBF-FD stencil weights that reproduce linear fields exactly. + +Inverse-distance (Shepard) weights are positive and sum to one, so they +reproduce a *constant* exactly but not a linear field: any field with a +gradient is smeared, and the error does not vanish as the stencil tightens. + +The weights built here do reproduce linears. For a target point +:math:`x^*` and its :math:`m` nearest source points :math:`x_j`, solve the +small saddle-point system + +.. math:: + + \\begin{bmatrix} A & P \\\\ P^{T} & 0 \\end{bmatrix} + \\begin{bmatrix} w \\\\ \\lambda \\end{bmatrix} + = \\begin{bmatrix} \\varphi(|x^* - x_j|) \\\\ 1,\\; x^* \\end{bmatrix} + +with :math:`A_{ij} = \\varphi(|x_i - x_j|)`, the polyharmonic (thin-plate) +kernel :math:`\\varphi(r) = r^2 \\log r`, and the affine tail +:math:`P = [1,\\; x_j]`. The lower block *is* the reproduction constraint +:math:`P^{T} w = [1,\\; x^*]`, so + +.. math:: + + \\sum_j w_j = 1, \\qquad \\sum_j w_j x_j = x^* + +hold to round-off by construction — constants and linears are exact. The +weights depend only on geometry, so one solve serves every data component. + +Cost is one dense :math:`(m + d + 1)^3` solve per target point, and the +result is sparse: :math:`m` non-zeros per row. + +This module is deliberately pure NumPy with no PETSc dependency; it is the +numeric core behind ``KDTree.rbf_interpolator_local(..., order=1)``. +""" + +import numpy as np + +# The kernel is not scale-invariant: under r -> s.r it gains an r^2 log(s) +# term, which is quadratic and therefore NOT absorbed by the affine tail. So +# every stencil is solved in a local frame (target at the origin, farthest +# neighbour at unit radius); conditioning is then independent of the absolute +# stencil size, and both tolerances below are dimensionless. +_RANK_TOL = 1.0e-8 # smallest/largest singular value of the affine block +_REPRODUCTION_TOL = 1.0e-6 # residual of the two reproduction identities + +# Peak working set per chunk of the stacked saddle matrices. Sized for memory, +# not speed: the full (N, m+d+1, m+d+1) stack for a million targets would be +# gigabytes. +_CHUNK_BYTES = 64 * 1024 * 1024 + + +def affine_trend(target_coords, neighbour_coords, values): + """Least-squares affine fit of ``values`` over each stencil. + + Returns the fit evaluated at the target and at each neighbour, which is + what a limiter needs in order to bound the *non-affine* part of an + interpolant without touching its linear part. + + Parameters + ---------- + target_coords : numpy.ndarray + Shape ``(n_targets, dim)``. + neighbour_coords : numpy.ndarray + Shape ``(n_targets, nnn, dim)``. + values : numpy.ndarray + Stencil data, shape ``(n_targets, nnn, n_components)``. + + Returns + ------- + at_target : numpy.ndarray + Shape ``(n_targets, n_components)``. + at_neighbours : numpy.ndarray + Shape ``(n_targets, nnn, n_components)``. + """ + target_coords = np.ascontiguousarray(target_coords, dtype=np.float64) + neighbour_coords = np.ascontiguousarray(neighbour_coords, dtype=np.float64) + + n_targets, nnn, _ = neighbour_coords.shape + + # Local frame again, so the constant column and the coordinate columns are + # comparably scaled and the fit is not dominated by the offset. + offsets = neighbour_coords - target_coords[:, None, :] + radius = np.linalg.norm(offsets, axis=2).max(axis=1) + y = offsets / np.where(radius > 0.0, radius, 1.0)[:, None, None] + + P = np.concatenate([np.ones((n_targets, nnn, 1)), y], axis=2) + # pinv rather than solve: a rank-deficient stencil gives the minimum-norm + # fit instead of raising, and those points are limited conservatively. + coefficients = np.linalg.pinv(P) @ values + + # The target sits at the origin of the local frame, so its affine value is + # the constant coefficient. + return coefficients[:, 0, :], P @ coefficients + + +def _polyharmonic(r): + """:math:`\\varphi(r) = r^2 \\log r`, with :math:`\\varphi(0) = 0`.""" + # The clip only keeps log() finite at coincident points; r**2 drives the + # product to zero there regardless of the value substituted. + safe = np.where(r > 0.0, r, 1.0) + return np.where(r > 0.0, r ** 2 * np.log(safe), 0.0) + + +def _solve_stencils(matrices, rhs): + """Batched saddle-point solve, tolerant of an exactly singular member. + + ``np.linalg.solve`` raises for the whole batch if any one matrix is + singular, and gives no way to tell which. The pre-screen in + :func:`linear_exact_weights` removes the common case; this catches the + residue by falling back to the pseudo-inverse for that chunk, which never + raises. The post-solve reproduction check then rejects anything the + pseudo-inverse merely papered over. + """ + try: + return np.linalg.solve(matrices, rhs[..., None])[..., 0] + except np.linalg.LinAlgError: + return (np.linalg.pinv(matrices) @ rhs[..., None])[..., 0] + + +def linear_exact_weights(target_coords, neighbour_coords): + """Stencil weights that reproduce constant and linear fields exactly. + + Parameters + ---------- + target_coords : numpy.ndarray + Points to interpolate to, shape ``(n_targets, dim)``. + neighbour_coords : numpy.ndarray + Source points of each target's stencil, shape + ``(n_targets, nnn, dim)``. Normally the ``nnn`` nearest neighbours + from a kd-tree query. + + Returns + ------- + weights : numpy.ndarray + Shape ``(n_targets, nnn)``. Interpolate with + ``(weights[:, :, None] * data[stencil]).sum(axis=1)``. + degenerate : numpy.ndarray + Boolean, shape ``(n_targets,)``. True where the stencil could not + support an affine fit — collinear neighbours in 2D, coplanar in 3D, + or a stencil that collapsed onto its target. + + Warnings + -------- + **Rows flagged ``degenerate`` are returned as ZEROS.** A caller that + ignores the mask gets a silent interpolation to zero, which is worse than + an inaccurate answer. This is a raw kernel: it has no kd-tree, so it can + neither widen the stencil nor fall back to inverse distance. + + Prefer :meth:`underworld3.ckdtree.KDTree.interpolation_matrix` unless you + genuinely need the bare weights. It builds the same weights, then widens + degenerate stencils and falls back to inverse distance with a warning, so + its rows are never empty. + + Notes + ----- + ``nnn`` must be at least ``dim + 2``. At exactly ``dim + 1`` the affine + block is square, the constraint alone determines the weights, and the + scheme degenerates to bare barycentric interpolation on the neighbour + simplex — which is singular whenever those points are collinear or + coplanar, a common situation near boundaries and in graded meshes. + """ + target_coords = np.ascontiguousarray(target_coords, dtype=np.float64) + neighbour_coords = np.ascontiguousarray(neighbour_coords, dtype=np.float64) + + n_targets, nnn, dim = neighbour_coords.shape + if target_coords.shape != (n_targets, dim): + raise ValueError( + f"target_coords has shape {target_coords.shape}, expected " + f"({n_targets}, {dim}) to match neighbour_coords " + f"{neighbour_coords.shape}." + ) + if nnn < dim + 2: + raise ValueError( + f"A linear-exact stencil needs at least dim + 2 = {dim + 2} " + f"neighbours, got nnn = {nnn}. At dim + 1 the affine tail is " + "exactly determined and the scheme reduces to barycentric " + "interpolation on the neighbour simplex, which is singular for " + "collinear (2D) or coplanar (3D) neighbours." + ) + + weights = np.zeros((n_targets, nnn), dtype=np.float64) + degenerate = np.zeros(n_targets, dtype=bool) + + # Budget the two largest per-row temporaries: the stacked saddle matrices + # (size x size) and the pairwise offset array (nnn x nnn x dim) that the + # kernel distances are formed from. The latter is easy to overlook and is + # the bigger of the two in 3D. + size = nnn + dim + 1 + bytes_per_row = 8 * (size * size + nnn * nnn * dim) + rows_per_chunk = max(1, _CHUNK_BYTES // bytes_per_row) + + for start in range(0, n_targets, rows_per_chunk): + stop = min(start + rows_per_chunk, n_targets) + _weights_for_chunk( + target_coords[start:stop], + neighbour_coords[start:stop], + weights[start:stop], + degenerate[start:stop], + ) + + return weights, degenerate + + +def _weights_for_chunk(targets, neighbours, weights_out, degenerate_out): + """Fill one chunk of ``weights``/``degenerate`` in place.""" + n_chunk, nnn, dim = neighbours.shape + + # Local frame: target at the origin, farthest neighbour at unit radius. + offsets = neighbours - targets[:, None, :] + radius = np.linalg.norm(offsets, axis=2).max(axis=1) + collapsed = radius <= 0.0 # every neighbour sits on the target + scale = np.where(collapsed, 1.0, radius) + y = offsets / scale[:, None, None] + + # Affine block, and its rank as the primary degeneracy test. + P = np.concatenate([np.ones((n_chunk, nnn, 1)), y], axis=2) + singular_values = np.linalg.svd(P, compute_uv=False) + rank_deficient = singular_values[:, -1] <= _RANK_TOL * singular_values[:, 0] + + healthy = ~(collapsed | rank_deficient) + degenerate_out[...] = ~healthy + if not healthy.any(): + return + + y_h = y[healthy] + P_h = P[healthy] + n_h = y_h.shape[0] + size = nnn + dim + 1 + + pair_distance = np.linalg.norm(y_h[:, :, None, :] - y_h[:, None, :, :], axis=3) + + matrices = np.zeros((n_h, size, size), dtype=np.float64) + matrices[:, :nnn, :nnn] = _polyharmonic(pair_distance) + matrices[:, :nnn, nnn:] = P_h + matrices[:, nnn:, :nnn] = np.transpose(P_h, (0, 2, 1)) + + # RHS: kernel from the target (the origin) to each neighbour, then the + # polynomial basis evaluated at the origin, [1, 0, ..., 0]. + rhs = np.zeros((n_h, size), dtype=np.float64) + rhs[:, :nnn] = _polyharmonic(np.linalg.norm(y_h, axis=2)) + rhs[:, nnn] = 1.0 + + w = _solve_stencils(matrices, rhs)[:, :nnn] + + # Validate the answer rather than guessing a condition number: these two + # identities are what the whole construction is for. + partition_error = np.abs(w.sum(axis=1) - 1.0) + linear_error = np.linalg.norm( + np.einsum("nm,nmd->nd", w, y_h), axis=1 + ) + reproduced = ( + np.isfinite(w).all(axis=1) + & (partition_error <= _REPRODUCTION_TOL) + & (linear_error <= _REPRODUCTION_TOL) + ) + + accepted = np.zeros(n_chunk, dtype=bool) + accepted[healthy] = reproduced + weights_out[accepted] = w[reproduced] + degenerate_out[...] = ~accepted diff --git a/tests/parallel/test_0776_linear_rbf_proxy_parallel.py b/tests/parallel/test_0776_linear_rbf_proxy_parallel.py new file mode 100644 index 00000000..a82430e2 --- /dev/null +++ b/tests/parallel/test_0776_linear_rbf_proxy_parallel.py @@ -0,0 +1,101 @@ +"""Linear-exact swarm proxy transfer under MPI. + +Proxy refresh is strictly rank-local: the kd-tree indexes only this rank's +particles, and there is no halo exchange (SWARM-15, recorded in +``docs/developer/design/SWARM_MODERNIZATION_DESIGN_2026-07.md`` §4). A proxy +node near a partition seam therefore gathers from a one-sided neighbourhood, +which is why proxy values have historically been np-dependent. + +A linear-exact stencil reproduces a linear field exactly from *any* +neighbourhood, one-sided or not. So for a linear field the seam error is +zero, and the proxy becomes np-independent. That is what this test pins. + +It does **not** claim the seam problem is fixed in general: for a field with +curvature the one-sided stencil still differs from a centred one, so np +dependence remains. The test asserts the linear case only, which is the part +that is genuinely exact. + +Run with: + mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_0776_linear_rbf_proxy_parallel.py + mpirun -n 4 python -m pytest --with-mpi tests/parallel/test_0776_linear_rbf_proxy_parallel.py +""" + +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.meshing import UnstructuredSimplexBox + +pytestmark = [ + pytest.mark.level_1, + pytest.mark.tier_b, + pytest.mark.mpi(min_size=2), + pytest.mark.timeout(120), +] + + +def _linear(coords): + return 0.5 + coords @ np.arange(1, coords.shape[1] + 1, dtype=float) + + +def _build(): + mesh = UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 12.0 + ) + swarm = uw.swarm.Swarm(mesh) + var = swarm.add_variable(name="f", size=1, proxy_degree=1) + swarm.populate(fill_param=4) + return mesh, swarm, var + + +@pytest.mark.mpi(min_size=2) +def test_proxy_is_linear_exact_on_every_rank(): + """Every rank's owned proxy nodes carry the analytic linear field.""" + mesh, swarm, var = _build() + + proxy = var._meshVar + var.data[:, 0] = _linear(swarm._particle_coordinates.data) + var._rbf_to_meshVar(proxy) + + node_coords = np.asarray(proxy.coords) + expected = _linear(node_coords) + local_error = ( + np.abs(np.asarray(proxy.data[:, 0]) - expected).max() / np.abs(expected).max() + if node_coords.shape[0] + else 0.0 + ) + + # Collective: a rank-local assertion would let one bad rank pass silently + # while the others reported success (the house failure mode). + global_error = uw.mpi.comm.allreduce(local_error, uw.MPI.MAX) + + assert global_error < 1.0e-12, ( + f"rank-local linear-exact transfer left a global max error of " + f"{global_error:.3e}; seam nodes should still be exact for a linear field" + ) + + del swarm + del mesh + + +@pytest.mark.mpi(min_size=2) +def test_proxy_constant_is_exact_and_collective_refresh_completes(): + """A constant is exact under either scheme; this checks the refresh itself + stays collective (every rank performs the same read-then-write).""" + mesh, swarm, var = _build() + + proxy = var._meshVar + var.data[:, 0] = 2.75 + var._rbf_to_meshVar(proxy) + + local_error = ( + np.abs(np.asarray(proxy.data[:, 0]) - 2.75).max() + if proxy.coords.shape[0] + else 0.0 + ) + global_error = uw.mpi.comm.allreduce(local_error, uw.MPI.MAX) + + assert global_error < 1.0e-12, f"constant field error {global_error:.3e}" + + del swarm + del mesh diff --git a/tests/test_0102_kdtree_linear_rbf.py b/tests/test_0102_kdtree_linear_rbf.py new file mode 100644 index 00000000..652861a8 --- /dev/null +++ b/tests/test_0102_kdtree_linear_rbf.py @@ -0,0 +1,586 @@ +"""Local RBF stencil weights that reproduce linear fields exactly. + +``KDTree.rbf_interpolator_local(..., order=1)`` builds polyharmonic +(:math:`r^2 \\log r`) weights with an affine tail, so both + + sum_j w_j = 1 (constants exact) + sum_j w_j x_j = x* (linears exact) + +hold by construction. ``order=0`` is inverse-distance weighting, which gets +the first identity but not the second — the property this whole module +exists to add. + +The weights themselves are tested through +``underworld3.utilities.rbf_stencil.linear_exact_weights``, because the +partition-of-unity and sparsity claims are about the *weights* and the +value-returning KDTree API cannot expose them. +""" + +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.utilities.rbf_stencil import linear_exact_weights + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +def _stencils(source, target, nnn): + """Brute-force kNN, so the test does not depend on the kd-tree it checks.""" + d = np.linalg.norm(target[:, None, :] - source[None, :, :], axis=2) + idx = np.argsort(d, axis=1)[:, :nnn] + return idx, source[idx] + + +def _linear(coords): + return 0.5 + coords @ np.arange(1, coords.shape[1] + 1, dtype=float) + + +# -------------------------------------------------------------------------- +# The two reproduction identities, at the weights level +# -------------------------------------------------------------------------- +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("nnn_rule", ["minimum", "default", "wide"]) +def test_weights_reproduce_constants_and_linears(dim, nnn_rule): + nnn = {"minimum": dim + 2, "default": 2 * (dim + 1), "wide": 14}[nnn_rule] + rng = np.random.default_rng(20260727 + dim) + source = rng.random((500, dim)) + target = rng.random((80, dim)) + + idx, stencil = _stencils(source, target, nnn) + weights, degenerate = linear_exact_weights(target, stencil) + + assert not degenerate.any(), ( + f"{degenerate.sum()} random-cloud stencils were rejected as degenerate" + ) + + partition = np.abs(weights.sum(axis=1) - 1.0).max() + assert partition < 1.0e-12, f"partition of unity violated by {partition:.3e}" + + reproduced = (weights[:, :, None] * stencil).sum(axis=1) + linear = np.abs(reproduced - target).max() + assert linear < 1.0e-12, f"linear reproduction violated by {linear:.3e}" + + values = _linear(source) + interpolated = (weights * values[idx]).sum(axis=1) + error = np.abs(interpolated - _linear(target)).max() + assert error < 1.0e-12, f"linear field interpolated to {error:.3e}, expected round-off" + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_weights_are_sparse(dim): + """Exactly nnn non-zeros per row: this is the point of a *local* scheme.""" + nnn = 2 * (dim + 1) + rng = np.random.default_rng(99) + source = rng.random((300, dim)) + target = rng.random((40, dim)) + + _, stencil = _stencils(source, target, nnn) + weights, _ = linear_exact_weights(target, stencil) + + assert weights.shape == (40, nnn) + assert np.count_nonzero(weights, axis=1).max() <= nnn + + +# -------------------------------------------------------------------------- +# Degeneracy: the affine block loses rank near boundaries and on graded meshes +# -------------------------------------------------------------------------- +def test_collinear_stencil_is_flagged_not_nan(): + stencil = np.zeros((1, 5, 2)) + stencil[0, :, 0] = np.linspace(0.0, 1.0, 5) + stencil[0, :, 1] = 0.3 + target = np.array([[0.5, 0.7]]) + + weights, degenerate = linear_exact_weights(target, stencil) + + assert degenerate[0], "collinear 2D stencil should be flagged degenerate" + assert np.isfinite(weights).all(), "degenerate stencil must not produce NaN" + + +def test_coplanar_stencil_is_flagged_not_nan(): + stencil = np.zeros((1, 6, 3)) + stencil[0, :, 0] = np.linspace(0.0, 1.0, 6) + stencil[0, :, 1] = np.linspace(0.0, 0.5, 6) + stencil[0, :, 2] = 0.2 + target = np.array([[0.4, 0.2, 0.9]]) + + weights, degenerate = linear_exact_weights(target, stencil) + + assert degenerate[0], "coplanar 3D stencil should be flagged degenerate" + assert np.isfinite(weights).all(), "degenerate stencil must not produce NaN" + + +def test_degenerate_and_healthy_stencils_in_one_batch(): + """A singular member must not take the whole batched solve down with it.""" + healthy = np.array([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]]) + collinear = np.array([[[0.0, 0.0], [0.25, 0.0], [0.5, 0.0], [1.0, 0.0]]]) + stencil = np.concatenate([healthy, collinear]) + target = np.array([[0.3, 0.3], [0.3, 0.3]]) + + weights, degenerate = linear_exact_weights(target, stencil) + + assert not degenerate[0] and degenerate[1] + assert np.isfinite(weights).all() + assert abs(weights[0].sum() - 1.0) < 1.0e-12 + + +def test_stencil_collapsed_onto_target_is_flagged(): + weights, degenerate = linear_exact_weights(np.zeros((1, 2)), np.zeros((1, 4, 2))) + assert degenerate[0] + assert np.isfinite(weights).all() + + +def test_degenerate_stencil_is_retried_on_a_wider_neighbourhood(): + """A locally-coplanar neighbourhood must be escaped, not surrendered to. + + Most source points lie on the plane z = 0, with a sparse off-plane set + further away. A target near the plane has all its nearest neighbours in + it, so the affine block is rank-deficient and the z-gradient of a linear + field is simply not representable from that stencil. Widening reaches the + off-plane points and recovers exactness. + + This matters because a single non-exact node sitting among exact ones is + an isolated spike, which damages derivatives far more than a smooth error + of the same magnitude. Measured on a 3D low-density swarm before the + retry existed: 2.2e-3 max error from 2 nodes in 28824. + """ + # A SMALL coplanar patch -- a local accident, which is what widening is + # for. (A globally coplanar cloud is a different situation, covered by + # test_locally_unrecoverable_cloud_falls_back_loudly.) + axis = np.linspace(0.45, 0.55, 5) + gx, gy = np.meshgrid(axis, axis, indexing="ij") + in_plane = np.stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)], axis=1) + + rng = np.random.default_rng(19) + off_plane = rng.random((80, 3)) + + source = np.vstack([in_plane, off_plane]) + target = np.array([[0.5, 0.5, 0.002], [0.49, 0.51, -0.001]]) + + nnn = 8 + _, stencil = _stencils(source, target, nnn) + _, degenerate = linear_exact_weights(target, stencil) + assert degenerate.all(), ( + "test setup no longer produces coplanar stencils, so it is not " + "exercising the retry path" + ) + + # A z-dependent linear field cannot be recovered from a coplanar stencil, + # so exactness here proves the wider neighbourhood was actually used. + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + + error = np.abs(got[:, 0] - _linear(target)).max() + assert error < 1.0e-12, ( + f"degenerate stencils were not recovered by widening: error {error:.3e}" + ) + + +def test_locally_unrecoverable_cloud_falls_back_loudly(): + """When widening cannot help, say so rather than pretending. + + A source cloud that is coplanar over a wide region is genuinely unable to + determine an affine fit in the third direction. Widening the stencil is + futile, and the honest outcome is a bounded inverse-distance answer plus a + warning -- not a silent loss of the guarantee (the failure mode of #424). + """ + axis = np.linspace(0.0, 1.0, 22) + gx, gy = np.meshgrid(axis, axis, indexing="ij") + source = np.stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)], axis=1) + target = np.array([[0.5, 0.5, 0.02]]) + + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + + with pytest.warns(UserWarning, match="could not support an affine fit"): + got = kdt.rbf_interpolator_local(target, data, 8, 2, False, order=1) + + assert np.isfinite(got).all(), "an unrecoverable stencil must not give NaN" + lo, hi = data.min(), data.max() + assert lo - 1e-12 <= got[0, 0] <= hi + 1e-12, ( + "the inverse-distance fallback should still be bounded" + ) + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_nnn_below_dim_plus_two_is_rejected(dim): + """At dim+1 the affine tail is exactly determined and the scheme reduces to + bare barycentric interpolation on the neighbour simplex — singular whenever + those points are collinear/coplanar, which is common.""" + with pytest.raises(ValueError, match=f"dim \\+ 2 = {dim + 2}"): + linear_exact_weights(np.zeros((1, dim)), np.zeros((1, dim + 1, dim))) + + +# -------------------------------------------------------------------------- +# Through the KDTree API +# -------------------------------------------------------------------------- +@pytest.mark.parametrize("dim", [2, 3]) +def test_kdtree_order1_beats_order0_on_a_linear_field(dim): + nnn = 2 * (dim + 1) + rng = np.random.default_rng(1234 + dim) + source = rng.random((2000, dim)) + target = rng.random((150, dim)) + data = _linear(source)[:, None] + expected = _linear(target) + + kdt = uw.kdtree.KDTree(source) + shepard = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=0) + exact = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + + exact_error = np.abs(exact[:, 0] - expected).max() + shepard_error = np.abs(shepard[:, 0] - expected).max() + + assert exact_error < 1.0e-12, f"order=1 gave {exact_error:.3e} on a linear field" + assert shepard_error > 1.0e-4, ( + "order=0 is expected to smear a linear field; if this fails the " + f"baseline has changed (got {shepard_error:.3e})" + ) + + +def test_kdtree_order0_is_unchanged_by_the_new_arguments(): + """Back-compatibility: the positional call must be bit-identical.""" + rng = np.random.default_rng(5) + source = rng.random((400, 2)) + target = rng.random((50, 2)) + data = _linear(source)[:, None] + + kdt = uw.kdtree.KDTree(source) + positional = kdt.rbf_interpolator_local(target, data, 4, 2, False) + explicit = kdt.rbf_interpolator_local(target, data, 4, 2, False, order=0, + monotone=False) + + assert np.array_equal(positional, explicit) + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_monotone_does_not_touch_a_linear_field(dim): + """The limiter must not cost the guarantee the scheme exists to provide. + + A naive clamp against the stencil's raw min/max does exactly that: a + target outside the convex hull of its own neighbours has a value outside + their range even for an exactly linear field, so such a clamp cannot tell + legitimate extrapolation from ringing. The limiter therefore bounds the + non-affine part only, leaving the linear reconstruction alone. + """ + nnn = 2 * (dim + 1) + rng = np.random.default_rng(77 + dim) + source = rng.random((1500, dim)) + target = rng.random((250, dim)) + data = _linear(source)[:, None] + expected = _linear(target) + + kdt = uw.kdtree.KDTree(source) + unlimited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + limited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1, + monotone=True) + + assert np.abs(limited[:, 0] - expected).max() < 1.0e-12, ( + "the limiter destroyed linear reproduction" + ) + assert np.abs(limited - unlimited).max() < 1.0e-12, ( + "the limiter should be a no-op on a field the scheme reproduces exactly" + ) + + +def test_monotone_bounds_the_correction_on_a_curved_field(): + """It must still do something: bound the RBF correction to the non-affine + variation actually present in the stencil.""" + dim = 2 + nnn = 2 * (dim + 1) + rng = np.random.default_rng(4242) + source = rng.random((800, dim)) + target = rng.random((200, dim)) + values = 0.5 + (source ** 2).sum(axis=1) + np.sin(6.0 * source[:, 0]) + data = values[:, None] + + kdt = uw.kdtree.KDTree(source) + unlimited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + limited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1, + monotone=True) + + assert np.abs(limited - unlimited).max() > 1.0e-6, ( + "the limiter had no effect on a curved field, so it is not limiting" + ) + + # The stated guarantee: the deviation from the local affine trend never + # exceeds the deviation the stencil itself shows. + from underworld3.utilities.rbf_stencil import affine_trend + + _, stencil = _stencils(source, target, nnn) + idx, _ = _stencils(source, target, nnn) + stencil_values = data[idx] + at_target, at_stencil = affine_trend(target, stencil, stencil_values) + residual = stencil_values - at_stencil + correction = limited - at_target + + assert (correction <= residual.max(axis=1) + 1.0e-12).all() + assert (correction >= residual.min(axis=1) - 1.0e-12).all() + + +def test_kdtree_rejects_bad_order_and_monotone_mode(): + rng = np.random.default_rng(11) + source = rng.random((200, 2)) + target = rng.random((10, 2)) + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + + with pytest.raises(ValueError, match="order must be 0"): + kdt.rbf_interpolator_local(target, data, 6, 2, False, order=2) + + with pytest.raises(ValueError, match="no meaning for a local kd-tree"): + kdt.rbf_interpolator_local(target, data, 6, 2, False, order=1, + monotone="pick") + + +# -------------------------------------------------------------------------- +# Oracle tests on a RANDOM field. +# +# Linear reproduction cannot detect wrong-neighbour selection: the identity +# sum_j w_j x_j = x* is evaluated against whatever points were handed to the +# solver, so it holds even if the kd-tree returned the wrong neighbours, and +# a linear field then interpolates exactly anyway. Only a field whose values +# vary independently of position -- i.e. random data -- makes the choice of +# neighbour observable. +# -------------------------------------------------------------------------- +def _reference_interpolate(source, target, values, nnn): + """An independent implementation, written the obvious slow way. + + Brute-force kNN, then one dense saddle solve per target point, with no + local rescaling. Deliberately does not share code with the library. + """ + out = np.empty(target.shape[0]) + dim = source.shape[1] + for i, x in enumerate(target): + d = np.linalg.norm(source - x, axis=1) + idx = np.argsort(d)[:nnn] + pts = source[idx] + + def phi(r): + r = np.where(r == 0.0, 1.0e-300, r) + return r ** 2 * np.log(r) + + rr = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=2) + P = np.hstack([np.ones((nnn, 1)), pts]) + M = np.block([[phi(rr), P], [P.T, np.zeros((dim + 1, dim + 1))]]) + rhs = np.concatenate([phi(np.linalg.norm(pts - x, axis=1)), [1.0], x]) + w = np.linalg.solve(M, rhs)[:nnn] + out[i] = w @ values[idx] + return out + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_matches_an_independent_implementation_on_random_data(dim): + """Random values, so wrong-neighbour selection cannot hide.""" + nnn = 2 * (dim + 1) + rng = np.random.default_rng(31337 + dim) + source = rng.random((400, dim)) + target = 0.15 + 0.7 * rng.random((40, dim)) + values = rng.standard_normal(source.shape[0]) + + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, values[:, None], nnn, 2, False, order=1) + want = _reference_interpolate(source, target, values, nnn) + + error = np.abs(got[:, 0] - want).max() + assert error < 1.0e-8, ( + f"disagrees with an independent implementation by {error:.3e} on " + "random data — check neighbour selection, not just reproduction" + ) + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_kdtree_neighbours_match_brute_force_on_random_data(dim): + """The value must change if the stencil is wrong, and it does not.""" + nnn = 2 * (dim + 1) + rng = np.random.default_rng(555 + dim) + source = rng.random((600, dim)) + target = 0.2 + 0.6 * rng.random((50, dim)) + values = rng.standard_normal(source.shape[0]) + + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, values[:, None], nnn, 2, False, order=1) + + idx, stencil = _stencils(source, target, nnn) + weights, degenerate = linear_exact_weights(target, stencil) + assert not degenerate.any() + want = (weights * values[idx]).sum(axis=1) + + assert np.abs(got[:, 0] - want).max() < 1.0e-10 + + # Control: a deliberately wrong stencil must give a different answer, or + # the test above proves nothing. + rolled = np.roll(idx, 1, axis=0) + wrong_weights, _ = linear_exact_weights(target, source[rolled]) + wrong = (wrong_weights * values[rolled]).sum(axis=1) + assert np.abs(wrong - want).max() > 1.0e-3, ( + "shuffling the stencils did not change the result, so this test " + "cannot detect wrong-neighbour selection" + ) + + +# -------------------------------------------------------------------------- +# The operator form +# -------------------------------------------------------------------------- +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("order", [0, 1]) +def test_interpolation_matrix_agrees_with_the_value_path(dim, order): + """`T @ data` must be what the value API returns, or they will drift.""" + rng = np.random.default_rng(808 + dim) + source = rng.random((500, dim)) + target = 0.1 + 0.8 * rng.random((60, dim)) + values = rng.standard_normal((source.shape[0], 2)) + + kdt = uw.kdtree.KDTree(source) + T = kdt.interpolation_matrix(target, order=order) + direct = kdt.rbf_interpolator_local(target, values, None, 2, False, order=order) + + assert T.shape == (target.shape[0], source.shape[0]) + assert np.abs(T @ values - direct).max() < 1.0e-12 + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_interpolation_matrix_rows_are_never_empty(dim): + """The raw-weights helper zeroes degenerate rows; the operator must not. + + A zero row would silently interpolate to zero, which is worse than the + inverse-distance fallback it replaces. + """ + axis = np.linspace(0.0, 1.0, 12) + grids = np.meshgrid(*([axis] * dim), indexing="ij") + source = np.stack([g.ravel() for g in grids], axis=1) + rng = np.random.default_rng(2) + target = rng.random((80, dim)) + + kdt = uw.kdtree.KDTree(source) + T = kdt.interpolation_matrix(target, order=1) + + per_row = np.diff(T.indptr) + assert per_row.min() > 0, "empty row in the interpolation operator" + row_sums = np.asarray(T.sum(axis=1)).ravel() + assert np.abs(row_sums - 1.0).max() < 1.0e-12 + + +# -------------------------------------------------------------------------- +# Contract: determinism and purity +# -------------------------------------------------------------------------- +def test_repeated_calls_are_bit_identical(): + """No global state, no RNG, no accumulation: the weights are a pure + function of the geometry.""" + rng = np.random.default_rng(64) + source = rng.random((300, 3)) + target = rng.random((40, 3)) + values = rng.standard_normal((300, 1)) + + kdt = uw.kdtree.KDTree(source) + first = kdt.rbf_interpolator_local(target, values, 8, 2, False, order=1) + second = kdt.rbf_interpolator_local(target, values, 8, 2, False, order=1) + assert np.array_equal(first, second) + + other = uw.kdtree.KDTree(source.copy()) + assert np.array_equal( + other.rbf_interpolator_local(target, values, 8, 2, False, order=1), first + ) + + +def test_equidistant_neighbours_resolve_deterministically(): + """Ties in the kNN search must not make the result run-dependent. + + Which of several equidistant source points is chosen is not specified — + but it must be the same choice every time, or results stop reproducing. + """ + # A target at the centre of a symmetric ring: all ring points are exactly + # equidistant, so the nnn selection is a tie. + angles = np.arange(12) * (2.0 * np.pi / 12.0) + ring = np.stack([np.cos(angles), np.sin(angles)], axis=1) + source = np.vstack([ring, 2.0 * ring]) + target = np.array([[0.0, 0.0]]) + values = np.arange(source.shape[0], dtype=float)[:, None] + + kdt = uw.kdtree.KDTree(source) + results = [ + kdt.rbf_interpolator_local(target, values, 6, 2, False, order=1) + for _ in range(5) + ] + for r in results[1:]: + assert np.array_equal(r, results[0]), "tied kNN selection is not deterministic" + + fresh = uw.kdtree.KDTree(source.copy()) + assert np.array_equal( + fresh.rbf_interpolator_local(target, values, 6, 2, False, order=1), results[0] + ) + + +def test_nearest_neighbour_path_keeps_its_shape(): + """`nnn=1` returns source rows directly, with no stencil axis. + + This is the checkpoint round-trip path (`read_timestep` uses `nnn=1`). + A refactor that swapped `query(k=1)` for `find_closest_n_points` silently + added an axis here, and only a snapshot test noticed. + """ + rng = np.random.default_rng(12) + source = rng.random((200, 2)) + target = rng.random((30, 2)) + data = rng.standard_normal((200, 1)) + + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, data, 1, 2, False) + + assert got.shape == (30, 1), f"nnn=1 returned shape {got.shape}, expected (30, 1)" + + # It really is the nearest neighbour's value, not an average. + d = np.linalg.norm(target[:, None, :] - source[None, :, :], axis=2) + assert np.array_equal(got[:, 0], data[np.argmin(d, axis=1), 0]) + + +@pytest.mark.parametrize("order", [0, 1]) +def test_default_nnn_works_in_both_dimensions(order): + """A default that raises is not a default (the fixed nnn=4 did, at + order=1 in 3D, where dim + 2 = 5).""" + for dim in (2, 3): + rng = np.random.default_rng(9) + source = rng.random((300, dim)) + target = rng.random((20, dim)) + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, data, order=order) + assert got.shape == (20, 1) + assert np.isfinite(got).all() + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_error_falls_with_point_spacing_faster_than_inverse_distance(dim): + """On a quadratic field, order=1 should converge under refinement and + order=0 should barely move — inverse-distance error is set by the stencil + geometry, not by how close the points are.""" + + def quadratic(coords): + return 0.5 + (coords ** 2).sum(axis=1) + + nnn = 2 * (dim + 1) + errors = {} + for n_per_side in (10, 20): + axis = (np.arange(n_per_side) + 0.5) / n_per_side + grid = np.meshgrid(*([axis] * dim), indexing="ij") + source = np.stack([g.ravel() for g in grid], axis=1) + rng = np.random.default_rng(3) + target = 0.2 + 0.6 * rng.random((100, dim)) + data = quadratic(source)[:, None] + expected = quadratic(target) + + kdt = uw.kdtree.KDTree(source) + for order in (0, 1): + got = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=order) + errors[(order, n_per_side)] = np.abs(got[:, 0] - expected).max() + + order1_ratio = errors[(1, 10)] / errors[(1, 20)] + order0_ratio = errors[(0, 10)] / errors[(0, 20)] + + assert order1_ratio > 3.0, ( + f"order=1 should converge at least ~O(h^2) under halving the spacing, " + f"got a factor of {order1_ratio:.2f}" + ) + assert order1_ratio > order0_ratio, ( + f"order=1 convergence ({order1_ratio:.2f}x) should beat order=0 " + f"({order0_ratio:.2f}x)" + ) diff --git a/tests/test_0506_swarm_proxy_linear_exact.py b/tests/test_0506_swarm_proxy_linear_exact.py new file mode 100644 index 00000000..b2ad40bb --- /dev/null +++ b/tests/test_0506_swarm_proxy_linear_exact.py @@ -0,0 +1,120 @@ +"""Swarm proxy variables transfer a linear field without smearing it. + +A proxy mesh variable is refreshed from the particles by local RBF +interpolation. That transfer used to be inverse-distance (Shepard) weighting, +which reproduces constants but not gradients; the proxy default is now +``order=1``, which reproduces both. + +Why a *linear* field is the right probe: it lies exactly inside both the P1 +and the P2 proxy space, so the finite element discretisation contributes +nothing at all to the measured error. Everything left is particle -> node +transfer error, which makes the assertion below unambiguous. + +Measured on this configuration when the default was changed (2026-07-27): +inverse distance gave a relative max error of ~5e-3 on the linear field and +~9e-3 on a quadratic one; the linear-exact scheme gives round-off and ~1e-4 +respectively. +""" + +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.meshing import UnstructuredSimplexBox + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +def _linear(coords): + return 0.5 + coords @ np.arange(1, coords.shape[1] + 1, dtype=float) + + +def _quadratic(coords): + return 0.5 + (coords ** 2).sum(axis=1) + + +@pytest.fixture +def proxied_swarm(): + mesh = UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8.0 + ) + swarm = uw.swarm.Swarm(mesh) + var = swarm.add_variable(name="f", size=1, proxy_degree=1) + swarm.populate(fill_param=4) + yield swarm, var + del swarm + del mesh + + +def _proxy_error(var, field, **refresh): + """Refresh the proxy from an analytic field and return its relative error.""" + proxy = var._meshVar + particle_coords = var.swarm._particle_coordinates.data.copy() + node_coords = np.asarray(proxy.coords) + + var.data[:, 0] = field(particle_coords) + var._rbf_to_meshVar(proxy, **refresh) + + expected = field(node_coords) + return np.abs(np.asarray(proxy.data[:, 0]) - expected).max() / np.abs(expected).max() + + +def test_proxy_default_reproduces_a_linear_field(proxied_swarm): + """The shipped default must be linear-exact — no arguments passed.""" + _, var = proxied_swarm + error = _proxy_error(var, _linear) + assert error < 1.0e-12, ( + f"proxy of an exactly linear field carries a relative error of " + f"{error:.3e}; the default transfer is not linear-exact" + ) + + +def test_inverse_distance_smears_the_same_linear_field(proxied_swarm): + """The control arm: without order=1 the error is real and much larger. + + This is what stops the test above passing for the wrong reason (e.g. a + tolerance that any scheme would meet on this mesh). + """ + _, var = proxied_swarm + exact = _proxy_error(var, _linear, order=1) + shepard = _proxy_error(var, _linear, order=0, nnn=3) + + assert shepard > 1.0e-4, ( + f"inverse distance is expected to smear a linear field, got {shepard:.3e}" + ) + assert shepard > 1.0e6 * max(exact, 1.0e-16) + + +def test_proxy_default_improves_a_quadratic_field(proxied_swarm): + """Linear exactness is not a trick that only helps linear fields.""" + _, var = proxied_swarm + exact = _proxy_error(var, _quadratic, order=1) + shepard = _proxy_error(var, _quadratic, order=0, nnn=3) + + assert exact < shepard / 10.0, ( + f"order=1 gave {exact:.3e} vs inverse distance {shepard:.3e} on a " + "quadratic field; expected at least a 10x improvement" + ) + + +def test_proxy_still_reproduces_a_constant_exactly(proxied_swarm): + """Constants were already exact under inverse distance — do not regress.""" + _, var = proxied_swarm + error = _proxy_error(var, lambda coords: np.full(coords.shape[0], 3.25)) + assert error < 1.0e-12, f"constant field carries error {error:.3e}" + + +def test_proxy_monotone_keeps_the_linear_field_exact(proxied_swarm): + """Turning the limiter on must not cost the proxy its exactness. + + Proxy nodes routinely sit outside the convex hull of their own particle + stencil — not only at the domain boundary — so a limiter that clipped + against the stencil's raw min/max would fire on the linear part. This one + limits only the non-affine correction. + """ + _, var = proxied_swarm + limited = _proxy_error(var, _linear, monotone=True) + assert limited < 1.0e-12, ( + f"proxy with the limiter on carries {limited:.3e} on an exactly linear " + "field; the limiter is clipping the linear reconstruction" + )