fix(interpolation): bypass DMLocatePoints when caller supplies a verified hint#203
Open
lmoresi wants to merge 1 commit into
Open
fix(interpolation): bypass DMLocatePoints when caller supplies a verified hint#203lmoresi wants to merge 1 commit into
lmoresi wants to merge 1 commit into
Conversation
…fied hint `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 — most notably 2-D simplex cells embedded in 3-D space, where the 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 verified cell-id per query coord, the PETSc re-verification adds no information and introduces this failure mode. This change skips DMLocatePoints when `owning_cell` is provided and treats the hint as authoritative. The hint source in `petsc_interpolate` switches from `Mesh.get_closest_cells` (kdtree-only) to `Mesh._get_closest_local_cells_internal`, which: - runs the in-cell test on the kdtree-nearest pick, - falls back to walking nearby cells when that pick fails, - returns `-1` for points not owned by this rank. `DMInterpolationSetUp_UW` reads the `-1` sentinel and only claims points owned locally; the existing `MPI_Allreduce(MIN, foundProcs)` then assigns each global point to its single owning rank, preserving the previous ownership convention. On volume meshes this is a behavioural no-op: the hint always passed PETSc's in-cell test under the old path and is the same cell now. On 2-manifold meshes (`SphericalManifold`, landing with PR #N [feature/parallel-point-eval]) the bypass eliminates the wrong-cell matches that previously forced SLCN advection through the `_evalf=True` RBF-Shepard workaround. Validated via `docs/examples/parallel_point_eval/probe_manual_p1_eval.py` on a `SphericalManifold`: wild-evaluation rate 16% → 0; `uw.function.evaluate` now matches manual barycentric P1 to FE accuracy. Underworld development team with AI support from Claude Code (claude.com/claude-code)
Contributor
There was a problem hiding this comment.
Pull request overview
This PR adjusts Underworld3’s PETSc-based point interpolation path to trust caller-supplied, already-verified cell hints (instead of re-validating via DMLocatePoints()), aiming to eliminate wrong-cell fallbacks that can produce extreme (“wild”) FE evaluations on geometrically subtle cell configurations (notably 2-manifold simplices embedded in 3D).
Changes:
- Switches Python-side hint generation in
petsc_interpolatefrom a KDTree-only closest-cell guess to a locally verified in-cell search (Mesh._get_closest_local_cells_internal, returning-1when not locally owned). - Updates
DMInterpolationSetUp_UWto bypassDMLocatePoints()entirely whenowning_cellis provided, claiming points based on the hint instead.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 2 comments.
| File | Description |
|---|---|
src/underworld3/function/_function.pyx |
Generates interpolation cell hints using an in-cell-tested local search, enabling safe trust of hints downstream. |
src/underworld3/function/petsc_tools.c |
Bypasses DMLocatePoints() when hints are provided, using the hint to assign point ownership and cell IDs. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
68
to
+79
| #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) { | ||
| /* |
Comment on lines
+104
to
+106
| /* owning_cell carries int64 reinterpreted as size_t; -1 → SIZE_MAX. | ||
| Casting back to (signed) PetscInt recovers the -1 sentinel. */ | ||
| if ((PetscInt)owning_cell[p] >= 0) foundProcs[p] = rank; |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
DMInterpolationSetUp_UW(UW3-side fork ofDMInterpolationSetUp) currently builds a hint PetscSF fromowning_celland hands it toDMLocatePoints(DM_POINTLOCATION_REMOVE), which re-verifies each hint viaDMPlexLocatePoint_Internal.DMLocatePointswhenowning_cellis provided and trusts the hint.What's in the patch
src/underworld3/function/_function.pyx— switch the hint source inpetsc_interpolatefromMesh.get_closest_cells(kdtree-only, no ownership filtering) toMesh._get_closest_local_cells_internal(runs the in-cell test, walks nearby cells on miss, returns-1for non-owned points).src/underworld3/function/petsc_tools.c— inDMInterpolationSetUp_UW, replace theDMLocatePoints + cellSF-as-hintpath with a direct loop that claimsfoundProcs[p] = rankiffowning_cell[p] >= 0. The existingMPI_Allreduce(MIN, foundProcs)then assigns each global point to its single owning rank, preserving the prior ownership convention.Net: +40 / -27 lines across the two files.
Why this PR targets feature/parallel-point-eval
The bypass requires
Mesh._get_closest_local_cells_internalto be reliable on this branch's mesh inventory — including the in-cell test rejecting only points genuinely not on the rank, not boundary vertices.On plain
development,_test_if_points_in_cells_internalis strict enough to reject ~43% of vertex points on a basic structured quad mesh (probed locally), which means the hint comes back as-1for legitimately-owned boundary vertices and the bypass drops them. This was not previously exposed becauseMesh.get_closest_cells(the prior hint source) is purely kdtree-nearest and never returns-1.feature/parallel-point-evalhas the manifold-aware perpendicular construction (cell_normal × edge) in_mark_faces_inside_and_out, and the in-cell test there handles the boundary-vertex case correctly. So this PR targets PPE and will land alongside it.Behaviour expectations
-1cell — the bypass produces the samectx->cellsas the previous DMLocatePoints-verified path. Behavioural no-op.SphericalManifold, on this branch): bypass eliminates wrong-cell matches that previously forced SLCN advection through the_evalf=TrueRBF-Shepard workaround.Test plan
./uw buildon this branchuw.function.evaluateon 2D P1, 2D P2, 3D P1, and 2-component vector — all match analytic valuespytest tests/test_0000_imports.py tests/test_0001_meshes.py tests/test_0004_pointwise_fns.py— 28/28 passdocs/examples/parallel_point_eval/probe_manual_p1_eval.pyon aSphericalManifold: wild-evaluation rate 16% → 0; DMInterpolation matches manual barycentric P1 to FE accuracydocs/examples/parallel_point_eval/probe_manifold_advection.py: SLCN advection runs without_evalf=True, peak rotates with the prescribed angular velocity, mass stays boundedKnown limitations / out of scope
owning_cell[p]access forp >= n_localinredundantPoints=FALSEmode (where N fromMPI_Allgathervcan exceed each rank's hint-array length) is a separate latent bug, not touched by this PR. Flagged for follow-up.DMInterpolationSetCellHintsto PETSc'sDMInterpolationAPI (the stock wrapper unconditionally drops caller hints — fixing that is useful for non-UW3 PETSc users but unrelated to this UW3 bypass).Underworld development team with AI support from Claude Code