fix(evaluate): recover boundary-face query points dropped by DMLocatePoints on quad/hex meshes (#390)

src/underworld3/function/_dminterp_wrapper.pyx

@@ -34,7 +34,7 @@ cdef extern from "petsc.h" nogil:
 
 # Import custom UW routines
 cdef extern from "petsc_tools.h" nogil:
-    PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ipInfo, PetscDM dm, int petscbool, int petscbool, size_t* owning_cell)
+    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)
 
 cdef class CachedDMInterpolationInfo:
@@ -141,31 +141,37 @@ cdef class CachedDMInterpolationInfo:
         # ignoreOutsideDomain=1: PETSc silently skips points it cannot locate
         # rather than crashing.
         #
-        # Hint / bypass policy lives in ONE place: mesh._eval_use_robust_location().
-        # When it is True (parallel simplex / manifold), pass the bulletproof
-        # barycentric hint as authoritative -> petsc_tools.c bypasses
-        # DMLocatePoints (ported from 17a5a8d) and evaluates in the hinted cell,
-        # avoiding PETSc's slow, on-face-rejecting, error-raising location.
-        # Otherwise pass NULL -> PETSc DMLocatePoints runs (serial bit-identical
-        # baseline; non-simplex volume meshes whose deformed faces need PETSc's
-        # rigorous search).
-        # The DMLocatePoints-BYPASS (use the supplied cell hint instead of
-        # PETSc's slow, on-face-rejecting search) is safe whenever the hint is
-        # *authoritative*: simplex cells (planar faces, affine reference map) or
-        # manifold meshes (dim != cdim). That property is independent of rank
-        # count, so the bypass applies in SERIAL too — restoring 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,
-        # which are a valid hint (the pre-merge serial behaviour, bit-for-bit).
-        # Non-simplex volume meshes (deformed quad/hex faces) keep NULL ->
-        # PETSc DMLocatePoints, where the kdtree-nearest hint can be wrong.
-        cdef bint use_hint = (bool(mesh.dm.isSimplex()) or (mesh.dim != mesh.cdim))
-        if n_points > 0 and use_hint:
+        # Hint policy — the hint is ALWAYS passed; what varies is its authority
+        # (the trailing hintAuthoritative flag):
+        #
+        # 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.
+        #
+        # 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:
             cells_view = np.ascontiguousarray(self.cells)
-            ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, <size_t*> &cells_view[0])
+            ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1,
+                                           <size_t*> &cells_view[0],
+                                           1 if hint_authoritative else 0)
         else:
-            ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, NULL)
+            ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, NULL, 0)
         if ierr != 0:
             DMInterpolationDestroy(&self._ipInfo)
             raise RuntimeError(f"DMInterpolationSetUp_UW failed with error {ierr}")

src/underworld3/function/_function.pyx

@@ -32,7 +32,7 @@ cdef extern from "petsc.h" nogil:
         pass
 
 cdef extern from "petsc_tools.h" nogil:
-    PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ipInfo, PetscDM dm, int petscbool, int petscbool, size_t* owning_cell)
+    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)
 
 cdef extern from "petsc.h" nogil:

src/underworld3/function/petsc_tools.c

@@ -16,7 +16,7 @@
 .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
 @*/
 
-PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain, size_t *owning_cell)
+PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain, size_t *owning_cell, PetscBool hintAuthoritative)
 {
   MPI_Comm           comm = ctx->comm;
   PetscScalar       *a;
@@ -68,10 +68,10 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool
   PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
   for (p = 0; p < N; ++p) foundProcs[p] = size;
   cellSF = NULL;
-  if (owning_cell) {
+  if (owning_cell && hintAuthoritative) {
     /*
-      Bypass DMLocatePoints when the caller supplies a hint (ported from
-      feature/dminterp-bypass-element-check, 17a5a8d).
+      Bypass DMLocatePoints when the caller supplies an AUTHORITATIVE hint
+      (ported from feature/dminterp-bypass-element-check, 17a5a8d).
 
       The caller is expected to call this path only on meshes where the
       hint is authoritative for cell containment — simplex cells (planar
@@ -97,7 +97,18 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool
       if (owning_cell[p] != (size_t)-1) foundProcs[p] = rank;
     }
   } else {
-    /* Build pointVec lazily — only the DMLocatePoints path needs it. */
+    /*
+      DMLocatePoints is authoritative. A non-authoritative `owning_cell`
+      hint (if supplied) is NOT used to bypass the search — it only
+      prefills `recovery_cells` below, rescuing points that PETSc's
+      _REMOVE flag silently drops (e.g. query points sitting exactly on
+      the domain's closed upper faces, where the grid-hash / in-cell
+      convention is half-open). Without that rescue such points fall
+      through to the explicit zero-fill in DMInterpolationEvaluate_UW —
+      the silent ψ*-corruption behind the VEP stability blow-up (#390).
+
+      Build pointVec lazily — only this DMLocatePoints path needs it.
+    */
     #if defined(PETSC_USE_COMPLEX)
     PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
     for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
@@ -183,9 +194,9 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool
 #else
   PetscCall(PetscFree2(foundProcs, globalProcs));
   PetscCall(PetscSFDestroy(&cellSF));
-  /* pointVec / globalPointsScalar are only constructed on the !owning_cell
-     (DMLocatePoints) branch above. */
-  if (!owning_cell) {
+  /* pointVec / globalPointsScalar are only constructed on the
+     DMLocatePoints branch above (pointVec stays NULL on the bypass). */
+  if (pointVec) {
     PetscCall(VecDestroy(&pointVec));
     if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
   }

src/underworld3/function/petsc_tools.h

@@ -5,5 +5,5 @@
 #include <petsc/private/petscimpl.h>
 #include <petsc/private/petscfeimpl.h>
 
-PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain, size_t* owning_cell);
+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);

tests/test_0503_evaluate.py

@@ -100,3 +100,48 @@ def test_single_vector_variable():
     assert np.allclose(np.array(((1.1, 1.2),)), result, rtol=1e-05, atol=1e-08)
 
     del mesh
+
+
+@pytest.mark.tier_a
+@pytest.mark.parametrize("simplex", [False, True], ids=["quad", "simplex"])
+def test_evaluate_on_domain_boundary_faces(simplex):
+    """Interpolation at points lying exactly ON the domain boundary must be
+    exact, on both the DMLocatePoints path (quad) and the hint-bypass path
+    (simplex).
+
+    PETSc's point location drops queries sitting exactly on the domain's
+    closed upper faces (half-open cell convention), and a dropped point is
+    zero-filled by the evaluator unless the kd-tree cell hint recovers it.
+    Losing that recovery on quad meshes silently corrupted semi-Lagrangian
+    stress-history trace-backs along the top boundary — the VEP stability
+    blow-up of issue #390.
+    """
+    if simplex:
+        mesh = uw.meshing.UnstructuredSimplexBox(
+            minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5), cellSize=0.25
+        )
+    else:
+        mesh = uw.meshing.StructuredQuadBox(
+            elementRes=(8, 4), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5)
+        )
+    var = uw.discretisation.MeshVariable("bdry_probe", mesh, 1, degree=2)
+    var.data[:, 0] = 1.0 + 2.0 * var.coords[:, 0] + 3.0 * var.coords[:, 1]
+
+    xs = np.linspace(-0.9, 0.9, 7)
+    ys = np.linspace(-0.4, 0.4, 5)
+    pts = np.array(
+        [[xi, 0.5] for xi in xs]        # top face (the face PETSc drops)
+        + [[xi, -0.5] for xi in xs]     # bottom face
+        + [[-1.0, yi] for yi in ys]     # left face
+        + [[1.0, yi] for yi in ys]      # right face
+        + [[-1.0, -0.5], [1.0, -0.5], [-1.0, 0.5], [1.0, 0.5]]  # corners
+    )
+    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"boundary-face interpolation error: "
+        f"max|err|={np.abs(vals - exact).max():.3e} "
+        f"worst at {pts[np.argmax(np.abs(vals - exact))]}"
+    )
+
+    del mesh