diff --git a/.github/workflows/build_uw3_and_test.yaml b/.github/workflows/build_uw3_and_test.yaml index 3918c0f3e..9fa0e8a23 100644 --- a/.github/workflows/build_uw3_and_test.yaml +++ b/.github/workflows/build_uw3_and_test.yaml @@ -44,5 +44,16 @@ jobs: # editable installs are project policy violation, see CLAUDE.md). run: pixi run -e dev build + # Import smoke test: a hard crash during `import underworld3` (segfault + # / abort in a compiled module, e.g. from a poisoned env cache) kills + # pytest BEFORE its first stdout flush, so every batch below fails with + # zero output and the run is undiagnosable from the log (this happened + # on PR #394 — a runner-side flake that cost a full debug round-trip). + # This 5-second step makes any import-time crash loud and localised: + # faulthandler prints the native traceback here rather than nothing + # there. + - name: Import smoke test (faulthandler) + run: pixi run -e dev python -u -X faulthandler -c "import faulthandler; faulthandler.enable(); import underworld3; print('import OK', underworld3.__version__)" + - name: Run tests run: pixi run -e dev ./scripts/test.sh --p 2 diff --git a/docs/developer/design/point-location-capability.md b/docs/developer/design/point-location-capability.md new file mode 100644 index 000000000..c03e722af --- /dev/null +++ b/docs/developer/design/point-location-capability.md @@ -0,0 +1,71 @@ +# Point-location capability and the evaluation fallback ladder + +*2026-07-23 — design record for issues #390 / #392 (PRs #391 and the +capability-ladder follow-up). The full discussion, probe numbers, and +acceptance criteria live on issue #392.* + +## The problem + +PETSc's `DMLocatePoints` uses a half-open cell convention: a query point +sitting exactly on the domain's closed upper boundary faces belongs to the +cell "above", which does not exist, so the point is silently dropped. The +evaluator used to zero-fill dropped points — a plausible-looking field value +that propagates silently. Semi-Lagrangian stress-history trace-backs slide +departure points along boundaries, so on quad meshes every step read a +corrupted history there (the VEP stability blow-up of issue #390). + +Underworld's own cell-wall estimator (`_test_if_points_in_cells_internal`, +a half-space intersection over face control-point planes) handles on-face +points correctly — but it is only *authoritative* on some geometries. + +## When the estimator is authoritative + +Its authority is governed by exactly two measurable properties: **face +planarity** (sagitta — max vertex deviation from the face's plane, relative +to face diameter) and **cell convexity** (no cell vertex outside another +face's plane). Measured regimes: + +| Capability | Meshes | Authority | +|------------|--------|-----------| +| `exact` | simplex; manifold (dim ≠ cdim); **any valid 2-D quad mesh, deformed included** (straight edges are always planar, convexity ⇔ untangled); planar-faced hexes (rectilinear boxes, affine images) | all field types — bypass `DMLocatePoints` | +| `continuous` | warped hexes within sagitta ≤ 5e-2 (cubed-sphere class: lateral faces exactly planar, spherical faces ≤ 1e-2) | continuous fields only — a misclassified point lands in the face-adjacent neighbour where continuous FE interpolants agree to O(sagitta × gradient); a face-aligned jump (e.g. layered viscosity on a spherical shell) would see O(jump) errors | +| `none` | badly warped or non-convex cells | none — PETSc locates; dropped points take the RBF fallback | + +The capability is a **measured mesh property**, not a type flag: +`Mesh._location_capability()` runs the geometry audit +(`_audit_cell_face_geometry`) and caches against +`(_mesh_version, _topology_version)`, so `deform()` and adaptation refresh +it automatically. It is per-rank (the evaluator runs on `COMM_SELF`; a +collective reduction here could deadlock since `petsc_interpolate` is +reached only by ranks holding points). + +## The fallback ladder + +1. **Authoritative** (`exact`, or `continuous` with all-continuous fields): + the estimator's owner is passed as an authoritative hint and + `DMLocatePoints` is bypassed (`petsc_tools.c`), with reference-coordinate + clamping for on-face queries. +2. **Not authoritative**: `DMLocatePoints` decides. Points it drops get + **NaN** in the interpolant (never zeros) and are reported in + `unlocated_mask`; the caller fills them per-variable via + `rbf_interpolate` — bounded, topology-free, honest. NaN survives only if + that plumbing is bypassed, which is exactly when it should be visible. +3. Genuinely out-of-domain points never reach this machinery — they are + classified out by `points_in_domain` and take the existing RBF + extrapolation path. + +The policy decision lives in ONE place — +`Mesh._hint_is_authoritative(all_fields_continuous)` — consumed by +`petsc_interpolate` (which also folds the policy into the DMInterpolation +cache key, since the same coords with a different field-continuity mix can +need both structures). + +## Sentinel tests + +- `tests/test_0503_evaluate.py::test_evaluate_on_domain_boundary_faces` — + boundary-face evaluation exact on quad and simplex boxes (the #390 class). +- `test_location_capability_measured` — capability grading of regular / + deformed quad and hex meshes. +- `test_evaluate_deformed_quad_boundary_exact` — the 2-D graduation. +- `test_evaluate_warped_hex_rbf_fallback_no_silent_values` — the RBF rung: + no NaN escapes, no silent zeros, values bounded by the data range. diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index f0c02025c..d26630c7b 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -5831,14 +5831,132 @@ def _robust_owning_cells(self, coords: numpy.ndarray) -> numpy.ndarray: dtype=numpy.int64, ) + def _audit_cell_face_geometry(self): + """Measure the two quantities that bound the cell-wall estimator's + authority on quad/hex meshes: worst face non-planarity (sagitta) and + worst outward half-space violation by a cell's own vertices, both + relative to the face diameter. + + The in-cell test (:meth:`_test_if_points_in_cells_internal`) is a + half-space intersection over the face planes the control points + define. It is exact when every face is planar and every cell is + convex; a warped face confines misclassification to a slab of + thickness ~sagitta at that face (and the misassigned cell is the + face-adjacent neighbour, where continuous FE interpolants agree). + The violation term catches non-convex / tangled cells, whose + half-space intersection no longer represents the cell at all. + + Returns ``(max_sagitta_rel, max_violation_rel)`` over local cells; + ``(0.0, 0.0)`` on an empty local mesh. + """ + nav_dm = self._nav_dm if self._nav_dm is not None else self.dm + nav_coords = self._nav_coords + + cStart, cEnd = nav_dm.getHeightStratum(0) + if cEnd == cStart: + return 0.0, 0.0 + cell_num_faces = self.element.entities[1] + cell_num_points = self.element.entities[self.dim] + face_num_points = self.element.face_entities[self.dim] + + max_sag = 0.0 + max_viol = 0.0 + for cell_id in range(cStart, cEnd): + cell_faces = nav_dm.getCone(cell_id) + cpoints = nav_dm.getTransitiveClosure(cell_id)[0][-cell_num_points:] + cell_point_coords = nav_coords[self._coord_rows_for_points(nav_dm, cpoints)] + cell_centroid = cell_point_coords.mean(axis=0) + for face in range(cell_num_faces): + fpoints = nav_dm.getTransitiveClosure(cell_faces[face])[0][-face_num_points:] + point_coords = nav_coords[self._coord_rows_for_points(nav_dm, fpoints)] + face_centroid = point_coords.mean(axis=0) + normal = self._facet_outward_unit_normal( + point_coords, face_centroid, cell_centroid, + cell_point_coords=cell_point_coords, + ) + diam = float(numpy.linalg.norm( + point_coords.max(axis=0) - point_coords.min(axis=0))) + if diam == 0.0: + continue + sag = float(numpy.abs((point_coords - face_centroid) @ normal).max()) + viol = float(((cell_point_coords - face_centroid) @ normal).max()) + max_sag = max(max_sag, sag / diam) + max_viol = max(max_viol, viol / diam) + return max_sag, max_viol + + # Capability thresholds. EXACT demands machine-planar faces and convex + # cells; CONTINUOUS admits face warp up to 5% of the face diameter (the + # cubed sphere measures ~1e-2 on its spherical faces; a smoothly deformed + # hex box ~1e-2). Beyond that the estimator's error bound is no longer + # small against per-cell field variation and it loses authority entirely. + _LOCATION_EXACT_TOL = 1.0e-9 + _LOCATION_CONTINUOUS_TOL = 5.0e-2 + + def _location_capability(self) -> str: + """Measured point-location capability of this mesh's cell-wall + estimator: ``"exact"``, ``"continuous"``, or ``"none"``. + + * ``"exact"`` — the estimator is authoritative for every field type: + simplex meshes, manifold meshes (dim != cdim, where PETSc's own + in-cell test is the unreliable party), and any quad/hex mesh whose + measured face sagitta and convexity violation are at machine level + (rectilinear boxes, affine images, and *all valid 2-D quad meshes* — + straight edges are always planar and convexity is equivalent to an + untangled mesh). + * ``"continuous"`` — authoritative for continuous fields only: warped + hexes within the sagitta tolerance (cubed-sphere class). A point in + the misclassification slab lands in the face-adjacent neighbour, + where continuous FE interpolants agree to O(sagitta x gradient) — + but a field with a face-aligned jump would see O(jump) errors, so + discontinuous evaluation must not trust the estimator here. + * ``"none"`` — badly warped or non-convex cells: the estimator carries + no authority and dropped points take the RBF fallback. + + The audit is computed per rank over local cells and cached against + ``(_mesh_version, _topology_version)`` — ``deform()`` and adaptation + both bump these. It is deliberately NOT reduced across ranks: the + evaluator runs on COMM_SELF, so a (pathological) capability split at + a threshold is rank-local and safe, while a collective reduction here + could deadlock (petsc_interpolate is reached only by ranks holding + points). + """ + key = (self._mesh_version, self._topology_version) + cached = getattr(self, "_location_capability_cache", None) + if cached is not None and cached[0] == key: + return cached[1] + + if bool(self.dm.isSimplex()) or (self.dim != self.cdim): + cap = "exact" + else: + sag_rel, viol_rel = self._audit_cell_face_geometry() + worst = max(sag_rel, viol_rel) + if worst < self._LOCATION_EXACT_TOL: + cap = "exact" + elif worst < self._LOCATION_CONTINUOUS_TOL: + cap = "continuous" + else: + cap = "none" + self._location_capability_cache = (key, cap) + return cap + + def _hint_is_authoritative(self, all_fields_continuous: bool = True) -> bool: + """Whether the cell-wall hint may bypass ``DMLocatePoints`` for an + evaluation involving the given field continuity. See + :meth:`_location_capability` for the regimes; ``"continuous"`` + capability is only authoritative when every interpolated variable is + continuous. + """ + cap = self._location_capability() + return cap == "exact" or (cap == "continuous" and bool(all_fields_continuous)) + def _eval_use_robust_location(self) -> bool: """Single switch for the parallel evaluation cell-location strategy. Returns True when ``uw.function`` evaluation should locate cells with the bulletproof barycentric hint (:meth:`_robust_owning_cells`) and the ``DMLocatePoints`` bypass (``petsc_tools.c``), rather than PETSc's - ``DMLocatePoints``. This is the *one place* the policy lives; the - evaluate_nd classifier, the petsc_interpolate hint, and the + ``DMLocatePoints``. This is the *one place* the parallel policy lives; + the evaluate_nd classifier, the petsc_interpolate hint, and the DMInterpolation wrapper all defer to it so the three stay consistent. Two conditions, both required: @@ -5848,16 +5966,15 @@ def _eval_use_robust_location(self) -> bool: are reliable with a single domain, so serial keeps the validated path bit-for-bit. The parallel-only failure is the rank-local RBF / wrong- region value at partition-seam node points. - * **hint is authoritative** — simplex cells (``dm.isSimplex()``: planar - faces, affine reference map, exact face containment) or manifold - meshes (``dim != cdim``: PETSc's in-cell test is unreliable near - 2-manifold simplex edges in 3-D). On non-simplex volume meshes - (quads/hexes) deformed faces can be non-planar and the kdtree-nearest - cell can be wrong, so those keep PETSc's DMLocatePoints search. - """ - return (uw.mpi.size > 1) and ( - bool(self.dm.isSimplex()) or (self.dim != self.cdim) - ) + * **the estimator carries geometric authority** — measured capability + ``"exact"`` or ``"continuous"`` (:meth:`_location_capability`). For + in/out *classification* the continuous regime is sufficient: the + only ambiguity is within the sagitta slab of the domain boundary, + which is exactly the tolerance the classifier already accepts. + Capability ``"none"`` (badly warped / non-convex quad meshes) keeps + PETSc's DMLocatePoints search. + """ + return (uw.mpi.size > 1) and (self._location_capability() != "none") def _get_mesh_sizes(self, verbose=False): """ diff --git a/src/underworld3/function/_dminterp_wrapper.pyx b/src/underworld3/function/_dminterp_wrapper.pyx index b208fae5d..6c10c7aef 100644 --- a/src/underworld3/function/_dminterp_wrapper.pyx +++ b/src/underworld3/function/_dminterp_wrapper.pyx @@ -36,6 +36,7 @@ cdef extern from "petsc.h" nogil: cdef extern from "petsc_tools.h" nogil: PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ipInfo, PetscDM dm, int petscbool, int petscbool, size_t* owning_cell, int petscbool) PetscErrorCode DMInterpolationEvaluate_UW(DMInterpolationInfo ipInfo, PetscDM dm, PetscVec x, PetscVec v) + PetscErrorCode DMInterpolationGetUnlocated_UW(DMInterpolationInfo ipInfo, PetscInt n, signed char* mask) cdef class CachedDMInterpolationInfo: """ @@ -62,7 +63,8 @@ cdef class CachedDMInterpolationInfo: cdef DMInterpolationInfo _ipInfo cdef public object coords # numpy array - Python keeps alive - cdef public object cells # numpy array - Python keeps alive + cdef public object cells # numpy array - Python keeps alive (may be None) + cdef public object unlocated_mask # bool ndarray: points with no owning cell cdef public int dofcount cdef public int dim cdef public bint is_valid @@ -73,11 +75,12 @@ cdef class CachedDMInterpolationInfo: """Initialize - structure not yet created.""" self.is_valid = False self.use_count = 0 + self.unlocated_mask = None import time self.creation_time = time.time() def create_structure(self, mesh, np.ndarray[double, ndim=2] coords, - np.ndarray[long, ndim=1] cells, int dofcount): + object cells, int dofcount, bint hint_authoritative=False): """ Create and set up the DMInterpolation structure. @@ -89,10 +92,18 @@ cdef class CachedDMInterpolationInfo: The mesh object coords : ndarray (n_points, dim) Coordinates to interpolate at - cells : ndarray (n_points,) - Cell hints for each coordinate + cells : ndarray (n_points,) or None + Cell hints for each coordinate. With ``hint_authoritative=True`` + the hint bypasses ``DMLocatePoints`` entirely; with ``None`` (or + ``hint_authoritative=False``) PETSc's search is authoritative and + unlocated points are reported in :attr:`unlocated_mask` for the + caller's RBF fallback. dofcount : int Total number of DOFs to interpolate + hint_authoritative : bool + Whether the supplied cells carry geometric authority — decided by + the mesh's measured capability + (``mesh._hint_is_authoritative(...)``), the single policy point. """ cdef PetscErrorCode ierr cdef int n_points = coords.shape[0] @@ -100,7 +111,7 @@ cdef class CachedDMInterpolationInfo: # Store references (Python keeps these alive) self.coords = coords.copy() # CRITICAL: keep alive! - self.cells = cells.copy() + self.cells = None if cells is None else np.asarray(cells).copy() self.dofcount = dofcount self.dim = dim @@ -141,31 +152,27 @@ cdef class CachedDMInterpolationInfo: # ignoreOutsideDomain=1: PETSc silently skips points it cannot locate # rather than crashing. # - # Hint policy — the hint is ALWAYS passed; what varies is its authority - # (the trailing hintAuthoritative flag): + # Hint policy — decided by the CALLER from the mesh's measured + # location capability (mesh._hint_is_authoritative, the single policy + # point) and passed in as hint_authoritative: # - # AUTHORITATIVE (bypass): simplex cells (planar faces, affine reference - # map) or manifold meshes (dim != cdim), where PETSc's own in-cell test - # is the unreliable party (pseudo-inverse chord-plane projection) and - # the barycentric hint correctly contains every query coord -> - # petsc_tools.c skips DMLocatePoints entirely (ported from 17a5a8d) - # and evaluates in the hinted cell. Independent of rank count, so it - # applies in serial too — the fast FE trace-back (PR #203). The hint - # *source* still follows mesh._eval_use_robust_location(): parallel -> - # _robust_owning_cells (correct owner across seams); serial -> the - # standard locator's cells. + # AUTHORITATIVE (bypass): the cell-wall estimator carries geometric + # authority — simplex/manifold meshes, and quad/hex meshes whose + # measured face planarity + convexity qualify ("exact", or + # "continuous" capability with all-continuous fields). petsc_tools.c + # skips DMLocatePoints entirely and evaluates in the hinted cell + # (with reference-coord clamping for on-face queries). # - # FALLBACK (non-simplex volume meshes — quad/hex, incl. deformed - # faces): PETSc DMLocatePoints remains authoritative (the - # kdtree-nearest hint can be wrong there), but the hint still prefills - # the recovery map so points DMLocatePoints DROPS — notably queries - # sitting exactly on the domain's closed upper faces, e.g. SL - # trace-back departure points sliding along the top boundary — are - # evaluated in the adjacent hinted cell (with reference-coord clamping) - # instead of silently zero-filled. Losing this recovery was the ψ* - # corruption behind the VEP stability blow-up (#390). - cdef bint hint_authoritative = (bool(mesh.dm.isSimplex()) or (mesh.dim != mesh.cdim)) - if n_points > 0: + # NOT AUTHORITATIVE: PETSc DMLocatePoints runs. Points it drops (e.g. + # queries exactly on the domain's closed upper faces — the ψ* + # corruption behind the VEP blow-up, #390) get NaN in the interpolant + # and are reported in unlocated_mask; the caller fills them via the + # RBF fallback. A supplied non-authoritative hint would prefill the + # C recovery map instead, but the capability policy passes cells=None + # here precisely because that hint is not trustworthy on such meshes. + cdef np.ndarray mask_arr = np.zeros(n_points, dtype=np.int8) + cdef signed char[::1] mask_view = mask_arr + if n_points > 0 and self.cells is not None: cells_view = np.ascontiguousarray(self.cells) ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, &cells_view[0], @@ -176,6 +183,16 @@ cdef class CachedDMInterpolationInfo: DMInterpolationDestroy(&self._ipInfo) raise RuntimeError(f"DMInterpolationSetUp_UW failed with error {ierr}") + # Which points ended up with no owning cell (dropped by + # DMLocatePoints, or hinted -1 by the robust locator)? Their + # interpolant slots will be NaN; the caller fills them via RBF. + if n_points > 0: + ierr = DMInterpolationGetUnlocated_UW(self._ipInfo, n_points, &mask_view[0]) + if ierr != 0: + DMInterpolationDestroy(&self._ipInfo) + raise RuntimeError(f"DMInterpolationGetUnlocated_UW failed with error {ierr}") + self.unlocated_mask = mask_arr.astype(bool) + self.is_valid = True def evaluate(self, mesh, np.ndarray[double, ndim=2] outarray): diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 607cb388c..ee13cb97c 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -1276,8 +1276,23 @@ def petsc_interpolate( expr, # Try to get cached structure first from underworld3.function._dminterp_wrapper import CachedDMInterpolationInfo + # Location policy: is the cell-wall hint authoritative for THIS + # evaluation? Decided by the mesh's measured capability (face + # planarity + convexity, mesh._location_capability) combined with the + # continuity of the variables being interpolated — "continuous" + # capability (small-sagitta warped hexes, e.g. the cubed sphere) is + # authoritative only when every field is continuous, because a + # face-aligned jump inside the misclassification slab would see + # 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) + authoritative = mesh._hint_is_authoritative(all_continuous) + location_policy = "auth" if authoritative else "locate" + # coords is already np.ndarray type in petsc_interpolate function signature - cached_info = mesh._dminterpolation_cache.get_structure(coords, dofcount) + cached_info = mesh._dminterpolation_cache.get_structure( + coords, dofcount, policy=location_policy) # Create output array cdef np.ndarray outarray = np.empty([len(coords), dofcount], dtype=np.float64) @@ -1295,27 +1310,34 @@ def petsc_interpolate( expr, # CACHE MISS - Create structure and cache it cached_info = CachedDMInterpolationInfo() - # Get cell hints. - # In PARALLEL use the bulletproof barycentric locator (the swarm- - # migration locator): get_closest_cells (first-pass kd-tree-nearest) - # can hand back a non-containing cell for on-face/seam node points, - # and that wrong cell is what the DMInterpolation recovery uses when - # DMLocatePoints drops the point -> a value from the wrong region. - # _robust_owning_cells returns the true containing cell (or a valid - # adjacent cell for on-face points). When the bypass is active - # (mesh._eval_use_robust_location()) this hint is trusted directly; - # otherwise it is the first-pass guess as before. Same single policy - # switch as the classifier and the DMInterpolation wrapper. - if mesh._eval_use_robust_location(): - cells = mesh._robust_owning_cells(coords) + # 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. + # 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) else: - cells = mesh.get_closest_cells(coords) + cells = None # Create and set up DMInterpolation structure (EXPENSIVE) # This calls DMLocatePoints which is COLLECTIVE — all ranks must enter. try: # coords is already np.ndarray type (function signature ensures this) - cached_info.create_structure(mesh, coords, cells, dofcount) + cached_info.create_structure(mesh, coords, cells, dofcount, + hint_authoritative=authoritative) except RuntimeError as e: # Handle DMInterpolationSetUp failures gracefully if "outside the domain" in str(e): @@ -1326,13 +1348,29 @@ def petsc_interpolate( expr, # Store in cache for reuse # coords is already np.ndarray type (function signature ensures this) - mesh._dminterpolation_cache.store_structure(coords, dofcount, cached_info) + mesh._dminterpolation_cache.store_structure( + coords, dofcount, cached_info, policy=location_policy) # Evaluate # swarm_sync=False: see the cache-hit branch above — only a # subset of ranks reaches petsc_interpolate. mesh.update_lvec(swarm_sync=False) cached_info.evaluate(mesh, outarray) + + # RBF fallback rung: points no cell owns (dropped by DMLocatePoints + # on a non-authoritative mesh, or hinted -1) hold NaN in outarray. + # Fill them per-variable with the bounded, topology-free RBF + # interpolant — the same machinery exterior points already use. NaN + # survives only if this plumbing is bypassed, which is exactly when + # it should be visible. + unlocated = getattr(cached_info, "unlocated_mask", None) + if unlocated is not None and unlocated.any(): + fallback_coords = coords[unlocated] + for var in vars: + var_start = var_start_index[var] + rbf_vals = np.asarray(var.rbf_interpolate(fallback_coords)) + rbf_vals = rbf_vals.reshape(len(fallback_coords), var.num_components) + outarray[unlocated, var_start:var_start + var.num_components] = rbf_vals # === END CACHING === # Create map between array slices and variable functions diff --git a/src/underworld3/function/dminterpolation_cache.py b/src/underworld3/function/dminterpolation_cache.py index e115b8b0b..fdb63f7c1 100644 --- a/src/underworld3/function/dminterpolation_cache.py +++ b/src/underworld3/function/dminterpolation_cache.py @@ -25,7 +25,7 @@ class DMInterpolationCache: """ Per-mesh cache for DMInterpolation structures. - Cache key: (coord_hash, dofcount) + Cache key: (coord_hash, dofcount, policy) Cache value: CachedDMInterpolationInfo object (Cython wrapper) Automatically tracks hits, misses, and invalidations. @@ -52,7 +52,7 @@ def _is_enabled(self) -> bool: # Check mesh flag (default: True) return getattr(self.mesh, 'enable_dminterpolation_cache', True) - def get_structure(self, coords: np.ndarray, dofcount: int): + def get_structure(self, coords: np.ndarray, dofcount: int, policy=None): """ Get cached CachedDMInterpolationInfo or None. @@ -62,6 +62,11 @@ def get_structure(self, coords: np.ndarray, dofcount: int): Evaluation coordinates dofcount : int Total DOF count for current variables + policy : hashable, optional + Location-policy token (e.g. ``"auth"`` / ``"locate"``). A + structure built under one policy must not serve an evaluation + under the other — same coords with a different field-continuity + mix can legitimately need both. Returns ------- @@ -73,7 +78,7 @@ def get_structure(self, coords: np.ndarray, dofcount: int): # Compute cache key coords_hash = self._hash_coords(coords) - key = (coords_hash, dofcount) + key = (coords_hash, dofcount, policy) # Check cache if key in self._cache: @@ -88,7 +93,7 @@ def get_structure(self, coords: np.ndarray, dofcount: int): self._stats['misses'] += 1 return None - def store_structure(self, coords: np.ndarray, dofcount: int, cached_info): + def store_structure(self, coords: np.ndarray, dofcount: int, cached_info, policy=None): """ Store CachedDMInterpolationInfo in cache. @@ -100,12 +105,14 @@ def store_structure(self, coords: np.ndarray, dofcount: int, cached_info): Total DOF count cached_info : CachedDMInterpolationInfo Cython wrapper object containing DMInterpolation structure + policy : hashable, optional + Location-policy token — must match the ``get_structure`` lookup. """ if not self._is_enabled(): return # Don't cache if disabled coords_hash = self._hash_coords(coords) - key = (coords_hash, dofcount) + key = (coords_hash, dofcount, policy) # LRU management: remove oldest entries if we exceed max_entries if len(self._cache) >= self.max_entries: diff --git a/src/underworld3/function/petsc_tools.c b/src/underworld3/function/petsc_tools.c index 2f0240166..96936d382 100644 --- a/src/underworld3/function/petsc_tools.c +++ b/src/underworld3/function/petsc_tools.c @@ -1,4 +1,5 @@ #include "petsc_tools.h" +#include /*@C DMInterpolationSetUp - Compute spatial indices for point location during interpolation @@ -267,17 +268,16 @@ PetscErrorCode DMInterpolationEvaluate_UW(DMInterpolationInfo ctx, DM dm, Vec x, if (ctx->cells[p] < 0) { // Point couldn't be located in any cell (DMLocatePoints - // returned -1 — typically happens for query points sitting - // exactly on inter-cell boundaries or just outside the - // domain). Skipping the FE evaluation leaves - // interpolant[p*dof + ...] holding whatever was in - // PETSc-allocated memory at v's creation, which presents to - // the caller as physics-violating outliers (values ~1e+246, - // -5e+92, etc. observed at degree=2 interior nodes during - // Phase H). Zero the slot explicitly so unlocatable points - // are visible as zeros rather than garbage. + // returned -1 and no hint recovered it). Fill with NaN, not + // zero: a zero here is a plausible-looking field value that + // propagates silently (the psi* corruption behind the VEP + // blow-up, #390), while NaN is loud. The caller retrieves the + // unlocated mask (DMInterpolationGetUnlocated_UW) and fills + // these slots via the RBF fallback; NaN survives only if that + // plumbing is bypassed — which is exactly when it should be + // visible. for (PetscInt fc = 0; fc < ctx->dof; ++fc) - interpolant[p * ctx->dof + fc] = 0.0; + interpolant[p * ctx->dof + fc] = (PetscScalar)NAN; continue; } for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]); @@ -369,3 +369,22 @@ PetscErrorCode DMInterpolationEvaluate_UW(DMInterpolationInfo ctx, DM dm, Vec x, } PetscFunctionReturn(PETSC_SUCCESS); } + +/* + DMInterpolationGetUnlocated_UW - report which of the context's owned points + have no owning cell (ctx->cells[p] < 0). These are the points DMLocatePoints + dropped and no hint recovered; their interpolant slots hold NaN. The caller + uses this mask to fill those slots via the RBF fallback. + + mask must have length n == ctx->n (the number of locally owned points, which + on the evaluator's COMM_SELF context equals the number of input points, in + input order). +*/ +PetscErrorCode DMInterpolationGetUnlocated_UW(DMInterpolationInfo ctx, PetscInt n, signed char *mask) +{ + PetscFunctionBegin; + PetscCheck(n == ctx->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, + "Mask length %" PetscInt_FMT " != %" PetscInt_FMT " owned points", n, ctx->n); + for (PetscInt p = 0; p < ctx->n; ++p) mask[p] = (ctx->cells[p] < 0) ? 1 : 0; + PetscFunctionReturn(PETSC_SUCCESS); +} diff --git a/src/underworld3/function/petsc_tools.h b/src/underworld3/function/petsc_tools.h index 2618f5bae..2dd3248ac 100644 --- a/src/underworld3/function/petsc_tools.h +++ b/src/underworld3/function/petsc_tools.h @@ -6,4 +6,5 @@ #include PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain, size_t* owning_cell, PetscBool hintAuthoritative); -PetscErrorCode DMInterpolationEvaluate_UW(DMInterpolationInfo ctx, DM dm, Vec x, Vec v); \ No newline at end of file +PetscErrorCode DMInterpolationEvaluate_UW(DMInterpolationInfo ctx, DM dm, Vec x, Vec v); +PetscErrorCode DMInterpolationGetUnlocated_UW(DMInterpolationInfo ctx, PetscInt n, signed char* mask); \ No newline at end of file diff --git a/tests/test_0503_evaluate.py b/tests/test_0503_evaluate.py index 022ea8575..4ddf8b77a 100644 --- a/tests/test_0503_evaluate.py +++ b/tests/test_0503_evaluate.py @@ -145,3 +145,140 @@ def test_evaluate_on_domain_boundary_faces(simplex): ) del mesh + + +@pytest.mark.tier_a +def test_location_capability_measured(): + """The point-location capability is a measured mesh property: + rectilinear and smoothly-deformed 2D quad meshes are 'exact' (straight + edges are always planar; convexity == mesh validity), a regular 3D hex + box is 'exact', and a warped hex box drops to 'continuous'. + """ + quad = uw.meshing.StructuredQuadBox( + elementRes=(8, 4), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5) + ) + assert quad._location_capability() == "exact" + + X = quad.X.coords.copy() + defo = X.copy() + defo[:, 0] += 0.03 * np.sin(np.pi * X[:, 0]) * np.cos(2 * np.pi * X[:, 1]) + defo[:, 1] += 0.03 * np.cos(2 * np.pi * X[:, 0]) * np.sin(np.pi * X[:, 1]) + quad.deform(defo) + assert quad._location_capability() == "exact", ( + "deformed 2D quads have straight edges and convex cells - still exact" + ) + + hexbox = uw.meshing.StructuredQuadBox( + elementRes=(4, 4, 4), minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0) + ) + assert hexbox._location_capability() == "exact" + + X3 = hexbox.X.coords.copy() + d3 = X3.copy() + d3 += 0.04 * np.sin(np.pi * X3[:, [1, 2, 0]]) * np.cos(np.pi * X3[:, [2, 0, 1]]) + hexbox.deform(d3) + assert hexbox._location_capability() == "continuous", ( + "warped hex faces are non-planar - estimator authoritative for " + "continuous fields only" + ) + + del quad, hexbox + + +@pytest.mark.tier_a +def test_evaluate_deformed_quad_boundary_exact(): + """A deformed (non-affine, convex) 2D quad mesh has 'exact' capability: + boundary-face and interior evaluation of a linear field must be exact + (isoparametric Q2 reproduces linears on bilinear cells). + """ + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 4), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5) + ) + var = uw.discretisation.MeshVariable("dq_probe", mesh, 1, degree=2) + # boundary node indices captured BEFORE deforming (node identity is + # stable under a coordinate move) + idx_bnd = np.nonzero( + (np.abs(var.coords[:, 1]) > 0.5 - 1e-12) + | (np.abs(var.coords[:, 0]) > 1.0 - 1e-12) + )[0] + + X = mesh.X.coords.copy() + defo = X.copy() + defo[:, 0] += 0.03 * np.sin(np.pi * X[:, 0]) * np.cos(2 * np.pi * X[:, 1]) + defo[:, 1] += 0.03 * np.cos(2 * np.pi * X[:, 0]) * np.sin(np.pi * X[:, 1]) + mesh.deform(defo) + + var.data[:, 0] = 1.0 + 2.0 * var.coords[:, 0] + 3.0 * var.coords[:, 1] + + # interior points plus points ON the (deformed) boundary surface + bnd = var.coords[idx_bnd[:: max(1, len(idx_bnd) // 16)]] + interior = np.column_stack([ + np.linspace(-0.8, 0.8, 9), np.linspace(-0.35, 0.35, 9) + ]) + pts = np.vstack([interior, bnd]) + vals = fn.evaluate(var.sym[0], pts).flatten() + exact = 1.0 + 2.0 * pts[:, 0] + 3.0 * pts[:, 1] + assert np.allclose(vals, exact, atol=1e-8), ( + f"max|err|={np.abs(vals - exact).max():.3e}" + ) + + del mesh + + +@pytest.mark.tier_a +def test_evaluate_warped_hex_rbf_fallback_no_silent_values(): + """On a warped hex box ('continuous' capability) a DISCONTINUOUS variable + must not trust the cell-wall hint: dropped boundary-face points take the + RBF fallback. The result must be finite (no NaN sentinel escaping, no + silent zeros) and bounded by the field's data range; a CONTINUOUS + variable on the same mesh stays on the authoritative path and evaluates + a linear field exactly (trilinear isoparametric reproduces linears). + """ + mesh = uw.meshing.StructuredQuadBox( + elementRes=(4, 4, 4), minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0) + ) + X3 = mesh.X.coords.copy() + d3 = X3.copy() + d3 += 0.04 * np.sin(np.pi * X3[:, [1, 2, 0]]) * np.cos(np.pi * X3[:, [2, 0, 1]]) + mesh.deform(d3) + assert mesh._location_capability() == "continuous" + + cvar = uw.discretisation.MeshVariable("wh_cont", mesh, 1, degree=1) + cvar.data[:, 0] = 1.0 + 2.0 * cvar.coords[:, 0] - 1.5 * cvar.coords[:, 2] + dvar = uw.discretisation.MeshVariable( + "wh_disc", mesh, 1, degree=1, continuous=False + ) + dvar.data[:, 0] = 1.0 + 2.0 * dvar.coords[:, 0] - 1.5 * dvar.coords[:, 2] + + # query points on the deformed top surface (actual boundary node coords) + zmax_nodes = cvar.coords[ + cvar.coords[:, 2] > np.percentile(cvar.coords[:, 2], 92) + ][:15] + interior = np.column_stack([ + np.linspace(0.15, 0.85, 8), + np.linspace(0.2, 0.8, 8), + np.linspace(0.25, 0.75, 8), + ]) + pts = np.vstack([interior, zmax_nodes]) + + # continuous: authoritative path, linear field reproduced exactly + cvals = fn.evaluate(cvar.sym[0], pts).flatten() + cexact = 1.0 + 2.0 * pts[:, 0] - 1.5 * pts[:, 2] + assert np.allclose(cvals, cexact, atol=1e-8), ( + f"continuous-field max|err|={np.abs(cvals - cexact).max():.3e}" + ) + + # discontinuous: locate + RBF fallback for dropped points — finite, + # in-range, and close to the (smooth) nodal data it interpolates + dvals = fn.evaluate(dvar.sym[0], pts).flatten() + dexact = 1.0 + 2.0 * pts[:, 0] - 1.5 * pts[:, 2] + assert np.all(np.isfinite(dvals)), "NaN escaped the RBF fallback" + lo, hi = dvar.data.min(), dvar.data.max() + assert np.all(dvals >= lo - 1e-12) and np.all(dvals <= hi + 1e-12), ( + "fallback values outside the field's data range" + ) + assert np.abs(dvals - dexact).max() < 0.25, ( + f"fallback error unexpectedly large: {np.abs(dvals - dexact).max():.3e}" + ) + + del mesh