diff --git a/src/underworld3/function/_dminterp_wrapper.pyx b/src/underworld3/function/_dminterp_wrapper.pyx index c6bfc2c8f..b208fae5d 100644 --- a/src/underworld3/function/_dminterp_wrapper.pyx +++ b/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, &cells_view[0]) + ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, + &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}") diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index c1d1bce3e..607cb388c 100644 --- a/src/underworld3/function/_function.pyx +++ b/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: diff --git a/src/underworld3/function/petsc_tools.c b/src/underworld3/function/petsc_tools.c index 473701d53..2f0240166 100644 --- a/src/underworld3/function/petsc_tools.c +++ b/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)); } diff --git a/src/underworld3/function/petsc_tools.h b/src/underworld3/function/petsc_tools.h index cab3cd695..2618f5bae 100644 --- a/src/underworld3/function/petsc_tools.h +++ b/src/underworld3/function/petsc_tools.h @@ -5,5 +5,5 @@ #include #include -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); \ No newline at end of file diff --git a/tests/test_0503_evaluate.py b/tests/test_0503_evaluate.py index fd6609add..022ea8575 100644 --- a/tests/test_0503_evaluate.py +++ b/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