diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index f141d859..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 @@ -5288,6 +5265,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 +5280,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,60 +5345,15 @@ def _build_kd_tree_index(self): centroids_list, dtype=numpy.float64).reshape(-1, self.cdim) self._centroid_index = uw.kdtree.KDTree(self._nav_centroids) - 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() + # Rejection radius for the lost-point walk. Rebuilt with the kd-tree + # (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 @@ -5867,6 +5807,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 @@ -5881,7 +5853,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 @@ -5889,7 +5864,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 @@ -5917,11 +5893,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 @@ -5931,11 +5913,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: @@ -5984,6 +5970,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 +5987,34 @@ 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. + + .. 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 @@ -6046,7 +6064,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 +6091,42 @@ 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. + # + # 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 + 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 +6136,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 @@ -6240,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 ee13cb97..6c785af0 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 @@ -1250,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 @@ -1286,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" @@ -1311,24 +1354,44 @@ 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. + # + # 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: - 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) + if cell_hints is not None and mesh is hinted_mesh: + # 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 " + 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.get_closest_cells(coords) + cells = mesh._robust_owning_cells(coords) else: cells = None 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_0761_point_locator.py b/tests/test_0761_point_locator.py new file mode 100644 index 00000000..874edd0c --- /dev/null +++ b/tests/test_0761_point_locator.py @@ -0,0 +1,673 @@ +"""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). + +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. +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. + +**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 +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 _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) + 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") + + +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) +# --------------------------------------------------------------------------- + +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() + + # 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 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() diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index c7fba4b7..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( @@ -165,6 +167,162 @@ 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, 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"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 _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. + + 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): """Rotated free-slip Stokes with viscosity given by a rampable UWexpression constant (the #416 idiom used by every continuation @@ -242,9 +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.""" + """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) @@ -252,24 +440,184 @@ 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) + # 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) - - s1.solve() + 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 + # 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") + + # 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"] 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() - 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}") + 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: {report_control}") + + # ------------------------------------------------------------------ + # (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" + f" locator reach : {reach_before:.8g} -> {reach_after:.8g}\n" + 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"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