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
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
| 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; |
…s when a kdtree hint is available
`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)
lmoresi
force-pushed
the
feature/dminterp-bypass-element-check
branch
from
97b5b8b to
17a5a8d
Compare
Member
Author
|
Force-pushed an updated design (commit The earlier push had a regression on quad meshes — Both Copilot review comments addressed in the new version (pointVec lazy-built, unsigned-sentinel comparison). |
6 tasks
Member
Author
|
Update — precursor needed Discussed offline: The right authoritative hint source is Once #207 merges, I'll rebase this PR onto it and switch the hint source in |
lmoresi
added a commit
that referenced
this pull request
May 25, 2026
…mantics `Mesh._test_if_points_in_cells_internal` used a strict `> 0` test on the squared-distance difference between mirrored inner/outer control points placed ±1e-3 along each face normal. A query exactly on a cell face has zero distance difference and was rejected. That meant mesh vertices — which sit on the faces of every cell containing them in their closure — failed the in-cell test for every candidate cell, and `_get_closest_local_cells_internal` returned -1 for them. For evaluation use cases this is the wrong semantics: the FE basis at a shared vertex / face is consistent across the adjacent cells (a DM consistency requirement of FE assembly), so "any cell whose closure contains the point" is correct. Picking one specific adjacent cell and returning its id lets downstream code evaluate correctly. Add an opt-in `on_boundary: bool = False` kwarg that flips the comparison from strict `> 0` to loose `>= -1e-12`. Default is False (preserves historical strict-inside behaviour); callers that want loose semantics must request them explicitly. The kwarg is forwarded through the call chain: - `_test_if_points_in_cells_internal(on_boundary=False)` — core test - `_get_closest_local_cells_internal(on_boundary=False)` — forwards - `test_if_points_in_cells(on_boundary=False)` — public wrapper, forwards - `get_closest_local_cells(on_boundary=False)` — public wrapper, forwards Backward-compatibility rationale: callers like `points_in_domain` (used in `petsc_interpolate` to split queries into interior / exterior for FE vs. RBF routing) depend on the strict semantics — a vertex on the domain boundary is classified as "outside" and routed to RBF extrapolation. Changing that default would change `uw.function.evaluate` output values for any code that evaluates at mesh vertices (broke `test_advDiff_boxmesh[mesh0]` in a first attempt). Opt-in avoids the regression while making the loose semantics available where they're genuinely needed. This unblocks the bypass design in PR #203 (`feature/dminterp-bypass-element-check`), which needs an authoritative per-rank cell-id source. With `on_boundary=True`, `_get_closest_local_cells_internal` returns a cell whose closure contains the query (or -1 if no local cell does) — exactly what the bypass requires. Also addresses Copilot review feedback: - docstring corrected to say `_test_if_points_in_cells_internal` expects model-coords input (the public wrapper handles units). - removed an unused `inside` initialisation left over from the strict/loose refactor. - public `test_if_points_in_cells` now coerces `cells` to a 1-D numpy array so list/tuple input works as the docstring promises. Tests `tests/test_0820_in_cell_test_loose_semantics.py` (7 tests) lock the contract: - default strict-inside rejects boundary-vertex queries - on_boundary=True accepts every vertex query on 2D simplex, 3D simplex, and 2D structured-quad meshes - strict and loose modes give a distinguishable result on the same vertex queries - cells returned with on_boundary=True genuinely contain the query in their closure - public `get_closest_local_cells` forwards the kwarg Underworld development team with AI support from Claude Code (claude.com/claude-code)
Member
Author
|
Update on the precursor: #207 now uses an opt-in design — When #207 merges, I'll update this PR's |
lmoresi
added a commit
that referenced
this pull request
May 25, 2026
…ries by default
`Mesh._test_if_points_in_cells_internal` used a strict `> 0` test on the
squared-distance difference between mirrored inner/outer control points
placed ±1e-3 along each face normal. A query exactly on a cell face has
zero distance difference and was rejected. That meant mesh vertices —
which sit on the faces of every cell containing them in their closure —
failed the in-cell test for every candidate cell, and
`_get_closest_local_cells_internal` returned -1 for them.
For evaluation use cases this is the wrong semantics: a closed domain Ω
includes ∂Ω, and the FE basis at a shared vertex / face is consistent
across the adjacent cells (a DM consistency requirement of FE assembly).
So "any cell whose closure contains the point" is the right notion.
Routing those queries through RBF Shepard extrapolation (the legacy
behaviour of `uw.function.evaluate` for points classified outside the
domain) is less accurate than evaluating via FE on any adjacent cell.
Add `on_boundary: bool = True` kwarg to the in-cell test and forward
through the call chain:
- `_test_if_points_in_cells_internal(on_boundary=True)` — core test
- `_get_closest_local_cells_internal(on_boundary=True)` — forwards
- `test_if_points_in_cells(on_boundary=True)` — public wrapper
- `get_closest_local_cells(on_boundary=True)` — public wrapper
With on_boundary=True (the default) the comparison is `>= -1e-12` (well
below the 1e-3 control-point offset, well above 64-bit float roundoff);
with on_boundary=False the historical `> 0` strict-inside semantics is
preserved. Callers that need uniqueness (a point claimed by exactly one
cell, never a shared-face point claimed by both adjacent cells) can opt
back into strict.
Observable consequences:
- `_get_closest_local_cells_internal` is now an authoritative cell-id
source: returns a cell whose closure contains the query, or -1 if no
local cell does.
- `points_in_domain` reports boundary-vertex queries as in the domain
(matches mathematical convention; a point on ∂Ω is in Ω).
- `uw.function.evaluate` now routes boundary-vertex queries through FE
instead of RBF Shepard — the more accurate path. This shifts the
output of `evaluate` at mesh vertices on the domain boundary.
Test impact:
- `tests/test_0820_in_cell_test_loose_semantics.py` (new, 7 tests):
locks the on_boundary=True default — every vertex of 2D simplex, 3D
simplex, and 2D quad meshes resolves to a containing cell; strict
mode still rejects on-face queries; loose mode returns cells whose
closure genuinely contains the query.
- `tests/test_0820_deform_mesh_solver_rebuild_regression.py` continues
to pass (it was the motivating regression for the PR #203 bypass
work that depends on this).
- `tests/test_1100_AdvDiffCartesian.py::test_advDiff_boxmesh[mesh0]`
marked xfail (strict=True). The file header already calls this "not
a great test"; its atol=0.05 tolerance previously aligned only
because boundary-vertex queries went through RBF Shepard's
smoothing, and the test was tuned to that result. Needs reworking
(smoother IC, larger transport distance) to pass under the more
accurate FE-evaluate path. The two simplex variants of the same
test (mesh1, mesh2) continue to pass.
Also addresses Copilot review feedback:
- docstring on `_test_if_points_in_cells_internal` corrected to say
model-coords input (the public wrapper handles units).
- removed an unused `inside` initialisation from the refactor.
- public `test_if_points_in_cells` now coerces `cells` to a 1-D numpy
array so list/tuple input works as the docstring promises.
This unblocks the bypass design in PR #203
(`feature/dminterp-bypass-element-check`), which needs an authoritative
per-rank cell-id source. With the new default,
`mesh._get_closest_local_cells_internal(coords)` directly gives the cell
id whose closure contains each query — exactly what the bypass requires.
Underworld development team with AI support from Claude Code (claude.com/claude-code)
Member
Author
|
Updated precursor #207 design: loose is now the default in # in _function.pyx
cells = mesh._get_closest_local_cells_internal(coords) # authoritative, no kwarg neededNo |
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
When the caller of
DMInterpolationSetUp_UW(UW3's fork of stockDMInterpolationSetUp) supplies a per-query cell hint, PETSc re-verifies the hint viaDMPlexLocatePoint_InternalinsideDMLocatePoints. On cell shapes where its in-cell test is geometrically subtle — notably 2-D simplex cells embedded in 3-D space — verification can reject a correct hint and fall back to a wrong-cell match, yielding wild FE evaluations.This PR skips
DMLocatePointswhen:mesh.dm.isSimplex()), where the face-containment test is exact for the affine reference map, ormesh.dim != mesh.cdim), where the kdtree-nearest cell on a closed surface correctly contains every query.Non-simplex volume meshes (quads, hexes) keep the original
DMLocatePointspath, since deformed cells can have non-planar faces where the kdtree-nearest cell can be wrong and PETSc's more rigorous search is needed.What's in the patch
src/underworld3/function/_function.pyx— hint source inpetsc_interpolateisMesh.get_closest_cells(kdtree-nearest, no in-cell test). Comment explains the role.src/underworld3/function/_dminterp_wrapper.pyx::create_structure— gates whether the hint is passed to PETSc onmesh.dm.isSimplex() or mesh.dim != mesh.cdim. On non-simplex volume meshes the hint is suppressed (NULL) and PETSc'sDMLocatePointsruns unchanged.src/underworld3/function/petsc_tools.c::DMInterpolationSetUp_UW— whenowning_cellis provided, skipDMLocatePointsentirely. TheMPI_Allreduce(MIN, foundProcs)after each rank claims its hinted points assigns each global point to the lowest-rank claimant, preserving the prior ownership convention. Lazily buildspointVec/globalPointsScalaronly on the DMLocatePoints path.Net: ~60 / -40 lines across the three files.
Behaviour
DMLocatePoints(unchanged)SphericalManifold, when PPE lands)_evalf=TrueRBF workaroundTest plan
pytest tests/test_0820_deform_mesh_solver_rebuild_regression.py::test_stokes_velocity_updates_after_second_deform— passes (this was the regression in the previous push; quad mesh now goes through DMLocatePoints unchanged)pytest tests/test_0000_imports.py tests/test_0001_meshes.py tests/test_0004_pointwise_fns.py— 28/28uw.function.evaluateon quad and simplex meshes — analytic matchfeature/parallel-point-eval:probe_manual_p1_eval.pyonSphericalManifold— wild-evaluation rate 16% → 0feature/parallel-point-eval:probe_manifold_advection.py— SLCN advection runs without_evalf=TrueCopilot review feedback (addressed)
pointVec/globalPointsScalarconstruction into theelsebranch so the bypass path skips them.(PetscInt)owning_cell[p] >= 0toowning_cell[p] != (size_t)-1to avoid implementation-defined unsigned→signed conversion.Known limitations / out of scope
get_closest_cellsdoesn't carry an ownership filter. This is symmetric with the pre-existing behaviour for any single-rank UW3 evaluate (typical use). The general parallel-point-eval architecture is being addressed separately onfeature/parallel-point-eval.owning_cell[p]access forp >= n_localinredundantPoints=FALSEmode (where N fromMPI_Allgathervcan exceed each rank's hint-array length) is a separate latent issue, not touched by this PR.DMInterpolationSetCellHintsto the stockDMInterpolationAPI (the stock wrapper unconditionally drops caller hints). Useful for non-UW3 PETSc users but unrelated to this UW3 bypass.Underworld development team with AI support from Claude Code