From 729123fce4b803ee90c03d71453b2cae7a866167 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sat, 24 Oct 2026 19:17:07 +1100 Subject: [PATCH] Preserve physical query positions near exact-cell faces Tighten full-dimensional exact-cell containment so neighbouring-cell hints cannot shift continuous-field probes through reference-coordinate clamping. Keep approximate warped-face and embedded-manifold tolerances unchanged. Validation: Gadi int32 installation at feature ef62a45b, 22 native regressions in serial and on each of 16 ranks, 272 annulus components, and exact submesh-transfer gates passed. Prior-runtime serial/MPI reproducer has 5.806e-6 error; tight tolerance reduces it below 8e-17. Four matched annulus workflows are being independently audited. --- .../discretisation/discretisation_mesh.py | 38 +++++++++----- .../test_1090_evaluation_near_interface.py | 50 +++++++++++++++++++ 2 files changed, 76 insertions(+), 12 deletions(-) create mode 100644 tests/parallel/test_1090_evaluation_near_interface.py diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index ccfca8c36..71cc80e0b 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -6157,14 +6157,11 @@ def _test_if_points_in_cells_internal(self, points, cells, roughly ``tol`` of the face (relative to the face-normal separation), while still rejecting points that lie inside a *different* cell (whose diff is strongly negative). The - parallel evaluation locator - (``Mesh._robust_owning_cells``) uses ``tol=1e-2`` to admit - on-face / partition-seam node queries that ``on_boundary=True``'s - absolute 1e-12 floor is too tight to accept — for query - coordinates slightly off the face (e.g. RBF-shifted node - points), the geometric scale of the test needs to match the - *mesh* spacing, not float roundoff. See memory - project_pr207_loose_boundary_clash. + evaluation locator uses a tight tolerance for full-dimensional, + geometrically exact cells. Its previous ``1e-2`` tolerance could + admit a point genuinely inside a neighbour; reference-coordinate + clamping then moved the query. Approximate warped-face and + embedded-manifold evaluation retain their previous tolerance. """ # Internal version - points assumed to already be in model units self._mark_faces_inside_and_out() @@ -6855,14 +6852,26 @@ def get_closest_local_cells( # point another rank owns. A point in the (closed) domain is then found by # exactly one rank (verified: 360/360, owners==1), and the migration routes # it to that owner. - _EVAL_FACE_TOL = 1.0e-2 + _EVAL_FACE_TOL = 1.0e-10 + _EVAL_APPROX_FACE_TOL = 1.0e-2 def _robust_owning_cells(self, coords: numpy.ndarray) -> numpy.ndarray: """Per-point owning cell for parallel evaluation (coords in model units). This is the strict barycentric/cell-wall locator - (:meth:`_get_closest_local_cells_internal`) with a *tight* face - tolerance (``_EVAL_FACE_TOL``): it returns the containing cell for + (:meth:`_get_closest_local_cells_internal`) with a roundoff-scale face + tolerance for full-dimensional, geometrically exact cells + (``_EVAL_FACE_TOL``). The previous ``1e-2`` tolerance relative to the + face control-point separation could claim a point in its neighbour; + reference-coordinate clamping then moves that point onto a face. Continuity alone does not + justify this: an exactly represented continuous field with a gradient + jump can acquire finite point errors. Tight containment returns the + actual containing cell, including near an internal density sheet. + + Approximate warped faces and embedded manifolds retain their previous + tolerance (``_EVAL_APPROX_FACE_TOL``); this change does not extend the + cell-wall estimator's geometric authority. The exact-cell locator returns + the containing cell for interior points and a valid sharing cell for genuinely-on-face points, and ``-1`` for points that lie inside a *different* cell or outside the local mesh. Crucially it does **not** fall back to a bounding-sphere @@ -6884,8 +6893,13 @@ def _robust_owning_cells(self, coords: numpy.ndarray) -> numpy.ndarray: coords = numpy.asarray(coords) if coords.shape[0] == 0: return numpy.array([], dtype=numpy.int64) + tolerance = ( + self._EVAL_FACE_TOL + if self.dim == self.cdim and self._location_capability() == "exact" + else self._EVAL_APPROX_FACE_TOL + ) return numpy.asarray( - self._get_closest_local_cells_internal(coords, tol=self._EVAL_FACE_TOL), + self._get_closest_local_cells_internal(coords, tol=tolerance), dtype=numpy.int64, ) diff --git a/tests/parallel/test_1090_evaluation_near_interface.py b/tests/parallel/test_1090_evaluation_near_interface.py new file mode 100644 index 000000000..02f839cb7 --- /dev/null +++ b/tests/parallel/test_1090_evaluation_near_interface.py @@ -0,0 +1,50 @@ +"""Exact-cell interpolation must preserve physical points near shared faces. + +A wide ownership tolerance can select a neighbouring cell, after which the +reference-coordinate clamp projects a near-face query onto that face. Both a +continuous kink and a global polynomial are exactly represented here, so their +point errors must be roundoff rather than a finite positional error. +""" + +import numpy as np +import pytest +import underworld3 as uw + + +def near_interface_geometry(): + mesh = uw.meshing.BoxInternalBoundary( + elementRes=(32,32), zelementRes=(16,16), zintCoord=.5, + minCoords=(0.,0.), maxCoords=(1.,1.), simplex=True, qdegree=2, + ) + x = np.random.default_rng(19).uniform(.02,.98,500) + points = np.concatenate([ + np.column_stack((x,np.full(len(x),.5+delta))) + for delta in (-1e-4,-1e-6,0.,1e-6,1e-4) + ]) + return mesh, points + + +@pytest.mark.parametrize("degree", (1,2)) +def test_continuous_gradient_jump_is_evaluated_in_the_containing_cell(degree): + mesh, points = near_interface_geometry() + field = uw.discretisation.MeshVariable("Kink",mesh,1,degree=degree) + field.data[:,0] = np.abs(np.asarray(field.coords)[:,1]-.5) + values, outside = uw.function.global_evaluate(field.sym[0],points,check_extrapolated=True) + assert not np.any(outside) + np.testing.assert_allclose(np.asarray(values).reshape(-1),np.abs(points[:,1]-.5),rtol=0,atol=1e-12) + + +@pytest.mark.parametrize("degree", (1,2)) +def test_polynomial_queries_are_not_projected_onto_nearby_faces(degree): + mesh, points = near_interface_geometry() + field = uw.discretisation.MeshVariable("Polynomial",mesh,1,degree=degree) + nodes = np.asarray(field.coords) + nodal = 1+nodes[:,0]+2*nodes[:,1] + expected = 1+points[:,0]+2*points[:,1] + if degree == 2: + nodal += nodes[:,0]**2-nodes[:,0]*nodes[:,1] + expected += points[:,0]**2-points[:,0]*points[:,1] + field.data[:,0] = nodal + values, outside = uw.function.global_evaluate(field.sym[0],points,check_extrapolated=True) + assert not np.any(outside) + np.testing.assert_allclose(np.asarray(values).reshape(-1),expected,rtol=0,atol=1e-12)