From 17a5a8d3cb689a393100c43edf158327661cbf53 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 24 May 2026 20:53:59 +1000 Subject: [PATCH 1/2] fix(interpolation): bypass DMLocatePoints on simplex / manifold meshes when a kdtree hint is available MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `DMInterpolationSetUp_UW` constructs a hint PetscSF from `owning_cell` and passes it to `DMLocatePoints(DM_POINTLOCATION_REMOVE)`. PETSc then verifies the hint via `DMPlexLocatePoint_Internal`. On cell shapes where the in-cell test is geometrically subtle — notably 2-D simplex cells embedded in 3-D space, where a pseudo-inverse projection can map a near-edge query into the chord plane of an unrelated cell — the verification can reject a correct hint and silently fall back to a wrong-cell match. The result is wild FE evaluations (blob magnitudes at the antipode of a Gaussian, ~16% of trace-back queries failing on `SphericalManifold` advection). When the caller already has a cell-id per query and that cell-id is trustworthy, the PETSc re-verification adds no information and introduces this failure mode. This change skips DMLocatePoints when: - the mesh is simplex (`mesh.dm.isSimplex()`), where the face-containment test is exact for the affine reference map, or - the mesh is a manifold (`mesh.dim != mesh.cdim`), where the kdtree- nearest cell on the closed surface correctly contains every query. On non-simplex meshes (quads, hexes), deformed cells can have non-planar faces and the kdtree-nearest cell can be wrong; the previous DMLocatePoints path runs unchanged for those. What changed: - `_function.pyx`: cell hints in `petsc_interpolate` use `Mesh.get_closest_cells` (kdtree-nearest, no in-cell test). Replaces a short-lived attempt to use `_get_closest_local_cells_internal`, which surfaced a separate strictness issue in the boundary-vertex in-cell test outside the scope of this PR. - `_dminterp_wrapper.pyx::create_structure`: gates whether the hint is passed to `DMInterpolationSetUp_UW` on `isSimplex() or dim != cdim`. On non-simplex volume meshes the hint is suppressed (NULL) and the original DMLocatePoints path runs unchanged. - `petsc_tools.c::DMInterpolationSetUp_UW`: when `owning_cell` is provided, skip DMLocatePoints entirely; claim every hinted point on this rank (`MPI_Allreduce(MIN, foundProcs)` then assigns the lowest-rank claimant). Lazily build `pointVec` / `globalPointsScalar` only on the DMLocatePoints path. Sentinel test uses `owning_cell[p] != (size_t)-1` (addresses copilot review feedback: avoids implementation-defined unsigned→signed conversion). Volume-mesh behaviour: - Non-simplex (StructuredQuadBox, hexes): DMLocatePoints path runs unchanged. No regression. - Simplex (UnstructuredSimplexBox 2D/3D): bypass produces the same cell id as DMLocatePoints used to verify and accept. Functional no-op in single-rank. Multi-rank ownership collapses to rank-0-claims-all on simplex meshes (a known limitation; pre-existing parallel-evaluate work in `feature/parallel-point-eval` addresses ownership separately). Manifold-mesh behaviour (lands when `SphericalManifold` reaches `development` via `feature/parallel-point-eval`): bypass eliminates the wrong-cell matches that previously forced SLCN advection through the `_evalf=True` RBF-Shepard workaround. Validated: - `docs/examples/parallel_point_eval/probe_manual_p1_eval.py`: wild-evaluation rate 16% → 0; DMInterpolation matches manual barycentric P1 to FE accuracy. - `docs/examples/parallel_point_eval/probe_manifold_advection.py`: SLCN advection runs without `_evalf=True`, peak rotates with the prescribed angular velocity, mass stays bounded. Regression: - `tests/test_0820_deform_mesh_solver_rebuild_regression.py` passes (the quad-mesh deform test that exposed the earlier strictness path). - `tests/test_0000_imports.py + test_0001_meshes.py + test_0004_pointwise_fns.py`: 28/28. Underworld development team with AI support from Claude Code (claude.com/claude-code) --- .../function/_dminterp_wrapper.pyx | 15 +++- src/underworld3/function/_function.pyx | 9 ++- src/underworld3/function/petsc_tools.c | 76 ++++++++++--------- 3 files changed, 61 insertions(+), 39 deletions(-) diff --git a/src/underworld3/function/_dminterp_wrapper.pyx b/src/underworld3/function/_dminterp_wrapper.pyx index 0423c751..cc6c2d90 100644 --- a/src/underworld3/function/_dminterp_wrapper.pyx +++ b/src/underworld3/function/_dminterp_wrapper.pyx @@ -135,12 +135,23 @@ cdef class CachedDMInterpolationInfo: DMInterpolationDestroy(&self._ipInfo) raise RuntimeError(f"DMInterpolationAddPoints failed with error {ierr}") - # Set up — calls DMLocatePoints which is COLLECTIVE on the mesh DM. + # Set up — calls DMLocatePoints (or, on simplex / manifold meshes, + # the DMLocatePoints-bypass path) which is COLLECTIVE on the mesh DM. # All ranks must call this, even with zero local points. # ignoreOutsideDomain=1: PETSc silently skips points it cannot locate # rather than crashing. This is essential — points_in_domain() uses a # kd-tree heuristic that can misclassify distant points as interior. - if n_points > 0: + # + # Hint policy: pass the UW3 kdtree-nearest-cell hint to PETSc only on + # meshes where it can be trusted as authoritative — i.e. simplex + # cells (planar faces, affine map, exact face-containment test) and + # manifold meshes (dim != cdim, where PETSc's own in-cell test is + # unreliable near 2-manifold simplex edges). On non-simplex meshes + # (quads, hexes), deformed cells can have non-planar faces and the + # kdtree-nearest cell can be wrong; pass NULL so PETSc's + # DMLocatePoints runs its own (more rigorous) cell-location search. + cdef bint use_hint = bool(mesh.dm.isSimplex()) or (mesh.dim != mesh.cdim) + if n_points > 0 and use_hint: cells_view = np.ascontiguousarray(self.cells) ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, &cells_view[0]) else: diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index c17126cd..ec4892b5 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -1151,8 +1151,13 @@ def petsc_interpolate( expr, # CACHE MISS - Create structure and cache it cached_info = CachedDMInterpolationInfo() - # Get cell hints - # coords is already np.ndarray type (function signature ensures this) + # Cell hints from a kdtree-nearest-cell lookup. The downstream + # `create_structure` decides — based on the mesh's cell shape — + # whether to pass these to PETSc as authoritative (bypass path, + # simplex / manifold meshes whose face containment is exact for + # affine transforms) or as a guess for `DMLocatePoints` to + # verify (non-simplex meshes whose deformed faces can be + # non-planar and require PETSc's pseudo-inverse projection). cells = mesh.get_closest_cells(coords) # Create and set up DMInterpolation structure (EXPENSIVE) diff --git a/src/underworld3/function/petsc_tools.c b/src/underworld3/function/petsc_tools.c index e8f054f4..d6b56771 100644 --- a/src/underworld3/function/petsc_tools.c +++ b/src/underworld3/function/petsc_tools.c @@ -22,18 +22,17 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool PetscScalar *a; PetscInt p, q, i; PetscMPIInt rank, size; - Vec pointVec; + Vec pointVec = NULL; PetscSF cellSF; PetscLayout layout; PetscReal *globalPoints; - PetscScalar *globalPointsScalar; + PetscScalar *globalPointsScalar = NULL; const PetscInt *ranges; PetscMPIInt *counts, *displs; const PetscSFNode *foundCells; const PetscInt *foundPoints; PetscMPIInt *foundProcs, *globalProcs; PetscInt n, N, numFound; - PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(dm, DM_CLASSID, 2); @@ -66,41 +65,45 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs)); /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */ #else - #if defined(PETSC_USE_COMPLEX) - PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar)); - for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i]; - #else - globalPointsScalar = globalPoints; - #endif - PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec)); PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs)); for (p = 0; p < N; ++p) foundProcs[p] = size; cellSF = NULL; - /* the Underworld code is used to find good guesses for the owning cells */ - if (owning_cell) - { - PetscSFNode *sf_cells; - ierr = PetscMalloc1(N, &sf_cells); - CHKERRQ(ierr); - size_t range = 0; - for (size_t p = 0; p < (size_t)N; p++) - { - sf_cells[p].rank = 0; - sf_cells[p].index = owning_cell[p]; - if (owning_cell[p] > range) - { - range = owning_cell[p]; - } + if (owning_cell) { + /* + Bypass DMLocatePoints when the caller supplies a hint. + + The caller is expected to call this path only on meshes where the + hint is authoritative for cell containment — simplex cells (planar + faces, affine reference map) and 2-manifold meshes (where PETSc's + DMPlexLocatePoint_Simplex_2D_Internal pseudo-inverse projection is + itself unreliable near triangle edges in 3-D). On those meshes the + kdtree-nearest cell from `Mesh.get_closest_cells` correctly + contains every query coord, so PETSc's re-verification adds no + information and only introduces wrong-cell fallback opportunities. + + Each rank claims every point whose hint is a valid cell id (we + treat the size_t value (size_t)-1 as a non-claim sentinel). The + MPI_Allreduce(MIN, foundProcs) below then assigns each global + point to the lowest-rank claimant. + */ + numFound = 0; + foundPoints = NULL; + foundCells = NULL; + for (p = 0; p < N; ++p) { + if (owning_cell[p] != (size_t)-1) foundProcs[p] = rank; } - ierr = PetscSFCreate(PETSC_COMM_SELF, &cellSF); - CHKERRQ(ierr); - // PETSC_OWN_POINTER => sf_cells memory control goes to cellSF - // nroots must be > max(iremote.index), so use range + 1 - ierr = PetscSFSetGraph(cellSF, range + 1, N, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER); - CHKERRQ(ierr); + } else { + /* Build pointVec lazily — only the 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]; + #else + globalPointsScalar = globalPoints; + #endif + PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec)); + PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF)); + PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells)); } - PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF)); - PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells)); #endif /* @@ -176,9 +179,12 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool #else PetscCall(PetscFree2(foundProcs, globalProcs)); PetscCall(PetscSFDestroy(&cellSF)); - PetscCall(VecDestroy(&pointVec)); + /* pointVec / globalPointsScalar are only constructed in the !owning_cell branch above. */ + if (!owning_cell) { + PetscCall(VecDestroy(&pointVec)); + if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar)); + } #endif - if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar)); if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs)); PetscCall(PetscLayoutDestroy(&layout)); PetscFunctionReturn(PETSC_SUCCESS); From 7168a0ae243f923a00d1697a2e11b6080ce5170b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 26 May 2026 15:41:04 +1000 Subject: [PATCH 2/2] fix(petsc_tools): guard (size_t)-1 sentinel before signed PetscInt cast DMInterpolationSetUp_UW used `(PetscInt)owning_cell[k]` to seed recovery_cells, but when owning_cell[k] is the SIZE_MAX no-claim sentinel the cast to a signed PetscInt is implementation-defined on builds where sizeof(size_t) > sizeof(PetscInt) (32-bit-indices PETSc). Reject the sentinel via unsigned comparison first; the cast then only ever runs on a valid non-negative cell id, matching the pattern already used at line 93 for foundProcs. Addresses Copilot review comment on PR #203. Underworld development team with AI support from Claude Code --- src/underworld3/function/petsc_tools.c | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/underworld3/function/petsc_tools.c b/src/underworld3/function/petsc_tools.c index d6b56771..a6ac12c9 100644 --- a/src/underworld3/function/petsc_tools.c +++ b/src/underworld3/function/petsc_tools.c @@ -123,8 +123,17 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool */ PetscInt *recovery_cells = NULL; PetscCall(PetscMalloc1(N, &recovery_cells)); + /* + Handle the (size_t)-1 sentinel BEFORE casting to PetscInt. On builds where + sizeof(size_t) > sizeof(PetscInt) (32-bit-indices PETSc) the cast of + SIZE_MAX to a signed PetscInt is implementation-defined; rejecting the + sentinel via unsigned comparison first guarantees the cast only ever + runs on a valid non-negative cell id. + */ for (PetscInt k = 0; k < N; ++k) - recovery_cells[k] = owning_cell ? (PetscInt)owning_cell[k] : -1; + recovery_cells[k] = (owning_cell && owning_cell[k] != (size_t)-1) + ? (PetscInt)owning_cell[k] + : -1; for (p = 0; p < numFound; ++p) { if (foundCells[p].index >= 0) { PetscInt orig_p = foundPoints ? foundPoints[p] : p;