From 4fc1e2be99060ac6d2729fb39039cb85f2e78a32 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 09:51:42 +1000 Subject: [PATCH 1/9] A point the local mesh cannot own should be rejected, not walked 50 times The lost-point walk in _get_closest_local_cells_internal tried the 50 nearest cell centroids and only stopped when every lost point had been found. A point no local cell could own therefore paid all 50 rounds - 51 containment tests against 1 for an owned point, measured - and every point already located was re-tested on each of them, so one unfindable point charged the whole batch. In parallel the fraction of points a rank does not own is exactly what grows with rank count. Two bounds. A point inside a cell is no further from its nearest kd-tree control point than from that cell's centroid, which for a convex cell is within the cell's vertex reach; the largest local reach is now recorded with the kd-tree and a lost point beyond twice that distance is rejected before the walk starts. The factor of two is slack for the in-cell test's face tolerance and for badly shaped cells. And the working set shrinks: a point leaves as soon as a cell claims it, or as soon as the sorted neighbour distances pass the rejection radius. Measured on the #551 dossier probes: a foreign point costs 1.0 containment tests instead of 51.0 (0.53 us/point instead of 7.99 serial, 0.79 instead of 9.6 at np=4), and one unfindable point in a batch of 1000 adds 3 point-tests instead of 991 (2-D) or 12251 (3-D). The nearest containing centroid now wins, where before the winner was the LAST of up to 50 rounds and so depended on whether some unrelated point in the same batch was findable. Across a nine-set battery in 2-D and 3-D at np=1/2/4, 15 of ~450000 located cells change; every one is a point that both the old and the new cell contain. Addresses item 1 of #551. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 80 +++++++++++++++++-- 1 file changed, 74 insertions(+), 6 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index f141d859..dd42d5cb 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -5288,6 +5288,12 @@ def _build_kd_tree_index(self): control_points_list = [] control_points_cell_list = [] centroids_list = [] + # Largest distance from a cell centroid to one of that cell's own + # vertices, maximised over local cells. A convex cell is the convex + # hull of its vertices, so every point of it lies within this distance + # of its centroid — which makes it the rejection radius the locator + # needs (see _get_closest_local_cells_internal). + cell_reach = 0.0 for cell, cell_id in enumerate(range(cStart, cEnd)): @@ -5297,6 +5303,8 @@ def _build_kd_tree_index(self): cell_point_coords = nav_coords[self._coord_rows_for_points(nav_dm, points)] cell_centroid = cell_point_coords.mean(axis=0) centroids_list.append(cell_centroid) + cell_reach = max(cell_reach, float(numpy.linalg.norm( + cell_point_coords - cell_centroid, axis=1).max())) # for face in range(cell_num_faces): @@ -5360,6 +5368,10 @@ def _build_kd_tree_index(self): centroids_list, dtype=numpy.float64).reshape(-1, self.cdim) self._centroid_index = uw.kdtree.KDTree(self._nav_centroids) + # Rejection radius for the lost-point walk. Rebuilt with the kd-tree + # (i.e. invalidated by deform / adapt along with _index). + self._local_cell_reach = cell_reach + return def _build_kd_tree_index_PIC(self): @@ -5984,6 +5996,10 @@ def get_closest_cells(self, coords: numpy.ndarray) -> numpy.ndarray: # CRITICAL: Must return 1D array, not 2D, for Cython buffer compatibility return numpy.array([], dtype=numpy.int64) + # Safety factor on the local cell reach used to reject a query point + # before the lost-point walk. See _get_closest_local_cells_internal. + _LOCATOR_REACH_MARGIN = 2.0 + def _get_closest_local_cells_internal( self, coords: numpy.ndarray, @@ -5997,6 +6013,12 @@ def _get_closest_local_cells_internal( is not guaranteed. Also compares the distance from the cell to the point - if this is larger than the "cell size" then returns -1 + A point the first containment test rejects is looked for among the + nearest cell centroids. Points too far from the local mesh to be in + any of its cells are rejected before that walk starts, and a point + leaves the walk as soon as a cell claims it, so the walk costs what is + still lost rather than what was asked for. + ``on_boundary`` and ``tol`` are forwarded to the in-cell containment test (see ``_test_if_points_in_cells_internal``). Default ``(on_boundary=True, tol=0.0)`` admits on-face queries @@ -6046,7 +6068,8 @@ def _get_closest_local_cells_internal( self._build_kd_tree_index() if len(coords) > 0: - dist, closest_points = self._index.query(coords, k=1, sqr_dists=False) + control_point_distance, closest_points = self._index.query( + coords, k=1, sqr_dists=False) # >= : valid indices are 0..n-1, and the empty-tree sentinel # (0 with n=0) must trip this guard, not index _indexMap (#399). if np.any(closest_points >= self._index.n): @@ -6072,7 +6095,32 @@ def _get_closest_local_cells_internal( cells[~inside] = -1 lost_points = np.where(inside == False)[0] - # Part 2 - try to find the lost points by walking nearby cells + if lost_points.shape[0] == 0: + return cells + + # Part 2 - try to find the lost points by walking nearby cells. + # + # Reject what cannot possibly be found, first. Every cell contributes + # its centroid to the control-point kd-tree, so a point lying in cell c + # is at most |p - centroid_c| from its NEAREST control point, and a + # convex cell puts that within the cell's vertex reach. A lost point + # whose nearest control point is beyond the largest local reach is in + # no local cell and the walk has nothing to find for it. Without this + # every genuinely foreign point pays the full 50-neighbour walk — 51 + # containment tests against 1 for an owned point — and the foreign + # fraction is exactly what grows with rank count (#551). + # + # The margin is deliberately loose. The in-cell test admits a thin + # slab outside each face (``tol``), and a badly shaped cell expands + # further under that slab than a well-shaped one; a factor of two on + # the reach covers both with room to spare while still rejecting + # anything more than about one cell away from the local mesh. + reach = getattr(self, "_local_cell_reach", None) + if reach is not None and reach > 0.0: + reach = self._LOCATOR_REACH_MARGIN * reach + lost_points = lost_points[control_point_distance[lost_points] <= reach] + if lost_points.shape[0] == 0: + return cells # Size by the nav-DM cell count, which is what _centroid_index # was built from (includes ghost cells on manifold meshes). @@ -6082,22 +6130,42 @@ def _get_closest_local_cells_internal( num_local_cells = nav_centroids.shape[0] num_testable_neighbours = min(num_local_cells, 50) - dist2, closest_centroids = self._centroid_index.query( + centroid_distance, closest_centroids = self._centroid_index.query( coords[lost_points], k=num_testable_neighbours, sqr_dists=False ) + # The kd-tree drops the neighbour axis at k == 1 (a rank owning a + # single cell); the walk indexes it either way. + centroid_distance = centroid_distance.reshape(lost_points.shape[0], -1) + closest_centroids = closest_centroids.reshape(lost_points.shape[0], -1) # This number is close to the point-point coordination value in 3D unstructured # grids (by inspection) + # The working set shrinks: a point drops out as soon as a neighbour + # claims it, or as soon as the neighbour distances (sorted, so + # monotonic in i) pass the rejection radius. The nearest containing + # centroid therefore wins, which also makes the answer independent of + # whether some OTHER point in the same batch is findable — previously + # a single unlocatable point kept every already-found point in the + # test set for all 50 rounds, and a shared-face point could be + # reassigned to a further cell in a later round. + working = np.arange(lost_points.shape[0]) for i in range(0, num_testable_neighbours): + if reach is not None and reach > 0.0: + working = working[centroid_distance[working, i] <= reach] + if working.shape[0] == 0: + break + + candidate_cells = closest_centroids[working, i] inside = self._test_if_points_in_cells_internal( - coords[lost_points], closest_centroids[:, i], + coords[lost_points[working]], candidate_cells, on_boundary=on_boundary, tol=tol, ) - cells[lost_points[inside]] = closest_centroids[inside, i] + cells[lost_points[working[inside]]] = candidate_cells[inside] - if np.count_nonzero(cells == -1) == 0: + working = working[~inside] + if working.shape[0] == 0: break return cells From 61de052325727d939bcc11a610b56490b500f13f Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 09:52:02 +1000 Subject: [PATCH 2/9] The cell hint that bypasses DMLocatePoints must contain the point Serial simplex meshes assert "exact" location capability, which makes the UW3 cell hint authoritative and skips PETSc's DMLocatePoints entirely - while the hint they handed over was get_closest_cells, a nearest-CONTROL-POINT kd-tree lookup with no containment test at all. On a tetrahedron nothing downstream can rescue that: the only remaining guard is a componentwise box clamp on the reference coordinates, and the reference tet is not the reference box. A query on a shared edge was answered by extrapolating the basis of a cell that does not contain it. That is #432, a recurrence of #390. The serial-simplex branch now takes _robust_owning_cells, the containment- checked locator every other authoritative path already used, so the three branches collapse to one call. Points it cannot place come back as -1, which the C bypass treats as a non-claim; they surface in unlocated_mask and take the RBF fallback that is already plumbed. The alternatives were worse. Restricting the "exact" assertion pushes serial simplex meshes back onto DMLocatePoints, which is slower and re-opens the #390 class of silent drops the bypass was added to close. A barycentric clamp in C only pins a wrong cell's reference coordinates to that cell's boundary, which gives the edge value of the WRONG cell - right only for continuous fields, and no help at all when the nominated cell is not adjacent to the point. Measured cost: one extra containment test per point, 1.0 tests per point for interior, on-face and on-edge queries. Measured benefit: 3-D P1 evaluation at quarter-points along cell edges was off by 3.0e-01 at 17 of 1318 points and is now exact. Addresses item 3 of #551 and closes #432. Underworld development team with AI support from Claude Code --- src/underworld3/function/_function.pyx | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index ee13cb97..32d2b0a9 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -1311,24 +1311,22 @@ def petsc_interpolate( expr, cached_info = CachedDMInterpolationInfo() # Cell hints, by policy: - # AUTHORITATIVE — the estimator owner is the answer. Parallel uses - # the bulletproof barycentric locator (correct owner across - # seams). Serial simplex keeps get_closest_cells (the validated - # bit-for-bit PR #203 path; with planar faces + ξ-clamp the - # nearest-cell hint evaluates exactly). Serial quad/hex meshes - # that qualify by MEASUREMENT use the estimator directly — the - # nearest-centroid guess is not containment-checked and these - # meshes only just graduated, so take the checked owner. + # AUTHORITATIVE — the hint bypasses DMLocatePoints, so it has to + # be a cell that CONTAINS the point. _robust_owning_cells is the + # containment-checked locator: it returns a cell whose walls the + # point is inside (any one of them, for a point on a shared face) + # and -1 when no local cell contains it. Every authoritative mesh + # takes the same route. Serial simplex meshes used to take the + # nearest-CONTROL-POINT lookup (get_closest_cells) with no + # containment test at all; on a tetrahedron the reference-coord + # clamp downstream is a box clamp and cannot rescue that, so a + # query on a shared edge was evaluated by extrapolating the basis + # of a cell that does not contain it (#432, a recurrence of #390). # NOT AUTHORITATIVE — no hint at all: DMLocatePoints decides, # dropped points surface in unlocated_mask and are filled by the # RBF fallback below. if authoritative: - if mesh._eval_use_robust_location(): - cells = mesh._robust_owning_cells(coords) - elif not bool(mesh.dm.isSimplex()) and mesh.dim == mesh.cdim: - cells = mesh._robust_owning_cells(coords) - else: - cells = mesh.get_closest_cells(coords) + cells = mesh._robust_owning_cells(coords) else: cells = None From 8c1faf4219d72df673f86562cb3ba6b974af8ff8 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 09:52:20 +1000 Subject: [PATCH 3/9] Classify and locate in one pass, so evaluation does not search twice points_in_domain located the points near the domain boundary, returned a boolean mask and threw the owning cells away; petsc_interpolate then located every interior point again. The near-boundary points were being searched for twice per evaluate call. Mesh._classify_points_in_domain returns both - the mask and the cells the classification actually looked up, with -1 meaning "not looked up" for an interior point and "not in the local mesh" for an exterior one. evaluate passes those to petsc_interpolate as cell_hints, which searches only for the entries still marked -1, and only on the DMInterpolation cache miss that needs them. A cache hit locates nothing at all, which the first version of this change got wrong: filling the whole hint array in the classifier made serial evaluate 29% slower because it searched on every call, cached or not. points_in_domain keeps its signature, its answer and its cost - it is now a one-line wrapper and does not search on the interpolator's behalf. Only the robust locator's answer is kept as a hint; the cell-wall test the serial classifier uses runs at a different face tolerance and its answer is a classification, not a hint. Measured on the #551 dossier probes: containment point-tests per located point inside evaluate fall from 2.28 to 1.10 at np=2, 3.35 to 1.20 at np=4 and 4.19 to 1.27 at np=8; evaluate wall time 0.0382 to 0.0361 s at np=4. Addresses item 2 of #551. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 56 +++++++++++++++++-- src/underworld3/function/_function.pyx | 43 ++++++++++++-- 2 files changed, 90 insertions(+), 9 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index dd42d5cb..3bf5747c 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -5879,6 +5879,38 @@ def points_in_domain(self, points, strict_validation=True): Whether to perform strict validation near boundaries """ + return self._classify_points_in_domain(points, strict_validation)[0] + + def _classify_points_in_domain(self, points, strict_validation=True): + """In/out classification, keeping the owning cells it had to locate. + + ``points_in_domain`` located the near-boundary points and threw the + owning cells away, leaving the interpolator to locate them a second + time. This hands them over instead, so evaluation locates each point + once (#551 item 2). + + Parameters + ---------- + points : array-like + Coordinate array in any physical unit system (will be + auto-converted). + strict_validation : bool + Whether to perform strict validation near boundaries. + + Returns + ------- + in_or_not : numpy.ndarray of bool + Exactly what :meth:`points_in_domain` returns. + cells : numpy.ndarray of int + Owning cell for the points the classification actually located, + at the evaluation face tolerance (:meth:`_robust_owning_cells`). + ``-1`` everywhere else — for an EXTERIOR point that means "not in + the local mesh"; for an INTERIOR point it means "not looked up", + because the boundary-sign test settled it without a search. A + caller that needs a cell for every interior point locates the + ``-1`` entries itself, and only when it needs them: nothing here + searches on the classifier's behalf. + """ # Convert points to model coordinates using the unified conversion function # This handles all coordinate formats: plain numbers, unit-aware coordinates, lists, tuples, arrays import underworld3 as uw @@ -5893,7 +5925,10 @@ def points_in_domain(self, points, strict_validation=True): max_radius = self.get_max_radius() if model_points.shape[0] == 0: - return numpy.array([], dtype=bool) + return (numpy.array([], dtype=bool), + numpy.array([], dtype=numpy.int64)) + + cells = numpy.full(model_points.shape[0], -1, dtype=numpy.int64) # Cd-1 surface mesh: no boundary-face control points exist # (see _mark_local_boundary_faces_inside_and_out). Per the @@ -5901,7 +5936,8 @@ def points_in_domain(self, points, strict_validation=True): # the manifold; the closest-local-cell test is the right # filter, not an inside/outside split. if self.boundary_face_control_points_kdtree is None: - return self._get_closest_local_cells_internal(model_points) != -1 + in_or_not = self._get_closest_local_cells_internal(model_points) != -1 + return in_or_not, cells dist2, closest_control_points_ext = self.boundary_face_control_points_kdtree.query( model_points, k=1, sqr_dists=True @@ -5929,11 +5965,17 @@ def points_in_domain(self, points, strict_validation=True): # cell (>= 0) for any point genuinely in/on the mesh and -1 only for # true exterior. Serial / non-simplex keep the cell-wall test # (bit-identical to the validated baseline). + # + # Only the robust locator's answer is kept as a cell hint: it is the + # same call the evaluation path makes, so keeping it saves a repeat. + # The cell-wall test runs at a different face tolerance and its answer + # is a classification, not a hint. near_boundary = numpy.where(dist2 < 2 * max_radius**2)[0] near_boundary_points = model_points[near_boundary] if self._eval_use_robust_location(): - in_or_not[near_boundary] = self._robust_owning_cells(near_boundary_points) >= 0 + cells[near_boundary] = self._robust_owning_cells(near_boundary_points) + in_or_not[near_boundary] = cells[near_boundary] >= 0 else: in_or_not[near_boundary] = ( self._get_closest_local_cells_internal(near_boundary_points) != -1 @@ -5943,11 +5985,15 @@ def points_in_domain(self, points, strict_validation=True): chosen_ones = numpy.where(in_or_not == True)[0] chosen_points = model_points[chosen_ones] if self._eval_use_robust_location(): - in_or_not[chosen_ones] = self._robust_owning_cells(chosen_points) >= 0 + cells[chosen_ones] = self._robust_owning_cells(chosen_points) + in_or_not[chosen_ones] = cells[chosen_ones] >= 0 else: in_or_not[chosen_ones] = self._get_closest_local_cells_internal(chosen_points) != -1 - return in_or_not + # A point demoted to exterior keeps no hint: it goes to RBF. + cells[~in_or_not] = -1 + + return in_or_not, cells @timing.routine_timer_decorator def get_closest_cells(self, coords: numpy.ndarray) -> numpy.ndarray: diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 32d2b0a9..7be624dd 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -1023,13 +1023,19 @@ def evaluate_nd( expr, # same fix that lets swarm migration claim them. Serial / non-simplex # keep the cell-wall test (bit-identical). See # parallel-repeated-solve-corruption.md. - in_or_not = mesh.points_in_domain(coords_array, strict_validation=False) + # + # The classification also hands back the cells it located on the way, + # so petsc_interpolate does not look those points up again (#551 + # item 2). + in_or_not, cell_hints = mesh._classify_points_in_domain( + coords_array, strict_validation=False) evaluation_interior = petsc_interpolate( expr, coords_array[in_or_not], coord_sys, mesh, simplify=simplify, - verbose=verbose, ) + verbose=verbose, + cell_hints=cell_hints[in_or_not], ) evaluation_interior = np.atleast_1d(evaluation_interior) # handle case where there is only 1 interior point @@ -1093,7 +1099,8 @@ def petsc_interpolate( expr, mesh=None, other_arguments=None, simplify=True, - verbose=False, ): + verbose=False, + cell_hints=None, ): """ Evaluate a given expression at a list of coordinates. @@ -1112,6 +1119,14 @@ def petsc_interpolate( expr, other_arguments: dict Dictionary of other arguments necessary to evaluate function. Not yet implemented. + cell_hints: numpy.ndarray, optional + One owning cell index per coordinate, as returned by + ``Mesh._classify_points_in_domain``: a cell the classification + already located, or ``-1`` for "not looked up", which this function + then locates itself. Supplying it means those points are not located + twice. Hints must have been located against ``mesh``; hints for any + other mesh in the expression are ignored and that mesh locates its + own. Notes ----- @@ -1218,6 +1233,10 @@ def petsc_interpolate( expr, # 2. Evaluate all mesh variables - there is no real # computational benefit in interpolating a subset. + # Any cell hints the caller supplied were located against THIS mesh; an + # expression spanning two meshes must locate the second one itself. + hinted_mesh = mesh + def interpolate_vars_on_mesh( varfns, np.ndarray coords ): """ This function performs the interpolation for the given variables @@ -1325,8 +1344,24 @@ def petsc_interpolate( expr, # NOT AUTHORITATIVE — no hint at all: DMLocatePoints decides, # dropped points surface in unlocated_mask and are filled by the # RBF fallback below. + # + # Cells the caller's classification already located are reused; + # only the ones it left at -1 are searched for, and only here, on + # the cache miss that actually needs them. That is what makes it + # one location per point per call rather than two. if authoritative: - cells = mesh._robust_owning_cells(coords) + if cell_hints is not None and mesh is hinted_mesh: + cells = np.ascontiguousarray(cell_hints, dtype=np.int64) + if cells.shape[0] != coords.shape[0]: + raise RuntimeError( + "cell_hints must carry one cell index per coordinate " + f"({cells.shape[0]} hints for {coords.shape[0]} points)." + ) + unhinted = np.where(cells < 0)[0] + if unhinted.shape[0] > 0: + cells[unhinted] = mesh._robust_owning_cells(coords[unhinted]) + else: + cells = mesh._robust_owning_cells(coords) else: cells = None From 82145992e20b236906f08aca07fd980de03a5aeb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 09:52:39 +1000 Subject: [PATCH 4/9] Pin the locator's two contracts: the hint contains the point, rejection is O(1) #432 and #390 have both been fixed and returned, and there was no regression test pinning on-edge or on-face queries. tests/test_0761_point_locator.py adds them in 2-D and 3-D. The oracle is the closed-form P1 value (1-t)*u_a + t*u_b along a cell edge, computed WITHOUT uw.function.evaluate so it cannot inherit the defect it is testing for - the same discipline test_0753's arbitrary-coarse-field reference uses, and for the same reason. Covers shared vertices, edge midpoints, edge quarter-points and 3-D face centroids. The performance guard counts containment tests per point rather than seconds. The count is what the algorithm does; a wall time is what the machine was doing at the time, and the house has been burned by timing tests before. Parallel wall-clock numbers belong in the PR, not in an assert. Negative controls, because a test that cannot fail proves nothing. The nodal field is MEASURED to vary by more than 0.1 between a cell edge and its midpoint, so a wrong-cell interpolant cannot pass by being smooth (a linear field cannot test neighbour selection). get_closest_cells is shown to nominate cells that do not contain the query, so the edge test is pinning something real - and only at 3-D quarter-points, since a midpoint is never misassigned and 2-D never is, which is why #432 is a 3-D report; the 2-D case skips with that count in its message. The containment counter is asserted to fire, so the bounds cannot pass by instrumenting nothing. Validated fail-before with the fixes stashed and the tree rebuilt: 3-D P1 evaluation at t=0.25 off by 3.036e-01 at 17 of 1318 points; a point outside the mesh costing 51.0 containment tests; one unlocatable point adding 991 (2-D) and 12251 (3-D) point-tests to a batch of 1000. Underworld development team with AI support from Claude Code --- tests/test_0761_point_locator.py | 363 +++++++++++++++++++++++++++++++ 1 file changed, 363 insertions(+) create mode 100644 tests/test_0761_point_locator.py diff --git a/tests/test_0761_point_locator.py b/tests/test_0761_point_locator.py new file mode 100644 index 00000000..9a27f909 --- /dev/null +++ b/tests/test_0761_point_locator.py @@ -0,0 +1,363 @@ +"""Point location: the cell hint must contain the point, and a point no local +cell can own must be rejected in O(1) (#551 items 1 and 3; #432, a recurrence +of #390). + +Two properties are pinned here, both of which have been broken before. + +**The hint contains the point.** ``petsc_interpolate`` may bypass PETSc's +``DMLocatePoints`` and evaluate the basis directly in a cell UW3 nominates. +Serial simplex meshes used to nominate the cell owning the nearest kd-tree +CONTROL POINT, with no containment test at all. On a tetrahedron nothing +downstream can rescue that: the reference-coordinate guard is a componentwise +box clamp and the reference tet is not the reference box, so a query on a +shared edge was answered by extrapolating the basis of a cell that does not +contain it. The oracle here is the closed-form P1 value ``(1-t)*u_a + t*u_b`` +along an edge — the same reference +``test_0753_nested_mg_prolongation.py::test_reproduces_an_arbitrary_coarse_field`` +uses, and for the same reason: it is computed WITHOUT ``uw.function.evaluate``, +so it cannot inherit the defect it is testing for. + +**Rejection is cheap.** A point outside the local mesh used to enter a +50-nearest-centroid walk whose only exit was every lost point being found — +51 containment tests against 1 for an owned point. The guard here counts +containment tests rather than seconds: an instrumented counter says what the +algorithm does, where a stopwatch says what the machine was doing at the time. +The parallel wall-clock table lives in the PR, not in a brittle timing assert. +""" +import numpy as np +import pytest +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +# High enough in frequency that neighbouring nodal values are uncorrelated: a +# smooth low-frequency field would be reproduced by an interpolant anchored to +# the WRONG cell almost as well as by the right one, and the test would pass +# while the locator was broken (the "linear field cannot test neighbour +# selection" trap). test_nodal_signal_discriminates_between_cells measures +# this rather than assuming it. +def _nodal_signal(coords): + coords = np.asarray(coords) + signal = np.sin(97.0 * coords[:, 0]) * np.cos(89.0 * coords[:, 1]) + if coords.shape[1] == 3: + signal = signal * np.sin(83.0 * coords[:, 2]) + return signal + + +def _box(dim, cell_size): + if dim == 2: + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=cell_size, regular=False, qdegree=2) + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), + cellSize=cell_size, regular=False, qdegree=2) + + +def _vertex_coords(mesh): + """Local mesh vertex coordinates, in DM point order.""" + pStart, pEnd = mesh.dm.getDepthStratum(0) + raw = mesh.dm.getCoordinatesLocal().array.reshape(-1, mesh.cdim) + return raw[: pEnd - pStart], pStart + + +def _interior(points, margin): + """Keep points at least ``margin`` inside the unit box, so the query is an + interior FE evaluation and not an RBF extrapolation of a boundary point.""" + points = np.asarray(points) + keep = np.all((points > margin) & (points < 1.0 - margin), axis=1) + return np.ascontiguousarray(points[keep]), keep + + +class _CountedContainment: + """Count the containment tests a locator call performs. + + ``points`` counts point-tests, which is the quantity that blew up: the walk + re-tested the whole lost set on every one of its 50 rounds. + """ + + def __init__(self, mesh): + self.mesh = mesh + self.points = 0 + self.calls = 0 + self._wrapped = mesh._test_if_points_in_cells_internal + + def __enter__(self): + def counted(points, cells, **kwargs): + self.calls += 1 + self.points += len(points) + return self._wrapped(points, cells, **kwargs) + + self.mesh._test_if_points_in_cells_internal = counted + return self + + def __exit__(self, *exc): + self.mesh._test_if_points_in_cells_internal = self._wrapped + return False + + +@pytest.fixture(scope="module", params=[2, 3], ids=["2d", "3d"]) +def located_box(request): + """Mesh plus the ONE P1 variable every test shares. + + Adding a mesh variable rebuilds the DM (#492), so a fixture that handed + out a bare mesh and let each test add its own variable invalidated the + mesh under the tests that ran before it. + """ + dim = request.param + mesh = _box(dim, 1.0 / 16 if dim == 2 else 1.0 / 8) + field = uw.discretisation.MeshVariable(f"u_p1_{dim}", mesh, 1, degree=1) + field.data[:, 0] = _nodal_signal(field.coords) + mesh._build_kd_tree_index() + mesh._mark_faces_inside_and_out() + mesh._mark_local_boundary_faces_inside_and_out() + return dim, mesh, field + + +def _edge_endpoints(mesh, limit=1500): + """Vertex-coordinate pairs of local mesh edges.""" + verts, pStart = _vertex_coords(mesh) + eStart, eEnd = mesh.dm.getDepthStratum(1) + a, b = [], [] + for edge in range(eStart, min(eEnd, eStart + limit)): + cone = mesh.dm.getCone(edge) + a.append(verts[cone[0] - pStart]) + b.append(verts[cone[1] - pStart]) + return np.array(a, dtype=np.float64), np.array(b, dtype=np.float64) + + +# --------------------------------------------------------------------------- +# The hint contains the point (#432) +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("t", [0.0, 0.5, 0.25]) +def test_p1_field_on_a_shared_edge_matches_the_closed_form(located_box, t): + """``evaluate`` at a point on a cell edge must return ``(1-t)u_a + t u_b``. + + ``t = 0`` is a shared VERTEX, ``t = 0.5`` an edge midpoint. Both are shared + by every incident cell, so a locator that returns any cell CONTAINING the + point gives the exact answer, and a locator that returns a cell merely + NEAR it does not. + """ + dim, mesh, u = located_box + + a, b = _edge_endpoints(mesh) + query = (1.0 - t) * a + t * b + query, keep = _interior(query, 1.0e-4) + assert query.shape[0] > 100, "not enough interior edges to test with" + + expected = (1.0 - t) * _nodal_signal(a)[keep] + t * _nodal_signal(b)[keep] + got = np.asarray(uw.function.evaluate(u.sym[0], query)).reshape(-1) + + worst = float(np.abs(got - expected).max()) + assert worst < 1.0e-11, ( + f"{dim}-D P1 evaluation at t={t} along cell edges is off by {worst:.3e} " + f"at {int(np.count_nonzero(np.abs(got - expected) > 1.0e-11))} of " + f"{query.shape[0]} points — the query was answered in a cell that does " + f"not contain it (#432)") + + +def test_p1_field_on_a_shared_face_matches_the_closed_form(located_box): + """3-D face centroids: the P1 value is the mean of the three face vertices. + + A face is shared by exactly two tets and lies in the interior of neither; + it is the case the box clamp on reference coordinates cannot express. + """ + dim, mesh, u = located_box + if dim != 3: + pytest.skip("faces are edges in 2-D and are covered by the edge test") + + verts, pStart = _vertex_coords(mesh) + fStart, fEnd = mesh.dm.getHeightStratum(1) + corners = [] + for face in range(fStart, min(fEnd, fStart + 1500)): + closure = mesh.dm.getTransitiveClosure(face)[0][-3:] + corners.append(verts[closure - pStart]) + corners = np.array(corners, dtype=np.float64) # (n, 3 vertices, 3 coords) + + query = corners.mean(axis=1) + query, keep = _interior(query, 1.0e-4) + assert query.shape[0] > 100 + + expected = np.mean( + [_nodal_signal(corners[:, k, :])[keep] for k in range(3)], axis=0) + got = np.asarray(uw.function.evaluate(u.sym[0], query)).reshape(-1) + + worst = float(np.abs(got - expected).max()) + assert worst < 1.0e-11, ( + f"P1 evaluation at tet face centroids is off by {worst:.3e} — the " + f"query was answered in a cell that does not contain it (#432)") + + +def test_nodal_signal_discriminates_between_cells(located_box): + """NEGATIVE CONTROL for the two tests above. + + If the nodal field were smooth on the cell scale, the P1 interpolant + anchored to a neighbouring cell would agree with the right one and the + oracle would pass on a broken locator. Measure that the field is NOT + smooth on that scale: the closed-form edge value must differ from the + signal evaluated at the same point by far more than the tolerance the + tests assert. + """ + dim, mesh, _ = located_box + a, b = _edge_endpoints(mesh) + midpoint = 0.5 * (a + b) + interpolated = 0.5 * (_nodal_signal(a) + _nodal_signal(b)) + pointwise = _nodal_signal(midpoint) + spread = float(np.abs(interpolated - pointwise).max()) + assert spread > 0.1, ( + f"the {dim}-D nodal field varies by only {spread:.3e} between a cell " + f"edge and its midpoint — too smooth to detect a wrong cell") + + +def test_the_unchecked_nearest_control_point_hint_misses(located_box): + """NEGATIVE CONTROL for the fix in ``petsc_interpolate``. + + ``get_closest_cells`` is a nearest-CONTROL-POINT lookup with no containment + test; it was the hint the serial simplex path handed to the + DMLocatePoints bypass. Show it really does nominate cells that do not + contain the query — otherwise the edge test above would be pinning + nothing — and that the containment-checked locator does not. + + Quarter-points, not midpoints. A midpoint is equidistant from both ends of + the edge and the nearest control point is (measured) always one belonging + to a cell that does contain the edge; a quarter-point leans towards one + vertex, and the nearest control point is then any cell around THAT vertex, + most of which do not contain the far end. In 2-D the vertex neighbourhood + is small enough that even the quarter-point never misses on these meshes, + which is why #432 is a 3-D report. + """ + dim, mesh, _ = located_box + a, b = _edge_endpoints(mesh) + query, _ = _interior(0.75 * a + 0.25 * b, 1.0e-4) + + unchecked = np.asarray(mesh.get_closest_cells(query)).reshape(-1) + contained = mesh._test_if_points_in_cells_internal(query, unchecked) + + checked = np.asarray(mesh._robust_owning_cells(query)).reshape(-1) + assert (checked >= 0).all(), "an interior edge point was not located at all" + assert mesh._test_if_points_in_cells_internal( + query, checked, tol=mesh._EVAL_FACE_TOL).all(), ( + "the containment-checked locator returned a cell that does not " + "contain the point") + + if dim == 2: + pytest.skip( + f"the unchecked hint contains {int(contained.sum())} of " + f"{query.shape[0]} 2-D queries — nothing to demonstrate here; the " + f"3-D case carries the control") + + assert not contained.all(), ( + "the unchecked nearest-control-point hint happens to contain every " + "query on this mesh, so it cannot demonstrate the defect — pick a " + "harsher query set before trusting the edge test above") + + +def test_every_located_cell_contains_its_point(located_box): + """The locator's contract, which the shrinking working set must preserve. + + Which of several qualifying cells is returned for a point on a shared face + is not defined (and used to depend on whether some unrelated point in the + same batch was findable). That the returned cell CONTAINS the point is. + """ + dim, mesh, _ = located_box + rng = np.random.default_rng(4) + query = np.ascontiguousarray(rng.uniform(0.02, 0.98, size=(2000, dim))) + + for tol in (0.0, mesh._EVAL_FACE_TOL): + cells = np.asarray( + mesh._get_closest_local_cells_internal(query, tol=tol)).reshape(-1) + found = cells >= 0 + assert found.any(), "nothing was located at all" + assert mesh._test_if_points_in_cells_internal( + query[found], cells[found], tol=tol).all(), ( + f"a cell returned at tol={tol} does not contain its point") + + +# --------------------------------------------------------------------------- +# Rejection is cheap (#551 item 1) +# --------------------------------------------------------------------------- + +def test_a_point_outside_the_mesh_costs_a_bounded_number_of_tests(located_box): + """A point no local cell can own must be rejected in O(1) containment + tests, not by walking 50 nearest centroids. + + Structural, not a stopwatch: the count is what the algorithm does. Before + the fix this was exactly 51 tests per point at every mesh size. + """ + dim, mesh, _ = located_box + rng = np.random.default_rng(11) + outside = np.ascontiguousarray(rng.uniform(1.5, 2.5, size=(1000, dim))) + owned = np.ascontiguousarray(mesh._centroids[:1000]) + + with _CountedContainment(mesh) as counter: + cells = np.asarray(mesh._robust_owning_cells(outside)).reshape(-1) + assert (cells < 0).all(), "a point well outside the mesh was located" + per_outside_point = counter.points / outside.shape[0] + + with _CountedContainment(mesh) as counter: + cells = np.asarray(mesh._robust_owning_cells(owned)).reshape(-1) + assert (cells >= 0).all(), "a cell centroid was not located in its own cell" + per_owned_point = counter.points / owned.shape[0] + + # The counter fires at all: an owned point costs exactly the one test that + # confirms the nearest control point's cell. A zero here would mean the + # instrumentation missed the call and the bound above proved nothing. + assert per_owned_point >= 1.0 + assert per_owned_point <= 2.0 + + assert per_outside_point <= 4.0, ( + f"a point outside the {dim}-D mesh costs {per_outside_point:.1f} " + f"containment tests per point (was 51 before the rejection path)") + + +def test_one_unlocatable_point_does_not_cost_the_whole_batch(located_box): + """The working set has to shrink. + + The walk used to re-test every already-located point on every round, and + to keep going until the LAST point was found — so a single point that + could never be found charged the whole batch 50 extra rounds. + """ + dim, mesh, _ = located_box + rng = np.random.default_rng(12) + interior = np.ascontiguousarray(rng.uniform(0.02, 0.98, size=(1000, dim))) + poisoned = np.ascontiguousarray( + np.vstack([interior, rng.uniform(1.5, 2.5, size=(1, dim))])) + + with _CountedContainment(mesh) as counter: + mesh._robust_owning_cells(interior) + clean = counter.points + + with _CountedContainment(mesh) as counter: + mesh._robust_owning_cells(poisoned) + with_poison = counter.points + + assert with_poison <= clean + 60, ( + f"one unlocatable point added {with_poison - clean} containment " + f"point-tests to a batch of {interior.shape[0]}") + + +def test_the_classifier_hands_over_the_cells_it_located(located_box): + """#551 item 2: the classification and the interpolation share one + location pass, so the cells the classifier looked up come back with the + in/out mask instead of being thrown away.""" + dim, mesh, _ = located_box + rng = np.random.default_rng(13) + query = np.ascontiguousarray( + np.vstack([rng.uniform(0.02, 0.98, size=(500, dim)), + rng.uniform(1.5, 2.5, size=(50, dim))])) + + in_or_not, cells = mesh._classify_points_in_domain( + query, strict_validation=False) + + assert np.array_equal( + in_or_not, mesh.points_in_domain(query, strict_validation=False)), ( + "splitting the classifier changed points_in_domain's answer") + assert cells.shape[0] == query.shape[0] + # An exterior point is never given a hint: it goes to RBF extrapolation. + assert (cells[~in_or_not] == -1).all() + # Every hint that IS offered must be a cell containing its point. + offered = in_or_not & (cells >= 0) + if offered.any(): + assert mesh._test_if_points_in_cells_internal( + query[offered], cells[offered], tol=mesh._EVAL_FACE_TOL).all() From 60ab3d630a385f3722225a564cfb98c84a15b7f4 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 11:03:05 +1000 Subject: [PATCH 5/9] Fill the points the locator loses instead of returning NaN The RBF fallback in petsc_interpolate has never run. It iterates mesh.vars.values(), and mesh.vars is a weakref.WeakValueDictionary whose .values() is a GENERATOR, not a view: the dofcount loop a few lines earlier consumes it, so the fallback loop iterates nothing and the NaN that DMInterpolationEvaluate_UW writes for an unplaced point is returned to the caller. That did not matter until the previous commit in this branch, which routed serial simplex evaluation through the containment-checked locator. That locator returns -1 - the older nearest-control-point hint never did - so the dead rung became load-bearing on the most common configuration in the library. Measured on a graded 3-D simplex box (cellSize 1/8, deformed x -> x**4), 3000 interior queries: the merge base returns 0 NaN, this branch returned 1. Materialising the list fixes it: same probe, 0 NaN, and the point takes the RBF value (0.2802) the rung was always supposed to write. The same exhausted generator was read again by the continuity gate, so all() over it was vacuously True and the gate has never bound. Un-breaking it is a behaviour change, so it is measured rather than assumed. The gate only does anything on meshes whose measured location capability is "continuous" - warped hexes, the cubed-sphere class; simplex, quad, annulus and rectilinear hex boxes are all "exact", where the gate is a no-op either way. On a warped hex box carrying a P1 and a P0, binding it over EVERY variable on the mesh would cost the CONTINUOUS field a factor of 15 in accuracy (linear field, max error 8.0e-3 -> 1.2e-1 at 62 of 1500 interior points) because one discontinuous variable elsewhere took every evaluation off the authoritative path. Scoped instead to the variables the call actually asks for - which is what the surrounding comment always said the policy was - the continuous field is bit-identical to before and only the P0 moves, by up to 1.71 on a field of range 2, which is the O(jump) correction the gate exists to make. Also here, from the same review: - petsc_interpolate copies cell_hints before filling its -1 entries. np.ascontiguousarray hands back the caller's own array when it is already int64 and contiguous, and cell_hints is a documented keyword, so the fill was writing through somebody else's array. - _build_kd_tree_index_PIC and _build_kd_tree_index_DS are deleted. Both set _index without the new _local_cell_reach, which is the one way a stale rejection radius could outlive the geometry it was measured on; both have had zero callers for a long time and are already on the readability review's delete list (READ-31). - The tie-break change is written down where it can be found: which containing cell a shared vertex/edge/face query returns has moved, and for a discontinuous field the cell is the answer. - The absolute face-tolerance slab and the mesh-relative rejection radius are noted as the different scales they are, with the domain size at which they would cross (about 5e-6 across, measured clean at 1e-4 and at 6371). Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 123 +++++++----------- src/underworld3/function/_function.pyx | 38 +++++- 2 files changed, 78 insertions(+), 83 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 3bf5747c..77822112 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -5198,29 +5198,6 @@ def _get_coords_for_basis(self, degree, continuous): return arrcopy - def _build_kd_tree_index_DS(self): - - if hasattr(self, "_index") and self._index is not None: - return - - # Build this from the PETScDS rather than the SWARM - - centroids = self._get_coords_for_basis(0, False) - index_coords = self._get_coords_for_basis(2, False) - - points_per_cell = index_coords.shape[0] // centroids.shape[0] - - cell_id = numpy.empty(index_coords.shape[0]) - for i in range(cell_id.shape[0]): - cell_id[i] = i // points_per_cell - - self._indexCoords = index_coords - self._index = uw.kdtree.KDTree(self._indexCoords) - # self._index.build_index() - self._indexMap = numpy.array(cell_id, dtype=numpy.int64) - - return - def _coord_rows_for_points(self, nav_dm, points): """Row indices into the navigation coordinate array (``_nav_coords``) for the given vertex plex points, via the coordinate PetscSection @@ -5369,66 +5346,17 @@ def _build_kd_tree_index(self): self._centroid_index = uw.kdtree.KDTree(self._nav_centroids) # Rejection radius for the lost-point walk. Rebuilt with the kd-tree - # (i.e. invalidated by deform / adapt along with _index). + # (i.e. invalidated by deform / adapt along with _index) — this is the + # ONE place it is set, and the only builder of ``_index``, so a stale + # reach cannot outlive the geometry it was measured on. Two other + # builders (``_build_kd_tree_index_PIC``, ``_build_kd_tree_index_DS``) + # set ``_index`` without the reach; both had zero callers and are + # deleted rather than taught the new invariant. Pinned by + # test_0761_point_locator.py::test_the_rejection_radius_is_rebuilt_when_the_mesh_moves. self._local_cell_reach = cell_reach return - def _build_kd_tree_index_PIC(self): - - if hasattr(self, "_index") and self._index is not None: - return - - ## Bootstrapping - the kd-tree is needed to build the index but - ## the index is also used in the kd-tree. - - from underworld3.swarm import Swarm, SwarmPICLayout - - # Create a temp swarm which we'll use to populate particles - # at gauss points. These will then be used as basis for - # kd-tree indexing back to owning cells. - - from petsc4py import PETSc - - tempSwarm = PETSc.DMSwarm().create() - tempSwarm.setDimension(self.dim) - tempSwarm.setCellDM(self.dm) - tempSwarm.setType(PETSc.DMSwarm.Type.PIC) - tempSwarm.finalizeFieldRegister() - - # 3^dim or 4^dim pop is used. This number may need to be considered - # more carefully, or possibly should be coded to be set dynamically. - - tempSwarm.insertPointUsingCellDM(PETSc.DMSwarm.PICLayoutType.LAYOUT_GAUSS, 3) - - # We can't use our own populate function since this needs THIS kd_tree to exist - # We will need to use a standard layout instead - - ## ?? is this required given no migration ?? - # tempSwarm.migrate(remove_sent_points=True) - - PIC_coords = tempSwarm.getField("DMSwarmPIC_coor").reshape(-1, self.dim) - PIC_cellid = tempSwarm.getField("DMSwarm_cellid") - - self._indexCoords = PIC_coords.copy() - self._index = uw.kdtree.KDTree(self._indexCoords) - self._indexMap = numpy.array(PIC_cellid, dtype=numpy.int64) - # self._index.build_index() - - # We don't need an indexMap for this one because there is only one point per cell - # and the returned kdtree value IS the index. - # Note: self._centroids is not yet defined: - - self._centroid_index = uw.kdtree.KDTree(self._get_coords_for_basis(0, False)) - # self._centroid_index.build_index() - - tempSwarm.restoreField("DMSwarmPIC_coor") - tempSwarm.restoreField("DMSwarm_cellid") # - - tempSwarm.destroy() - - return - # Note - need to add this to the mesh rebuilding triggers def _facet_outward_unit_normal( self, facet_point_coords, facet_centroid, cell_centroid, @@ -6065,6 +5993,28 @@ def _get_closest_local_cells_internal( leaves the walk as soon as a cell claims it, so the walk costs what is still lost rather than what was asked for. + .. note:: **Which containing cell you get changed (#551).** + + A point on a shared vertex, edge or face is contained by several + cells and this routine returns one of them; *which* one has never + been part of the contract. It used to be the last cell to claim the + point across up to 50 rounds of the walk — an order that depended + on whether some unrelated point in the same batch was still lost. + It is now the containing cell with the nearest centroid, which is + batch-independent. Measured on a uniform 3-D simplex box, 35% of + near-vertex queries (the population that actually enters the walk; + exact vertices and centroids are answered before it) come back in a + different — equally containing — cell. + + For a CONTINUOUS field that is invisible: the interpolants of the + containing cells agree at the shared point. For a DISCONTINUOUS + field (P0, or the P2/P0-discontinuous pressure space the fault work + uses) the cell *is* the answer, so the evaluated value moves by + O(jump) at such points — measured max 1.935 on a P0 field of range + 2. Both values are legitimate: each is the value of a cell that + contains the query. Code that needs a specific side of a jump must + say which side, not rely on the locator's tie-break. + ``on_boundary`` and ``tol`` are forwarded to the in-cell containment test (see ``_test_if_points_in_cells_internal``). Default ``(on_boundary=True, tol=0.0)`` admits on-face queries @@ -6161,6 +6111,16 @@ def _get_closest_local_cells_internal( # further under that slab than a well-shaped one; a factor of two on # the reach covers both with room to spare while still rejecting # anything more than about one cell away from the local mesh. + # + # Note the two scales are set differently: the slab the containment + # test admits is ``tol`` times the face control-point separation, + # which _mark_faces_inside_and_out fixes at an ABSOLUTE 1e-3 in model + # units, while the radius here is a fraction of the LOCAL cell size. + # They only cross over when the largest local cell reach falls below + # about 5e-6 in model units — a whole domain a few microns across, at + # which scale the containment test's own absolute floors have already + # gone. Measured: a mesh 1e-4 across (reach 6.9e-6) and one 6371 + # across both reject nothing they should have kept. reach = getattr(self, "_local_cell_reach", None) if reach is not None and reach > 0.0: reach = self._LOCATOR_REACH_MARGIN * reach @@ -6354,6 +6314,11 @@ def _robust_owning_cells(self, coords: numpy.ndarray) -> numpy.ndarray: Never calls PETSc ``DMLocatePoints`` (slow, raises out-of-domain), and is purely kd-tree / Euclidean — manifold-safe, no manifold branch. + + For a point several cells share, the cell returned is the containing + one with the nearest centroid — see the tie-break note on + :meth:`_get_closest_local_cells_internal` for what that changed and + why it is visible only to discontinuous fields. """ coords = numpy.asarray(coords) if coords.shape[0] == 0: diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 7be624dd..6c785af0 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -1269,8 +1269,16 @@ def petsc_interpolate( expr, mesh._evaluation_interpolated_results = None - # For now, eval over all vars - vars = mesh.vars.values() + # For now, eval over all vars. + # + # MATERIALISE THE LIST. ``mesh.vars`` is a weakref.WeakValueDictionary + # and its ``.values()`` is a GENERATOR, not a view: the dofcount loop + # below consumes it, and every later reader (the continuity gate, the + # RBF fallback rung) then iterates an empty sequence. That silently + # disabled the fallback — points the locator returns -1 for kept the + # NaN written by DMInterpolationEvaluate_UW and handed it back to the + # caller — and it silently pinned the continuity gate at True. + vars = list(mesh.vars.values()) cdef DM dm = mesh.dm @@ -1305,7 +1313,23 @@ def petsc_interpolate( expr, # O(jump) wrong-side errors. The policy participates in the cache key # so the same coords evaluated with a different field mix cannot # reuse a structure built under the other policy. - all_continuous = all(getattr(var, "continuous", True) for var in vars) + # + # The continuity test runs over the variables THIS CALL ASKED FOR, + # not every variable on the mesh. The structure carries all of them + # (dofcount above), but only the requested slices are read, and the + # gate exists to protect a field whose jump sits on a cell face. Over + # the whole mesh instead, one discontinuous variable anywhere would + # take every evaluation off the authoritative path: measured on a + # warped hex box (capability "continuous") carrying a P1 and a P0, + # that costs the CONTINUOUS field a factor 15 in accuracy (linear + # field, max error 8.0e-3 -> 1.2e-1 at 62 of 1500 interior points) + # for no correctness gain. Scoped to the request, the continuous + # field is bit-identical and only the P0 moves. + # + # NOTE this gate has never bound before: `vars` was an exhausted + # generator (see above) so `all()` was vacuously True. + all_continuous = all( + getattr(varfn.meshvar(), "continuous", True) for varfn in varfns) authoritative = mesh._hint_is_authoritative(all_continuous) location_policy = "auth" if authoritative else "locate" @@ -1351,7 +1375,13 @@ def petsc_interpolate( expr, # one location per point per call rather than two. if authoritative: if cell_hints is not None and mesh is hinted_mesh: - cells = np.ascontiguousarray(cell_hints, dtype=np.int64) + # COPY: the unhinted entries are filled in below, and + # ascontiguousarray hands back the caller's own array when + # it is already int64 and contiguous. petsc_interpolate + # takes cell_hints as a documented keyword, so writing + # through it would mutate somebody else's array. + cells = np.array(cell_hints, dtype=np.int64, copy=True, + order="C") if cells.shape[0] != coords.shape[0]: raise RuntimeError( "cell_hints must carry one cell index per coordinate " From 61326ee66b07769e9d06249f23f9f2acd95632c3 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 11:03:43 +1000 Subject: [PATCH 6/9] Test the direction of the rejection radius that would be silent Nothing asserted that a point a local cell really does contain is never rejected by the new radius, which is the failure mode that would not announce itself: the point vanishes into the RBF fallback in serial, or is claimed by nobody in parallel. Four tests, all with their own negative control. A brute-force oracle - every query against every local cell, with the same containment predicate the locator uses - says which points are genuinely owned, and the reach margin is swept to prove the probe fires. Serial: 2-D {2.0: 0, 1.0: 0, 0.5: 5, 0.25: 16}, 3-D {2.0: 0, 1.0: 0, 0.5: 4, 0.25: 168} out of 800 interior points. The bound is tight at about 1.0, so the shipped 2.0 carries a factor of two over the first observable loss. 0 at the shipped margin on every rank at np2 and np4 as well. A P0 discontinuous field evaluated at shared mesh vertices pins the tie-break: the value must belong to one of the cells that contains the query, not a specific cell (that was never the contract, and pinning it would break on the next legitimate change). Every vertex in the set is contained by 18 to 44 cells and their values differ by more than 0.5, so a wrong cell is visible. The RBF fallback rung gets a fault-injection test: five points are forced to -1 and the result must be finite AND equal the RBF value AND leave the other 395 points bit-identical. Waiting for a graded mesh to lose a point naturally works, but how many it loses is a property of whatever gmsh produced that day. Confirmed fail-before by rebuilding with the one-line generator fix reverted: exactly the five injected points come back NaN. The reach invalidation was asserted in prose only. Deform a 2-D box by 50 and the stored reach must follow exactly, and every point in the expanded domain must still be located. Finally, test_the_classifier_hands_over_the_cells_it_located ended with its real assertion inside `if offered.any():`, and offered is ALWAYS empty in serial - 0 hints for 1465 interior points, against 462 of 462 at np4 - so the strongest claim in item 2's only test was a no-op in the default test run. It now asserts the serial reality (no hints, by design) and requires hints under MPI. Underworld development team with AI support from Claude Code --- tests/test_0761_point_locator.py | 320 ++++++++++++++++++++++++++++++- 1 file changed, 315 insertions(+), 5 deletions(-) diff --git a/tests/test_0761_point_locator.py b/tests/test_0761_point_locator.py index 9a27f909..874edd0c 100644 --- a/tests/test_0761_point_locator.py +++ b/tests/test_0761_point_locator.py @@ -2,7 +2,7 @@ cell can own must be rejected in O(1) (#551 items 1 and 3; #432, a recurrence of #390). -Two properties are pinned here, both of which have been broken before. +Five properties are pinned here, all of which have been broken before. **The hint contains the point.** ``petsc_interpolate`` may bypass PETSc's ``DMLocatePoints`` and evaluate the basis directly in a cell UW3 nominates. @@ -23,6 +23,26 @@ containment tests rather than seconds: an instrumented counter says what the algorithm does, where a stopwatch says what the machine was doing at the time. The parallel wall-clock table lives in the PR, not in a brittle timing assert. + +**Rejection never costs a point the mesh owns.** The cheap direction of the +rejection radius is easy to see; the expensive one is silent. A brute-force +oracle (every point against every local cell) says which points a local cell +really contains, and a sweep over the reach margin is the negative control: +0 false rejections at the shipped margin and at half of it, and a growing +count below that, so a clean result is a measurement and not a coincidence. + +**Which containing cell you get is allowed to change; a NON-containing cell +is not.** #551 moved the tie-break for a point several cells share from "last +cell to claim it across up to 50 walk rounds" (which depended on the rest of +the batch) to "containing cell with the nearest centroid". For a discontinuous +field the cell is the answer, so this is user-visible; the test pins the +property that survives — the value belongs to a cell containing the query — +rather than a specific cell, which was never part of the contract. + +**A point no cell owns gets a value, not a NaN.** The RBF fallback rung in +``petsc_interpolate`` used to iterate an exhausted generator and fill nothing. +The loss is injected rather than waited for, so the rung is exercised on every +mesh and every platform. """ import numpy as np import pytest @@ -114,6 +134,25 @@ def located_box(request): return dim, mesh, field +def _cells_containing(mesh, points, tol): + """Brute-force oracle: for every point, which local cells contain it? + + Tests every point against every local cell with the SAME containment + predicate the locator uses, so a disagreement is a locator disagreement + and not a difference of geometric opinion. Returns a boolean + ``(n_points, n_cells)`` table. + """ + nav_dm = mesh._nav_dm if mesh._nav_dm is not None else mesh.dm + cStart, cEnd = nav_dm.getHeightStratum(0) + table = np.zeros((points.shape[0], cEnd - cStart), dtype=bool) + column = np.empty(points.shape[0], dtype=np.int64) + for cell in range(cEnd - cStart): + column[:] = cell + table[:, cell] = mesh._test_if_points_in_cells_internal( + points, column, tol=tol) + return table + + def _edge_endpoints(mesh, limit=1500): """Vertex-coordinate pairs of local mesh edges.""" verts, pStart = _vertex_coords(mesh) @@ -274,6 +313,258 @@ def test_every_located_cell_contains_its_point(located_box): f"a cell returned at tol={tol} does not contain its point") +def test_a_point_a_local_cell_contains_is_never_rejected(located_box): + """The dangerous direction of the rejection radius. + + Rejecting a foreign point early is the point of the radius; rejecting a + point that a local cell really does contain would be silent — the point + disappears into the RBF fallback (serial) or is claimed by nobody + (parallel), and no other test in this file would notice. + + The oracle is brute force: every point against every local cell, with the + same containment predicate the locator uses. The sweep over + ``_LOCATOR_REACH_MARGIN`` is the negative control — it shows the probe + detects false rejections when the radius is deliberately too tight, so a + clean result at the shipped margin means something. + """ + dim, mesh, _ = located_box + rng = np.random.default_rng(21) + query = np.ascontiguousarray(rng.uniform(0.02, 0.98, size=(800, dim))) + tol = mesh._EVAL_FACE_TOL + + # Rank-local: in parallel a rank only holds its share of the query, so the + # bar scales with the partition (measured 179-434 of 800 per rank at + # np2/np4, 800 of 800 in serial). + contained = _cells_containing(mesh, query, tol).any(axis=1) + assert contained.sum() > 0.4 * query.shape[0] / uw.mpi.size, ( + f"only {int(contained.sum())} of {query.shape[0]} interior points are " + f"contained by a cell of this rank — the oracle is not exercising the " + f"local mesh") + + shipped = mesh._LOCATOR_REACH_MARGIN + counts = {} + try: + for margin in (shipped, 1.0, 0.5, 0.25, 0.1): + mesh._LOCATOR_REACH_MARGIN = margin + got = np.asarray(mesh._get_closest_local_cells_internal( + query, tol=tol)).reshape(-1) + counts[margin] = int(np.count_nonzero(contained & (got < 0))) + finally: + del mesh._LOCATOR_REACH_MARGIN + assert mesh._LOCATOR_REACH_MARGIN == shipped + + # Measured serial: 2-D {2.0: 0, 1.0: 0, 0.5: 5, 0.25: 16, 0.1: 16}, + # 3-D {2.0: 0, 1.0: 0, 0.5: 4, 0.25: 168, 0.1: 237} — the bound is tight + # at about 1.0 and the shipped 2.0 keeps a factor of two over the first + # observable loss. 0 at 2.0 and 1.0 on every rank at np2 and np4 too. + assert counts[shipped] == 0, ( + f"the shipped reach margin ({shipped}) rejected {counts[shipped]} of " + f"{int(contained.sum())} {dim}-D points that a local cell contains " + f"(sweep: {counts})") + assert counts[1.0] == 0 + + # NEGATIVE CONTROL. A quarter of the true reach must lose points, or the + # clean result above is the probe failing to fire rather than the radius + # being right. + assert counts[0.25] > 0, ( + f"a reach margin of 0.25 rejected nothing in {dim}-D, so this test " + f"cannot tell a correct radius from a broken one (sweep: {counts})") + + +def test_the_rejection_radius_is_rebuilt_when_the_mesh_moves(): + """A stale SMALL reach is the one way the rejection silently loses points. + + ``_local_cell_reach`` is measured in :meth:`_build_kd_tree_index` and + invalidated with ``_index``. That is an invariant, not a comment: deform + the mesh so its cells grow, and both the stored reach and the locations it + admits must follow. + """ + mesh = _box(2, 1.0 / 10) + mesh._build_kd_tree_index() + mesh._mark_faces_inside_and_out() + before = mesh._local_cell_reach + assert before > 0.0 + + mesh.deform(np.asarray(mesh.X.coords, dtype=np.float64) * 50.0) + mesh._build_kd_tree_index() + mesh._mark_faces_inside_and_out() + after = mesh._local_cell_reach + + assert after / before == pytest.approx(50.0, rel=1e-9), ( + f"reach {before:.6g} -> {after:.6g} after deform(x50): the rejection " + f"radius did not follow the geometry") + + # And the expanded domain is actually usable: every point a cell contains + # is still located, none rejected by a radius measured on the old mesh. + rng = np.random.default_rng(22) + query = np.ascontiguousarray(rng.uniform(1.0, 49.0, size=(400, 2))) + tol = mesh._EVAL_FACE_TOL + contained = _cells_containing(mesh, query, tol).any(axis=1) + got = np.asarray(mesh._get_closest_local_cells_internal( + query, tol=tol)).reshape(-1) + assert contained.any(), "no query landed in the deformed mesh" + assert not (contained & (got < 0)).any(), ( + f"{int((contained & (got < 0)).sum())} points inside the deformed mesh " + f"were rejected — the reach did not grow with the cells") + + del mesh + + +def test_a_discontinuous_field_at_a_shared_vertex_takes_a_containing_cell(): + """#551 changed WHICH containing cell a shared point gets, and for a + discontinuous field the cell IS the answer. + + A mesh vertex is contained by every cell around it. The locator returns + one of them, and the tie-break moved from "last cell to claim it across up + to 50 walk rounds" (batch-dependent) to "containing cell with the nearest + centroid". Pinning a SPECIFIC cell would over-pin: it is not part of the + contract, and the ambiguity is genuine. What is pinned is that the value + comes from a cell that CONTAINS the point — so a future change that + answers from a merely nearby cell fails here. + + The P2/P0-discontinuous pressure space the fault work uses is exactly this + configuration, which is why it is worth a test rather than a paragraph. + """ + if uw.mpi.size > 1: + pytest.skip( + "serial only: the oracle is the LOCAL cell table, and in parallel " + "the rank that answers a seam vertex need not be the rank whose " + "cells this rank can see") + + mesh = _box(3, 1.0 / 6) + p0 = uw.discretisation.MeshVariable( + "p0_jump", mesh, 1, degree=0, continuous=False) + mesh._build_kd_tree_index() + mesh._mark_faces_inside_and_out() + mesh._mark_local_boundary_faces_inside_and_out() + + rng = np.random.default_rng(23) + p0.data[:, 0] = rng.uniform(-1.0, 1.0, size=p0.data.shape[0]) + + # Map nav-DM cell index -> P0 degree of freedom by centroid identity, and + # check the map is a bijection before trusting it. + centroids = np.asarray(mesh._nav_centroids, dtype=np.float64) + dof_coords = np.asarray(p0.coords, dtype=np.float64) + separation = np.linalg.norm( + centroids[:, None, :] - dof_coords[None, :, :], axis=-1) + dof_of_cell = separation.argmin(axis=1) + assert separation.min(axis=1).max() < 1.0e-10 + assert len(set(dof_of_cell.tolist())) == centroids.shape[0] + + verts, _ = _vertex_coords(mesh) + query, _ = _interior(verts, 0.05) + query = np.ascontiguousarray(query[:200]) + assert query.shape[0] > 50, "not enough interior vertices to test with" + + tol = mesh._EVAL_FACE_TOL + table = _cells_containing(mesh, query, tol) + shared = table.sum(axis=1) + assert (shared >= 2).mean() > 0.9, ( + f"only {(shared >= 2).mean():.2f} of the vertex queries are contained " + f"by more than one cell — nothing here is ambiguous, so the test " + f"cannot see a tie-break at all") + + got = np.asarray(uw.function.evaluate(p0.sym[0], query)).reshape(-1) + assert np.isfinite(got).all() + + # NEGATIVE CONTROL: the containing cells must actually disagree, otherwise + # "the value is one of them" is satisfied by any locator at all. + spread = np.array([ + np.ptp(p0.data[dof_of_cell[np.flatnonzero(row)], 0]) for row in table]) + assert spread.max() > 0.5, ( + f"the containing cells of every vertex agree to within " + f"{spread.max():.3e} — a wrong cell would be undetectable here") + + for i, row in enumerate(table): + allowed = p0.data[dof_of_cell[np.flatnonzero(row)], 0] + assert np.abs(allowed - got[i]).min() < 1.0e-12, ( + f"P0 evaluation at a shared vertex returned {got[i]:.6g}, which is " + f"not the value of any of the {int(row.sum())} cells containing it " + f"({np.sort(allowed)}) — the query was answered in a cell that does " + f"not contain it") + + del mesh + + +def test_evaluate_fills_a_point_no_cell_owns_instead_of_returning_nan(): + """The RBF fallback rung has to be reachable. + + ``petsc_interpolate`` may hand ``DMInterpolation`` a cell hint of ``-1`` + for a point it could not place; ``DMInterpolationEvaluate_UW`` writes NaN + there, and the rung below is what replaces it with the bounded, topology- + free RBF value. The rung iterated ``mesh.vars.values()``, a WeakValueDict + GENERATOR that an earlier loop had already exhausted, so it filled nothing + and the NaN was returned to the caller. + + The loss is injected rather than waited for: a graded mesh does lose the + occasional interior point to the 50-neighbour cap, but how many is a + property of whatever gmsh produced today. Injecting it makes the test + exercise the rung on every mesh, every platform, every run. + """ + if uw.mpi.size > 1: + pytest.skip( + "serial only: with more than one rank the injected loss re-routes " + "the point to whichever rank still claims it, so the value is not " + "this rank's RBF interpolant. Serial is where #551 item 3 opened " + "the door onto this rung") + + mesh = _box(3, 1.0 / 8) + u = uw.discretisation.MeshVariable("u_fallback", mesh, 1, degree=1) + coords = np.asarray(mesh.X.coords, dtype=np.float64).copy() + coords[:, 0] = coords[:, 0] ** 4 # strong grading, as #551 reports + mesh.deform(coords) + mesh._build_kd_tree_index() + mesh._mark_faces_inside_and_out() + mesh._mark_local_boundary_faces_inside_and_out() + u.data[:, 0] = _nodal_signal(np.asarray(u.coords, dtype=np.float64)) + + rng = np.random.default_rng(24) + query = np.ascontiguousarray(rng.uniform(0.02, 0.98, size=(400, 3))) + intact = np.asarray(uw.function.evaluate(u.sym[0], query)).reshape(-1) + assert np.isfinite(intact).all(), ( + f"{int((~np.isfinite(intact)).sum())} of {query.shape[0]} interior " + f"points came back NaN with the locator untouched") + + victims = np.array([3, 17, 61, 128, 349]) + victim_coords = query[victims] + located = mesh._robust_owning_cells + mesh._dminterpolation_cache.invalidate_all("fault injection") + + def loses_the_victims(points): + cells = np.asarray(located(points), dtype=np.int64).copy() + for victim in victim_coords: + cells[np.all(points == victim, axis=1)] = -1 + return cells + + mesh._robust_owning_cells = loses_the_victims + try: + got = np.asarray(uw.function.evaluate(u.sym[0], query)).reshape(-1) + finally: + del mesh._robust_owning_cells + mesh._dminterpolation_cache.invalidate_all("fault injection") + + assert np.isfinite(got).all(), ( + f"{int((~np.isfinite(got)).sum())} points the locator could not place " + f"came back NaN: the RBF fallback rung did not run") + + expected = np.asarray( + u.rbf_interpolate(np.ascontiguousarray(victim_coords))).reshape(-1) + assert np.abs(got[victims] - expected).max() < 1.0e-12, ( + "the unlocated points were filled with something other than the RBF " + "interpolant") + + # NEGATIVE CONTROL: the fallback value must be distinguishable from the FE + # value, or "no NaN" could be passing on an untouched array. + assert np.abs(got[victims] - intact[victims]).max() > 1.0e-6, ( + "the injected loss changed nothing — the fault injection missed and " + "the rung was never asked to run") + others = np.setdiff1d(np.arange(query.shape[0]), victims) + assert np.abs(got[others] - intact[others]).max() < 1.0e-12, ( + "injecting a loss at five points changed the answer elsewhere") + + del mesh + + # --------------------------------------------------------------------------- # Rejection is cheap (#551 item 1) # --------------------------------------------------------------------------- @@ -356,8 +647,27 @@ def test_the_classifier_hands_over_the_cells_it_located(located_box): assert cells.shape[0] == query.shape[0] # An exterior point is never given a hint: it goes to RBF extrapolation. assert (cells[~in_or_not] == -1).all() - # Every hint that IS offered must be a cell containing its point. + + # Which points get a hint is decided by _eval_use_robust_location(), and + # that is False in serial by design: the serial classifier keeps the + # validated cell-wall path bit-for-bit and offers NO cells at all. Assert + # that reality rather than wrapping the real check in `if offered.any():` + # — measured 0 hints for 1465 interior points in serial against 462 of 462 + # at np4, so the guarded version asserted nothing in the default test run. offered = in_or_not & (cells >= 0) - if offered.any(): - assert mesh._test_if_points_in_cells_internal( - query[offered], cells[offered], tol=mesh._EVAL_FACE_TOL).all() + if uw.mpi.size == 1: + assert not mesh._eval_use_robust_location() + assert not offered.any(), ( + f"the serial classifier handed over {int(offered.sum())} cells; it " + f"is supposed to hand over none, and the evaluator locates the " + f"points itself") + return + + assert mesh._eval_use_robust_location() + assert offered.any(), ( + f"the parallel classifier located points in the domain " + f"({int(in_or_not.sum())} of {query.shape[0]}) but handed over no " + f"cells, so item 2 saves nothing") + # Every hint that IS offered must be a cell containing its point. + assert mesh._test_if_points_in_cells_internal( + query[offered], cells[offered], tol=mesh._EVAL_FACE_TOL).all() From de87c7e7cfdc447f7517aefa1da91e587d3537d4 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 11:32:13 +1000 Subject: [PATCH 7/9] Make the deform-invalidation test say what went wrong when it fails PR #556 failed CI on test_rotated_workspace_deform_invalidates with "post-deform solve differs from fresh control by 1.16e-01", and the failure does not reproduce here. Not in the amr-dev toolchain (custom AMR PETSc), not in the dev toolchain CI itself uses (conda PETSc), not as a single test, not as the whole file, not in the exact CI batch shape (pytest tests/test_101*py tests/test_102*py: 137 passed in both environments), with or without CI's environment variables. Locally err is not small, it is exactly 0.0: one Krylov iteration, |r| = 8.267546030330813e-12 bit-identical between the post-deform solve and the fresh control. The locator work in this PR is measurably inert in this test: - the mesh is a StructuredQuadBox, and for a quad/hex mesh the merge base ALREADY located through _robust_owning_cells - item 3 changed the serial SIMPLEX path only, so the evaluate route here is unchanged; - there are 8 point-location calls during the test (from swarm.migrate -> points_in_domain and from global_evaluate, not from the assembly). Running every one of them twice, with the new rejection radius and with it disabled, gives identical cells: 0 differences out of 2270 points; - the deformed quad box measures location capability "exact" before and after the deform, so _hint_is_authoritative is True whatever the field continuity and the continuity gate cannot bind; - in serial the classifier hands over no cell hints at all, so item 2 is inert too; - the locator reach follows the deform (0.070710678 -> 0.073041055). Rather than guess at a fix for something that cannot be reproduced, the test now distinguishes the two things it could be. Both solvers assemble the same system on the same deformed mesh from a zero guess, so any deterministic solver owes them the same answer: if both report converged and the answers still differ, the OPERATORS differ and something survived the deform; if the convergence reports differ, the linear solves did. The failure message carries both solves' reason, iteration counts and residuals plus the locator reach before and after, and the reach invalidation is now asserted here as well as in test_0761 (the review's minor m1) because this is the only test in the suite that deforms a mesh between two solves. Underworld development team with AI support from Claude Code --- tests/test_1018_rotated_freeslip.py | 56 +++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index c7fba4b7..aea2b930 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -165,6 +165,19 @@ def test_rotated_linear_workspace_reuses_unchanged_operator(): assert time_refresh_error < 1.0e-6 +def _solve_report(info): + """The rotated solve's own verdict, for failure messages: an assertion + that says only "the answer moved" cannot tell a stale operator from a + linear solve that stopped early.""" + return (f"converged={info.get('converged')} " + f"ksp_reason={info.get('ksp_reason')} " + f"ksp_its={info.get('ksp_its')} " + f"newton_its={info.get('nonlinear_iterations')} " + f"|r|={info.get('rnorm')} |r0|={info.get('rnorm0')} " + f"rotation_reused={info.get('rotation_reused')} " + f"workspace_reused={info.get('workspace_reused')}") + + def _rampable_rotated_stokes(mesh, k_expr, tag, forcing=None): """Rotated free-slip Stokes with viscosity given by a rampable UWexpression constant (the #416 idiom used by every continuation @@ -244,7 +257,18 @@ def test_rotated_workspace_constant_ramp_invalidates(): def test_rotated_workspace_deform_invalidates(): """mesh.deform between solves: geometry changed, so the whole workspace must be rebuilt (rotation_reused False) and the answer must match a fresh - solver on the deformed mesh.""" + solver on the deformed mesh. + + Both solvers assemble the same system on the same deformed mesh from the + same (zero) initial guess, so any deterministic linear solver has to give + them the same answer: a non-zero ``err`` means one of the two systems is + not what it should be, not that a Krylov path drifted. The checks below + say WHICH, because this test has failed in CI on a platform where it + cannot be reproduced locally (macOS arm64, both the conda-PETSc `dev` and + the AMR-PETSc `amr-dev` toolchains, single test / whole file / CI batch: + err is exactly 0.0, one Krylov iteration, |r| = 8.3e-12 identical between + the two solves). Reporting the solver's own diagnostics and the mesh's + geometry invariants turns the next occurrence into an attribution.""" mesh = uw.meshing.StructuredQuadBox( elementRes=(10, 10), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) @@ -252,24 +276,52 @@ def test_rotated_workspace_deform_invalidates(): s1, v1 = _rampable_rotated_stokes(mesh, k, "Dfm") s1.solve() assert s1._rotated_linear_cache is not None + info_1 = dict(s1._rotated_freeslip_info) + reach_before = mesh._local_cell_reach # bump the top boundary (the free-surface pattern) coords = mesh.X.coords.copy() coords[:, 1] += 0.02 * coords[:, 1] * np.sin(np.pi * coords[:, 0]) mesh.deform(coords) + # Mesh-side invariant: the kd-tree index and the point locator's rejection + # radius (#551) are measured together in _build_kd_tree_index, which + # deform() drops and rebuilds eagerly. A stale SMALL radius is the one way + # point location silently loses points, and this is the only test in the + # suite that deforms a mesh between two solves, so it is the one that + # would see it. + reach_after = mesh._local_cell_reach + assert reach_after != reach_before, ( + f"the locator reach did not follow the deform " + f"({reach_before:.8g} -> {reach_after:.8g}) — it is measured with the " + f"kd-tree index and must be rebuilt with it") + s1.solve() info = s1._rotated_freeslip_info assert not info["rotation_reused"], ( "workspace survived a mesh.deform — stale rotation Q in use") assert not info["workspace_reused"] + assert info["converged"], ( + f"the post-deform solve did not converge: {_solve_report(info)}") k_c = uw.function.expression(r"k_dc", 1.0, "control viscosity") s_c, v_c = _rampable_rotated_stokes(mesh, k_c, "DfC") s_c.solve() + info_c = s_c._rotated_freeslip_info + assert info_c["converged"], ( + f"the fresh control solve did not converge: {_solve_report(info_c)}") + err = np.linalg.norm(v1.data - v_c.data) / np.linalg.norm(v_c.data) assert err < 1e-6, ( - f"post-deform solve differs from fresh control by {err:.2e}") + f"post-deform solve differs from fresh control by {err:.2e}\n" + f" first solve : {_solve_report(info_1)}\n" + f" post-deform : {_solve_report(info)}\n" + f" fresh control : {_solve_report(info_c)}\n" + f" locator reach : {reach_before:.8g} -> {reach_after:.8g}\n" + f"Both solves assemble the same system from a zero guess, so equal " + f"convergence reports with a non-zero err means the two OPERATORS " + f"differ — something survived the deform. Different |r| or iteration " + f"counts means the linear solves themselves diverged.") @pytest.mark.level_2 From e8ec403b1a11a8eecd4de7ba2ea4605f73c52b8b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 12:58:57 +1000 Subject: [PATCH 8/9] Record which rigid-body modes the rotated gauge admits, and decompose the difference CI's first diagnosis-grade report on test_rotated_workspace_deform_invalidates said both solves converged from an IDENTICAL initial residual (|r0| = 0.02785134866629431 to every digit, |r| ~ 1e-11 for both) and still disagreed by 7.02e-02. Two solves that both drive the residual to machine zero on the same system can only differ in the null space, so this round measures the null space instead of arguing about it. Solver side: the rigid-rotation gauge decision was invisible. Whether a mode is admitted decides whether a component is projected out of the answer, and _mode_satisfies_constraints made that call silently, per mode, from the boundary normals. It now optionally records the constraint violation, the operator violation and the verdict; _finalize_rotated_solution collects one record per offered mode and the solve result carries it as "rotation_gauge". Default arguments unchanged, so nothing else moves. Test side, all of it reported in the failure message and all of it computed EAGERLY, because an instrument that only runs when the test fails is an instrument that has never been run: - the difference field is decomposed onto the rigid-body span (the same modes the solver considers, plus translations, Gram-Schmidt'd in the same order, built from nodal coordinates so it does not inherit the machinery it is measuring). If the difference lives in the span, one solve admitted a mode the other rejected and this is a #543 gauge bug; if it lives off the span, the two operators differ and something survived the deform; - per-solver constrained-row count, distinct-row count and boundary list; - a point-location tally around the deform, the post-deform solve and the control, answering every call TWICE - as shipped and with the rejection radius set aside - and reporting how many answers the radius changed. That is the number that decides whether #556 is implicated, taken on whatever mesh the machine actually built rather than argued from a local run. Three negative controls, so none of those numbers is unchecked: a pure rigid rotation must decompose with 4e-18 of it off the span and a random field with 0.9985 of it off (both asserted on the mesh under test), and squeezing the reach margin to 0.05 must make the radius comparison see answers change. Measured here: 84 constrained rows for both solvers, one rotation mode offered and rejected by both with identical violation 3.045e-01, gauge_removed False for both, and during the deform 8 location calls over 2270 points with 38 returning -1 and radius_changed = 0. The locator changes nothing on this platform. If CI reports otherwise, that is the finding. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/rotated_bc.py | 46 ++++- tests/test_1018_rotated_freeslip.py | 227 ++++++++++++++++++++++-- 2 files changed, 251 insertions(+), 22 deletions(-) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 61ae92a3..0882f854 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -593,10 +593,18 @@ def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge) # three 3D modes are not mutually orthogonal on a general mesh. # COLLECTIVE: all ranks walk the same mode list, same order. removed = False + # What the gauge decision actually was, recorded for the caller. Two solves + # of the same system that both converge can still differ by an admitted + # gauge mode, and until this was recorded there was no way to see which of + # them admitted what. + gauge = {"offered": 0, "modes": [], "orthonormal": 0, "removed": False, + "requested": bool(remove_rotation_gauge)} if remove_rotation_gauge: live = [] for tg in _rigid_rotation_modes(solver): - if _mode_satisfies_constraints(solver, Q, normal_rows, tg): + gauge["offered"] += 1 + if _mode_satisfies_constraints(solver, Q, normal_rows, tg, + report=gauge["modes"]): live.append(tg.copy()) # owned copy — pool vec goes back below dm.restoreGlobalVec(tg) ortho = [] @@ -609,10 +617,13 @@ def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge) ortho.append(w) else: w.destroy() + gauge["orthonormal"] = len(ortho) for q in ortho: U.axpy(-U.dot(q), q) q.destroy() removed = True + gauge["removed"] = removed + solver._rotated_gauge_report = gauge # scatter U → velocity/pressure fields, completing each field's essential # (Dirichlet) DOFs. Those are absent from the global vector, so a plain @@ -769,6 +780,11 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, * ``"normal_rows"`` — global rows of the constrained normal components; * ``"boundaries"`` — the boundary specs the rotation was built from; * ``"rotation_gauge_removed"`` — whether a rigid-rotation gauge was projected out; + * ``"rotation_gauge"`` — how that was decided: how many rigid-body modes + were offered, and per mode the constraint violation, the operator + violation and the verdict. Two solves of the same system that both + converge can still differ by an admitted gauge mode, and this is the + only place that decision is visible; * ``"ksp_reason"`` — converged-reason of the LAST Newton increment's KSP; * ``"ksp_its"`` — list of linear iteration counts, one per Newton iteration; * ``"nonlinear_iterations"``, ``"converged"`` — outer-loop count (the number @@ -1307,7 +1323,11 @@ def rotated_residual(uvec, keep_cartesian=False): return {"Q": Q, "Qt": Qt, "reaction": reaction, "U": u, "normal_rows": normal_rows, "boundaries": list(boundaries), - "rotation_gauge_removed": removed, "ksp_reason": last_reason, + "rotation_gauge_removed": removed, + # per-mode record of the gauge decision (offered / accepted / the + # violations that decided it) — see _finalize_rotated_solution + "rotation_gauge": getattr(solver, "_rotated_gauge_report", None), + "ksp_reason": last_reason, "nonlinear_iterations": newton_its, "converged": converged, "ksp_its": lin_its, "rnorm": rnorm, "rnorm0": r0, "vel_its_last": vel_its_last, "pres_its_last": pres_its_last, @@ -1726,7 +1746,8 @@ def _rotated_nullspace(solver, Q, normal_rows): return PETSc.NullSpace().create(constant=False, vectors=ortho, comm=dm.comm) -def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): +def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8, + report=None): """True iff the rigid-body mode ``tg`` is a genuine null mode of the constrained problem: it satisfies all rotated v_n=0 constraints (Q·tg ~0 on every constrained normal row) AND it is a null vector of the assembled @@ -1748,7 +1769,14 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): COLLECTIVE: every rank runs the same global-vector ops. Do NOT early-return on a per-rank ``not normal_rows`` — in parallel a rank may own no boundary node (empty normal_rows) while others do, and an early return there would desync the - collective norms below and deadlock.""" + collective norms below and deadlock. + + ``report``, if given, is a list this appends one dict to per call: the + constraint violation, the operator violation (``None`` when the constraint + test already rejected the mode) and the verdict. Whether a mode is admitted + decides whether a gauge component is projected out of the answer, so when + two solves of the same system disagree this is the first thing to look at + — but nothing in the returned dict was observable before.""" tr = tg.duplicate() Q.mult(tg, tr) full = tr.norm() # parallel norm @@ -1766,6 +1794,9 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): tr.destroy() trc.destroy() if viol >= tol: + if report is not None: + report.append({"viol": float(viol), "op_viol": None, + "accepted": False, "rejected_by": "constraint"}) return False # operator-nullity: rigid rotation has exactly zero strain in the discrete # space (P2 contains linear fields, affine cells integrate the form exactly), @@ -1777,8 +1808,15 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): jn = Jm.norm() Jm.destroy() if jn == 0.0 and J.norm() == 0.0: # J never assembled → cannot verify + if report is not None: + report.append({"viol": float(viol), "op_viol": None, + "accepted": True, "rejected_by": "unassembled-J"}) return True op_viol = jn / (_velocity_diag_scale(J, solver) * (tg.norm() + 1e-30)) + if report is not None: + report.append({"viol": float(viol), "op_viol": float(op_viol), + "accepted": bool(op_viol < tol), + "rejected_by": None if op_viol < tol else "operator"}) return op_viol < tol diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index aea2b930..92c91e80 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -168,14 +168,140 @@ def test_rotated_linear_workspace_reuses_unchanged_operator(): def _solve_report(info): """The rotated solve's own verdict, for failure messages: an assertion that says only "the answer moved" cannot tell a stale operator from a - linear solve that stopped early.""" + linear solve that stopped early, or either of those from a null-space + gauge that one solve removed and the other did not.""" + rows = info.get("normal_rows") + gauge = info.get("rotation_gauge") or {} + modes = "; ".join( + f"[{i}] viol={m['viol']:.3e} op_viol=" + f"{'n/a' if m['op_viol'] is None else format(m['op_viol'], '.3e')} " + f"accepted={m['accepted']}" + f"{'' if m['accepted'] else ' (' + str(m['rejected_by']) + ')'}" + for i, m in enumerate(gauge.get("modes", []))) return (f"converged={info.get('converged')} " f"ksp_reason={info.get('ksp_reason')} " f"ksp_its={info.get('ksp_its')} " f"newton_its={info.get('nonlinear_iterations')} " f"|r|={info.get('rnorm')} |r0|={info.get('rnorm0')} " f"rotation_reused={info.get('rotation_reused')} " - f"workspace_reused={info.get('workspace_reused')}") + f"workspace_reused={info.get('workspace_reused')} " + f"constrained_rows={None if rows is None else len(rows)} " + f"distinct_rows={None if rows is None else len(set(rows))} " + f"boundaries={info.get('boundaries')} " + f"gauge_offered={gauge.get('offered')} " + f"gauge_orthonormal={gauge.get('orthonormal')} " + f"gauge_removed={info.get('rotation_gauge_removed')} " + f"modes=({modes})") + + +class _LocatorTally: + """Count point-location work, how much came back ``-1``, and how much of + that the #551 rejection radius is responsible for. + + Every call is answered twice — once as shipped, once with + ``_local_cell_reach`` set aside so the walk cannot reject anything early — + and the cells are compared. ``radius_changed=0`` means the rejection radius + did not alter a single answer, which is the measurement that decides + whether the locator is implicated in a downstream disagreement. It is not + an argument, it is a count, and it is taken on whatever mesh the machine + running the test actually built. + """ + + def __init__(self): + self.calls = 0 + self.points = 0 + self.rejected = 0 + self.radius_changed = 0 + self._mesh_cls = uw.discretisation.discretisation_mesh.Mesh + self._wrapped = self._mesh_cls._get_closest_local_cells_internal + + def __enter__(self): + outer = self + + def counted(mesh_self, coords, **kwargs): + got = outer._wrapped(mesh_self, coords, **kwargs) + saved = getattr(mesh_self, "_local_cell_reach", None) + try: + mesh_self._local_cell_reach = None # no early rejection + reference = outer._wrapped(mesh_self, coords, **kwargs) + finally: + mesh_self._local_cell_reach = saved + got_a = np.asarray(got).reshape(-1) + outer.calls += 1 + outer.points += len(coords) + outer.rejected += int(np.count_nonzero(got_a < 0)) + outer.radius_changed += int(np.count_nonzero( + got_a != np.asarray(reference).reshape(-1))) + return got + + self._mesh_cls._get_closest_local_cells_internal = counted + return self + + def __exit__(self, *exc): + self._mesh_cls._get_closest_local_cells_internal = self._wrapped + return False + + def __str__(self): + return (f"calls={self.calls} points={self.points} " + f"returned_-1={self.rejected} " + f"radius_changed={self.radius_changed}") + + +def _rigid_body_decomposition(coords, diff): + """Split a velocity difference field into rigid-body modes and the rest. + + Two solves that both drive the residual to machine zero on a system with + the same initial residual can only differ in the null space, so the + question "is the difference a rigid-body mode?" has a yes/no answer. The + modes are the ones the solver itself considers (``_rigid_rotation_modes``: + ``(-y, x)`` in 2-D, ``e_k x r`` in 3-D) plus the translations, built here + from the nodal coordinates so this is independent of the solver's own + machinery. Nodal (not PETSc-global) inner products — exact in serial, + which is where this test runs. + + Returns ``(|d|, per-direction shares, |d| off the rigid-body span)``. The + modes are NOT mutually orthogonal — the rotation about the origin has a + large translation component on a box in the first quadrant — so they are + Gram-Schmidt'd in order, exactly as the solver does before projecting, and + the per-direction shares are therefore shares on the k-th ORTHOGONALISED + direction, not on the named mode. The number that answers the question is + the off-span norm: a pure rigid rotation gives 4e-18 of it, a random field + gives 0.9985 of it (both measured). + """ + coords = np.asarray(coords, dtype=np.float64) + diff = np.asarray(diff, dtype=np.float64) + dim = coords.shape[1] + modes = {} + for axis in range(dim): + e = np.zeros_like(coords) + e[:, axis] = 1.0 + modes[f"translation_{'xyz'[axis]}"] = e + if dim == 2: + modes["rotation_z"] = np.column_stack([-coords[:, 1], coords[:, 0]]) + else: + x, y, z = coords[:, 0], coords[:, 1], coords[:, 2] + zero = np.zeros_like(x) + modes["rotation_x"] = np.column_stack([zero, -z, y]) + modes["rotation_y"] = np.column_stack([z, zero, -x]) + modes["rotation_z"] = np.column_stack([-y, x, zero]) + + total = np.linalg.norm(diff) + residual = diff.copy() + shares = {} + basis = [] + for name, mode in modes.items(): + w = mode.copy() + for q in basis: # Gram-Schmidt, as the solver does + w -= np.vdot(q, w) * q + n = np.linalg.norm(w) + if n <= 1.0e-14: + continue + w /= n + basis.append(w) + component = np.vdot(w, diff) + shares[name] = abs(component) / (total + 1.0e-300) + residual -= component * w + return total, shares, np.linalg.norm(residual) def _rampable_rotated_stokes(mesh, k_expr, tag, forcing=None): @@ -279,10 +405,13 @@ def test_rotated_workspace_deform_invalidates(): info_1 = dict(s1._rotated_freeslip_info) reach_before = mesh._local_cell_reach - # bump the top boundary (the free-surface pattern) + # bump the top boundary (the free-surface pattern). The deform is where the + # point locator is actually exercised — measured, the solves themselves make + # no location calls at all — so the tally goes here. coords = mesh.X.coords.copy() coords[:, 1] += 0.02 * coords[:, 1] * np.sin(np.pi * coords[:, 0]) - mesh.deform(coords) + with _LocatorTally() as tally_deform: + mesh.deform(coords) # Mesh-side invariant: the kd-tree index and the point locator's rejection # radius (#551) are measured together in _build_kd_tree_index, which @@ -296,32 +425,94 @@ def test_rotated_workspace_deform_invalidates(): f"({reach_before:.8g} -> {reach_after:.8g}) — it is measured with the " f"kd-tree index and must be rebuilt with it") - s1.solve() + # NEGATIVE CONTROL for that tally: squeeze the reach and the comparison + # must see answers change. Otherwise a "radius_changed=0" report is the + # instrument failing to fire, not the radius being inert. + probe = np.ascontiguousarray( + np.random.default_rng(5).uniform(0.02, 0.98, size=(400, 2))) + with _LocatorTally() as tally_tight: + mesh._LOCATOR_REACH_MARGIN = 0.05 + try: + mesh._get_closest_local_cells_internal( + probe, tol=mesh._EVAL_FACE_TOL) + finally: + del mesh._LOCATOR_REACH_MARGIN + assert tally_tight.radius_changed > 0, ( + "a reach margin of 0.05 changed no locator answer, so the " + "radius_changed counts reported below cannot tell an inert rejection " + "radius from an instrument that is not looking") + + with _LocatorTally() as tally_deformed: + s1.solve() info = s1._rotated_freeslip_info assert not info["rotation_reused"], ( "workspace survived a mesh.deform — stale rotation Q in use") assert not info["workspace_reused"] - assert info["converged"], ( - f"the post-deform solve did not converge: {_solve_report(info)}") k_c = uw.function.expression(r"k_dc", 1.0, "control viscosity") - s_c, v_c = _rampable_rotated_stokes(mesh, k_c, "DfC") - s_c.solve() + with _LocatorTally() as tally_control: + s_c, v_c = _rampable_rotated_stokes(mesh, k_c, "DfC") + s_c.solve() info_c = s_c._rotated_freeslip_info + + # Build every diagnostic EAGERLY, not inside the assertion messages: an + # instrument that only runs when the test fails is an instrument that has + # never been run. + report_1 = _solve_report(info_1) + report_deformed = _solve_report(info) + report_control = _solve_report(info_c) + diff = np.asarray(v1.data) - np.asarray(v_c.data) + err = np.linalg.norm(diff) / np.linalg.norm(v_c.data) + total, shares, off_mode = _rigid_body_decomposition(v1.coords, diff) + share_text = " ".join(f"{n}={f:.4f}" for n, f in shares.items()) + off_fraction = off_mode / (total + 1.0e-300) + rigid_fraction = np.sqrt(max(0.0, 1.0 - off_fraction ** 2)) + + # NEGATIVE CONTROL for the decomposition, run on the mesh CI actually + # built: a pure rigid rotation must come out entirely inside the span, a + # random field almost entirely outside it. Without this the off-span + # fraction reported above would be a number nobody had checked. + node_coords = np.asarray(v1.coords, dtype=np.float64) + pure = np.column_stack([-node_coords[:, 1], node_coords[:, 0]]) + _, _, pure_off = _rigid_body_decomposition(v1.coords, pure) + assert pure_off / np.linalg.norm(pure) < 1e-10, ( + f"the rigid-body decomposition does not recognise a pure rotation " + f"({pure_off / np.linalg.norm(pure):.3e} of it off the span)") + noise = np.random.default_rng(0).normal(size=pure.shape) + _, _, noise_off = _rigid_body_decomposition(v1.coords, noise) + assert noise_off / np.linalg.norm(noise) > 0.9, ( + f"the rigid-body decomposition absorbs a random field " + f"({noise_off / np.linalg.norm(noise):.3f} of it off the span), so a " + f"small off-span fraction above would prove nothing") + + assert info["converged"], ( + f"the post-deform solve did not converge: {report_deformed}") assert info_c["converged"], ( - f"the fresh control solve did not converge: {_solve_report(info_c)}") + f"the fresh control solve did not converge: {report_control}") - err = np.linalg.norm(v1.data - v_c.data) / np.linalg.norm(v_c.data) assert err < 1e-6, ( f"post-deform solve differs from fresh control by {err:.2e}\n" - f" first solve : {_solve_report(info_1)}\n" - f" post-deform : {_solve_report(info)}\n" - f" fresh control : {_solve_report(info_c)}\n" + f" first solve : {report_1}\n" + f" post-deform : {report_deformed}\n" + f" fresh control : {report_control}\n" f" locator reach : {reach_before:.8g} -> {reach_after:.8g}\n" - f"Both solves assemble the same system from a zero guess, so equal " - f"convergence reports with a non-zero err means the two OPERATORS " - f"differ — something survived the deform. Different |r| or iteration " - f"counts means the linear solves themselves diverged.") + f" locator work : deform {tally_deform} | post-deform solve " + f"{tally_deformed} | control {tally_control}\n" + f" difference : |d|={total:.6e}\n" + f" IN the rigid-body span: {rigid_fraction:.6f} of |d|\n" + f" OFF it: {off_fraction:.6f} of |d| " + f"({off_mode:.6e})\n" + f" per orthogonalised direction (NOT per named mode, " + f"they are not orthogonal): {share_text}\n" + f"Both solves assemble the same system from a zero guess and both " + f"converge, so a difference can only live in the NULL SPACE. If the " + f"rigid-body span accounts for it (off-fraction near 0; a random field " + f"gives 0.9985), read the per-mode gauge decisions above — one solve " + f"admitted a mode the other rejected, and that is a #543 gauge bug " + f"this branch merely exposes. If the difference is OFF the span, the " + f"two OPERATORS differ and something survived the deform. If the " + f"constrained row counts or the locator tallies differ between the two " + f"solvers, point location is implicated and the radius is ours.") @pytest.mark.level_2 From 1b3efa5377cc0a40a99a780cfee7902b38393592 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 14:26:39 +1000 Subject: [PATCH 9/9] Compare the deform-invalidation solves only where the system determines them The test made two claims in one assertion. One is well posed and is what #543 wrote it for: the deform invalidated the workspace. The other is not: that a post-deform solve matches a fresh control to 1e-6, on a system that does not determine its own answer in one direction (#560). They are now separate. (a) stays hard: rotation_reused and workspace_reused both False, the locator's rejection radius followed the deform, both solves converged. (b) is narrowed, not relaxed. Rotated free-slip on a curved boundary loses the constant-pressure gauge, and the solution acquires a component along one unpinned direction whose amplitude is round-off. Measured: a coordinate change of two machine epsilons (4.44e-16) moves the velocity by 1.33e-01, and the move does not scale with the perturbation - 4.4e-16, 2.2e-15, 1e-14, 1e-12 and 1e-9 all give between 4e-2 and 2e-1. That is why this assertion was green on macOS (err exactly 0.0 in 81 consecutive runs across two PETSc toolchains and nine PYTHONHASHSEEDs) and intermittently red on CI: it passes only where the two assemblies agree bitwise. The unpinned subspace is exactly one-dimensional - five different perturbations move the answer along the same direction to cosine 1.000000, the normalised difference set has singular values [2.236, 4.4e-9, 3.4e-9, 2.3e-9, 1.7e-9], and removing the leading direction leaves 2e-9 of each difference. So the test measures that direction with one extra perturbed solve and requires the two solutions to agree in every OTHER direction, at the same 1e-6 it always used. The tolerance is untouched; it is the claim that is made honest. Two things keep it from becoming a rubber stamp. A negative control injects a 1e-3 discrepancy orthogonal to the unpinned direction and asserts it survives the projection (measured 1.000e-03 against a 1e-4 floor), so the projection cannot absorb a real disagreement. And the test branches on what it measures: if a 2-eps perturbation stops moving the answer - i.e. when #560 is fixed - the projection becomes a no-op and the solutions are compared directly again, with the branch reported in the failure message. Fixing #560 strengthens this test instead of breaking it. All the instrumentation stays: the gauge decisions, the constrained-row counts, the locator tallies and the rigid-body decomposition are what turned an unreadable CI failure into a filed defect, and they are the diagnostic for the next one. Verified: this test 10/10 in amr-dev and 10/10 in dev, the CI batch shape (tests/test_101*py tests/test_102*py) 137 passed in dev, full level_1/tier_a gate 627 passed 0 failed. Underworld development team with AI support from Claude Code --- tests/test_1018_rotated_freeslip.py | 155 +++++++++++++++++++++++----- 1 file changed, 130 insertions(+), 25 deletions(-) diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index 92c91e80..6aa789ab 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -13,6 +13,8 @@ pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] +_EPS = np.finfo(float).eps + def _wrap(dm, m0): return uw.discretisation.Mesh( @@ -247,6 +249,23 @@ def __str__(self): f"radius_changed={self.radius_changed}") +def _deformed_box_rotated_solve(tag, coord_scale=1.0): + """A rotated free-slip solve on the deformed box, optionally with every + coordinate multiplied by ``coord_scale``. Used to MEASURE the direction the + system does not determine (#560): the physical response to a coordinate + change of a few machine epsilons is of that order, so anything the solution + does beyond that is the unpinned mode.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(10, 10), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + c = mesh.X.coords.copy() + c[:, 1] += 0.02 * c[:, 1] * np.sin(np.pi * c[:, 0]) + mesh.deform(c * coord_scale) + k = uw.function.expression(r"k_probe", 1.0, "viscosity") + s, v = _rampable_rotated_stokes(mesh, k, tag) + s.solve() + return np.asarray(v.data).copy() + + def _rigid_body_decomposition(coords, diff): """Split a velocity difference field into rigid-body modes and the rest. @@ -381,20 +400,39 @@ def test_rotated_workspace_constant_ramp_invalidates(): def test_rotated_workspace_deform_invalidates(): - """mesh.deform between solves: geometry changed, so the whole workspace - must be rebuilt (rotation_reused False) and the answer must match a fresh - solver on the deformed mesh. - - Both solvers assemble the same system on the same deformed mesh from the - same (zero) initial guess, so any deterministic linear solver has to give - them the same answer: a non-zero ``err`` means one of the two systems is - not what it should be, not that a Krylov path drifted. The checks below - say WHICH, because this test has failed in CI on a platform where it - cannot be reproduced locally (macOS arm64, both the conda-PETSc `dev` and - the AMR-PETSc `amr-dev` toolchains, single test / whole file / CI batch: - err is exactly 0.0, one Krylov iteration, |r| = 8.3e-12 identical between - the two solves). Reporting the solver's own diagnostics and the mesh's - geometry invariants turns the next occurrence into an attribution.""" + """mesh.deform between solves: the geometry changed, so the whole workspace + must be rebuilt and the answer must match a fresh solver on the deformed + mesh — in every direction the deformed system actually determines. + + Two claims, deliberately separated, because only one of them is well posed: + + **(a) the workspace was invalidated.** ``rotation_reused`` and + ``workspace_reused`` both False, the locator's rejection radius followed + the deform, and both solves converged. This is the contract #543 wrote the + test for and it is asserted hard. + + **(b) the two answers agree.** Rotated free-slip on a CURVED boundary loses + the constant-pressure gauge (#560), leaving one direction the operator does + not pin: a coordinate change of two machine epsilons moves the velocity by + 1.3e-01, and the move does not scale with the perturbation. So "the two + solves agree" is false as stated — a plain comparison passes only where the + two assemblies happen to agree bitwise, which is why this test was green on + macOS (err exactly 0.0 in 81 consecutive runs across two PETSc toolchains + and nine PYTHONHASHSEEDs) and intermittently red on CI (err 7e-2 to 1.2e-1, + the size of the unpinned component rather than a drift). + + The unpinned subspace is one-dimensional (five perturbations move the + answer along the same direction to cosine 1.000000), so the claim is + narrowed rather than dropped: measure that direction with one extra + perturbed solve and require agreement in every other direction. The + tolerance is NOT relaxed — 1e-6, as before — it is the claim that is made + honest. When #560 is fixed the probe stops finding a direction, the + projection becomes a no-op, and the test compares the solutions directly + again with no further edit. + + The instrumentation below (gauge decisions, constrained-row counts, locator + tallies, rigid-body decomposition) stays: it is what turned an unreadable + CI failure into #560, and it is the diagnostic for the next one.""" mesh = uw.meshing.StructuredQuadBox( elementRes=(10, 10), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) @@ -490,8 +528,78 @@ def test_rotated_workspace_deform_invalidates(): assert info_c["converged"], ( f"the fresh control solve did not converge: {report_control}") - assert err < 1e-6, ( - f"post-deform solve differs from fresh control by {err:.2e}\n" + # ------------------------------------------------------------------ + # (b) the answers agree, in the directions the system actually determines + # ------------------------------------------------------------------ + # Rotated free-slip on a CURVED boundary loses the constant-pressure gauge + # (#560): the pressure level runs to ~1e4 against a pressure variation of + # 6.6e-2, and the solution acquires a component along one unpinned + # direction whose amplitude is set by round-off. Measured: a coordinate + # change of two machine epsilons (4.44e-16) moves the velocity by 1.33e-01, + # and the size of the move does not track the size of the perturbation + # (4.4e-16, 2.2e-15, 1e-14, 1e-12 and 1e-9 all give 4e-2 to 2e-1). The same + # solve on a STRAIGHT-walled box moves by 5.9e-15, and native essential + # free-slip on this same deformed mesh moves by 4.8e-13 — so it is the + # rotated path on a curved boundary, and it is present at the merge base. + # + # So "the two solves agree" is not a property this system has, and a plain + # comparison passes only where the two assemblies happen to agree bitwise + # (macOS/arm64: err is exactly 0.0 in 81 consecutive runs across two PETSc + # toolchains and nine PYTHONHASHSEEDs; CI's Linux build: err ~7e-2 to + # 1.2e-1, which is the SIZE of the unpinned component, not a drift). + # + # The unpinned subspace is exactly ONE-dimensional — five different + # perturbations move the answer along the same direction to cosine + # 1.000000, the normalised difference set has singular values + # [2.236, 4.4e-9, ...], and removing the leading direction leaves 2e-9 of + # each difference. So the comparison can be made well posed rather than + # abandoned: measure that direction with one extra perturbed solve and + # assert the two solutions agree in every OTHER direction. + probe = _deformed_box_rotated_solve("Prb", coord_scale=1.0 + 2 * _EPS) + unpinned = np.asarray(probe) - np.asarray(v_c.data) + unpinned_size = np.linalg.norm(unpinned) / np.linalg.norm(v_c.data) + + if unpinned_size > 1.0e-3: + # #560 is present (the expected branch today). Project it out. + direction = (unpinned / np.linalg.norm(unpinned)).ravel() + flat = diff.ravel() + residual = flat - float(np.dot(flat, direction)) * direction + constrained_err = np.linalg.norm(residual) / np.linalg.norm(v_c.data) + branch = (f"#560 present: the unpinned direction carries " + f"{unpinned_size:.3e} of the solution, projected out") + else: + # #560 has been fixed — there is no unpinned direction to remove, so + # compare directly and let this test go back to its full strength. + constrained_err = err + branch = (f"#560 appears FIXED (a 2-eps perturbation moves the answer " + f"by only {unpinned_size:.3e}) — the projection below is " + f"now a no-op and this test is comparing solutions directly. " + f"Delete the projection and the _deformed_box_rotated_solve " + f"probe.") + + # NEGATIVE CONTROL: the projection must not absorb a genuine discrepancy. + # Inject a difference orthogonal to the unpinned direction and check it + # survives, or "constrained_err is small" would be true of anything. + injected = np.random.default_rng(6).normal(size=diff.shape).ravel() + if unpinned_size > 1.0e-3: + d_hat = (unpinned / np.linalg.norm(unpinned)).ravel() + injected -= float(np.dot(injected, d_hat)) * d_hat + injected *= 1.0e-3 * np.linalg.norm(v_c.data) / np.linalg.norm(injected) + poisoned = diff.ravel() + injected + if unpinned_size > 1.0e-3: + d_hat = (unpinned / np.linalg.norm(unpinned)).ravel() + poisoned = poisoned - float(np.dot(poisoned, d_hat)) * d_hat + poisoned_err = np.linalg.norm(poisoned) / np.linalg.norm(v_c.data) + assert poisoned_err > 1.0e-4, ( + f"a deliberate 1e-3 discrepancy orthogonal to the unpinned direction " + f"survives the projection as only {poisoned_err:.3e}, so the " + f"constrained comparison below would not notice a real disagreement") + + assert constrained_err < 1e-6, ( + f"post-deform solve differs from fresh control by {constrained_err:.2e} " + f"OUTSIDE the direction the system leaves undetermined " + f"(raw difference {err:.2e})\n" + f" branch : {branch}\n" f" first solve : {report_1}\n" f" post-deform : {report_deformed}\n" f" fresh control : {report_control}\n" @@ -504,15 +612,12 @@ def test_rotated_workspace_deform_invalidates(): f"({off_mode:.6e})\n" f" per orthogonalised direction (NOT per named mode, " f"they are not orthogonal): {share_text}\n" - f"Both solves assemble the same system from a zero guess and both " - f"converge, so a difference can only live in the NULL SPACE. If the " - f"rigid-body span accounts for it (off-fraction near 0; a random field " - f"gives 0.9985), read the per-mode gauge decisions above — one solve " - f"admitted a mode the other rejected, and that is a #543 gauge bug " - f"this branch merely exposes. If the difference is OFF the span, the " - f"two OPERATORS differ and something survived the deform. If the " - f"constrained row counts or the locator tallies differ between the two " - f"solvers, point location is implicated and the radius is ours.") + f"This is the assertion that survives #560: the two solves must agree " + f"in every direction the operator determines. A failure here is NOT " + f"the pressure gauge. If the constrained row counts or the locator " + f"tallies differ between the two solvers, point location is " + f"implicated; if they match, the two OPERATORS differ and something " + f"survived the deform.") @pytest.mark.level_2