Point locator: reject foreign points cheaply, locate once, and containment-check the cell hint (items 1-3 of #551, fixes the #432 class) - #556
Conversation
…imes The lost-point walk in _get_closest_local_cells_internal tried the 50 nearest cell centroids and only stopped when every lost point had been found. A point no local cell could own therefore paid all 50 rounds - 51 containment tests against 1 for an owned point, measured - and every point already located was re-tested on each of them, so one unfindable point charged the whole batch. In parallel the fraction of points a rank does not own is exactly what grows with rank count. Two bounds. A point inside a cell is no further from its nearest kd-tree control point than from that cell's centroid, which for a convex cell is within the cell's vertex reach; the largest local reach is now recorded with the kd-tree and a lost point beyond twice that distance is rejected before the walk starts. The factor of two is slack for the in-cell test's face tolerance and for badly shaped cells. And the working set shrinks: a point leaves as soon as a cell claims it, or as soon as the sorted neighbour distances pass the rejection radius. Measured on the #551 dossier probes: a foreign point costs 1.0 containment tests instead of 51.0 (0.53 us/point instead of 7.99 serial, 0.79 instead of 9.6 at np=4), and one unfindable point in a batch of 1000 adds 3 point-tests instead of 991 (2-D) or 12251 (3-D). The nearest containing centroid now wins, where before the winner was the LAST of up to 50 rounds and so depended on whether some unrelated point in the same batch was findable. Across a nine-set battery in 2-D and 3-D at np=1/2/4, 15 of ~450000 located cells change; every one is a point that both the old and the new cell contain. Addresses item 1 of #551. Underworld development team with AI support from Claude Code
Serial simplex meshes assert "exact" location capability, which makes the UW3 cell hint authoritative and skips PETSc's DMLocatePoints entirely - while the hint they handed over was get_closest_cells, a nearest-CONTROL-POINT kd-tree lookup with no containment test at all. On a tetrahedron nothing downstream can rescue that: the only remaining guard is a componentwise box clamp on the reference coordinates, and the reference tet is not the reference box. A query on a shared edge was answered by extrapolating the basis of a cell that does not contain it. That is #432, a recurrence of #390. The serial-simplex branch now takes _robust_owning_cells, the containment- checked locator every other authoritative path already used, so the three branches collapse to one call. Points it cannot place come back as -1, which the C bypass treats as a non-claim; they surface in unlocated_mask and take the RBF fallback that is already plumbed. The alternatives were worse. Restricting the "exact" assertion pushes serial simplex meshes back onto DMLocatePoints, which is slower and re-opens the #390 class of silent drops the bypass was added to close. A barycentric clamp in C only pins a wrong cell's reference coordinates to that cell's boundary, which gives the edge value of the WRONG cell - right only for continuous fields, and no help at all when the nominated cell is not adjacent to the point. Measured cost: one extra containment test per point, 1.0 tests per point for interior, on-face and on-edge queries. Measured benefit: 3-D P1 evaluation at quarter-points along cell edges was off by 3.0e-01 at 17 of 1318 points and is now exact. Addresses item 3 of #551 and closes #432. Underworld development team with AI support from Claude Code
points_in_domain located the points near the domain boundary, returned a boolean mask and threw the owning cells away; petsc_interpolate then located every interior point again. The near-boundary points were being searched for twice per evaluate call. Mesh._classify_points_in_domain returns both - the mask and the cells the classification actually looked up, with -1 meaning "not looked up" for an interior point and "not in the local mesh" for an exterior one. evaluate passes those to petsc_interpolate as cell_hints, which searches only for the entries still marked -1, and only on the DMInterpolation cache miss that needs them. A cache hit locates nothing at all, which the first version of this change got wrong: filling the whole hint array in the classifier made serial evaluate 29% slower because it searched on every call, cached or not. points_in_domain keeps its signature, its answer and its cost - it is now a one-line wrapper and does not search on the interpolator's behalf. Only the robust locator's answer is kept as a hint; the cell-wall test the serial classifier uses runs at a different face tolerance and its answer is a classification, not a hint. Measured on the #551 dossier probes: containment point-tests per located point inside evaluate fall from 2.28 to 1.10 at np=2, 3.35 to 1.20 at np=4 and 4.19 to 1.27 at np=8; evaluate wall time 0.0382 to 0.0361 s at np=4. Addresses item 2 of #551. Underworld development team with AI support from Claude Code
…on is O(1) #432 and #390 have both been fixed and returned, and there was no regression test pinning on-edge or on-face queries. tests/test_0761_point_locator.py adds them in 2-D and 3-D. The oracle is the closed-form P1 value (1-t)*u_a + t*u_b along a cell edge, computed WITHOUT uw.function.evaluate so it cannot inherit the defect it is testing for - the same discipline test_0753's arbitrary-coarse-field reference uses, and for the same reason. Covers shared vertices, edge midpoints, edge quarter-points and 3-D face centroids. The performance guard counts containment tests per point rather than seconds. The count is what the algorithm does; a wall time is what the machine was doing at the time, and the house has been burned by timing tests before. Parallel wall-clock numbers belong in the PR, not in an assert. Negative controls, because a test that cannot fail proves nothing. The nodal field is MEASURED to vary by more than 0.1 between a cell edge and its midpoint, so a wrong-cell interpolant cannot pass by being smooth (a linear field cannot test neighbour selection). get_closest_cells is shown to nominate cells that do not contain the query, so the edge test is pinning something real - and only at 3-D quarter-points, since a midpoint is never misassigned and 2-D never is, which is why #432 is a 3-D report; the 2-D case skips with that count in its message. The containment counter is asserted to fire, so the bounds cannot pass by instrumenting nothing. Validated fail-before with the fixes stashed and the tree rebuilt: 3-D P1 evaluation at t=0.25 off by 3.036e-01 at 17 of 1318 points; a point outside the mesh costing 51.0 containment tests; one unlocatable point adding 991 (2-D) and 12251 (3-D) point-tests to a batch of 1000. Underworld development team with AI support from Claude Code
There was a problem hiding this comment.
Pull request overview
This PR hardens and speeds up point location for uw.function.evaluate by (1) adding an O(1) rejection path and shrinking working set in the centroid-walk locator, (2) reusing point-classification location results as interpolation cell hints to avoid double location, and (3) ensuring any “authoritative” cell hint is containment-checked (fixing the #432 mis-evaluation class on serial simplex meshes).
Changes:
- Add
_local_cell_reachand use it to cheaply reject definitely-foreign points, while shrinking the centroid-walk working set as points are located. - Split
points_in_domaininto_classify_points_in_domainthat returns both the in/out mask and reusable owning-cell hints; thread hints intopetsc_interpolate. - Add a new
level_1/tier_aregression test suite covering containment-correct hints and bounded containment-test counts.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 2 comments.
| File | Description |
|---|---|
src/underworld3/discretisation/discretisation_mesh.py |
Adds local reach computation and improves _get_closest_local_cells_internal; introduces _classify_points_in_domain to return both mask + cell hints. |
src/underworld3/function/_function.pyx |
Threads classification-provided cell_hints into petsc_interpolate and enforces containment-checked authoritative hints via _robust_owning_cells. |
tests/test_0761_point_locator.py |
Adds regression tests for #432 and for bounded containment-test counts (rejection + shrinking working set), plus classifier/hint handoff checks. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| self._mark_local_boundary_faces_inside_and_out() | ||
|
|
||
| max_radius = self.get_max_radius() | ||
|
|
||
| if model_points.shape[0] == 0: | ||
| return numpy.array([], dtype=bool) | ||
| return (numpy.array([], dtype=bool), | ||
| numpy.array([], dtype=numpy.int64)) | ||
|
|
||
| cells = numpy.full(model_points.shape[0], -1, dtype=numpy.int64) | ||
|
|
||
| # Cd-1 surface mesh: no boundary-face control points exist | ||
| # (see _mark_local_boundary_faces_inside_and_out). Per the | ||
| # surface-mesh contract, query points are assumed to lie on | ||
| # the manifold; the closest-local-cell test is the right | ||
| # filter, not an inside/outside split. | ||
| if self.boundary_face_control_points_kdtree is None: | ||
| return self._get_closest_local_cells_internal(model_points) != -1 | ||
| in_or_not = self._get_closest_local_cells_internal(model_points) != -1 | ||
| return in_or_not, cells |
| for face in range(fStart, min(fEnd, fStart + 1500)): | ||
| closure = mesh.dm.getTransitiveClosure(face)[0][-3:] | ||
| corners.append(verts[closure - pStart]) | ||
| corners = np.array(corners, dtype=np.float64) # (n, 3 vertices, 3 coords) |
Adversarial review — PR #556, point locator (items 1-3 of #551, #432 class)Reviewed at We are independent of the implementer. Verdict: request changes. One merge-blocker, with a one-line fix we have MERGE-BLOCKERSB1.
|
| build | interior points classified in-domain | locator returns -1 | NaN out of evaluate |
|---|---|---|---|
merge base e475246f |
2992 | 1 | 0 |
PR head 82145992 |
2992 | 1 | 1 |
The point is at (0.5254, 0.2065, 0.3302), points_in_domain says it is
inside, and V.rbf_interpolate at that exact coordinate returns 0.2802 — the
value the rung was supposed to write.
This is a silent NaN in the return value of the single most-used function in
the library, on the most common configuration (serial simplex), reachable on
any mesh where the 50-neighbour walk can miss — i.e. deformed, graded and
adapted meshes, which is the free-surface / adaptivity / fault workload.
Fix (measured):
vars = list(mesh.vars.values()) # _function.pyx:1273With that one line applied and rebuilt, the same probe returns 0 NaN and
the point takes the RBF value. Nothing else changed.
Two riders on the fix, both of which belong in this PR because it is this PR
that makes the rung load-bearing:
- The same exhausted generator is read again at
_function.pyx:1308,
all_continuous = all(getattr(var, "continuous", True) for var in vars).
all()over an empty generator isTrue, so the continuity gate on the
"continuous"location capability has never bound. Fixingvarsmakes it
bind — correctly, but for the first time, on cubed-sphere-class warped-hex
meshes. On"exact"meshes it is a no-op either way
(_hint_is_authoritative(True) == _hint_is_authoritative(False) == True,
measured on a quad box carrying a P1 and a P0 variable), so the exposure is
bounded, but it needs its own line of measurement before it lands. - Add a test. The natural one is the probe above: a graded serial 3-D mesh, a
random interior batch, assertnp.isfinite(evaluate(...)).all(). It fails at
PR head and passes with the fix.
Note the pre-existing half of this: in parallel the rung was already
reachable before the PR, because parallel already used _robust_owning_cells.
Same deformed 3-D mesh, node-point queries:
| np | NaN at merge base | NaN at PR head |
|---|---|---|
| 1 | 0 | 0 |
| 2 | 18 | 18 |
| 4 | 120 | 120 |
So the dead rung is not this PR's bug. Opening the serial door onto it is.
Major
M1. The answer-identity claim is measured on the point sets that bypass the walk
"15 of roughly 450 000 located cells differ ... nothing changed in 2-D, nothing
at tol = 0, and nothing for vertices, edges, faces, centroids, boundary or
exterior points."
We reproduced the 15 exactly, from the artifacts, independently of their
harness (before/ vs after/ cells_d{2,3}_np{1,2,4}.json, every set, every
rank, every tolerance): 15 differences over 661 962 compared entries, all 3-D,
all at _EVAL_FACE_TOL, all in scatter/mixed, 0 found/lost flips, and
only 4 distinct query points (indices 1255, 2439, 3122, 3652) appearing
across ranks and batches. The claim is true as stated.
It is also not the interesting population. Exact mesh vertices, edge midpoints
and face centroids are answered by the first-pass control-point test
(control points are 0.99*vertex + 0.01*centroid), so they never enter the
walk and cannot change. The points that do enter the walk are the ones sitting
near a shared vertex/edge/face. On a plain uniform 3-D simplex box
(cellSize=1/8, the PR's own mesh class), 2615 near-vertex barycentric queries:
| set | tol | cells differ (new vs old walk) | not double-contained |
|---|---|---|---|
| near-vertex barycentric | 0.0 | 922 / 2615 (35%) | 0 |
| near-vertex barycentric | _EVAL_FACE_TOL |
411 / 2615 (16%) | 0 |
| uniform random interior | _EVAL_FACE_TOL |
1 / 3000 | 0 |
and on a deformed (sliver / graded) mesh even uniform random interior points
move (338 / 3000 on the sliver mesh at the eval tolerance). Every single
difference we found, on every mesh, is a point both cells contain under the
same test — the contract holds, and the new rule (nearest containing centroid)
is strictly better than the old one (last of up to 50 rounds, dependent on
whether an unrelated point in the batch was findable). We are not disputing the
change; we are disputing the size claimed for it.
It matters because the cell choice is the answer for a discontinuous field.
Measured on that same box with a P0 (continuous=False) variable, values drawn
uniform on [-1, 1]:
uw.function.evaluatereturns the new cell's value for 2615/2615 points
and the old cell's for 2208/2615 — the located cell fully determines the
answer, as expected;- over the 407 points whose cell changed,
|P0(new) - P0(old)|has
max 1.935, mean 0.646 on a field of range 2.
P2/P0-discontinuous is the pressure space the fault work runs on. The fault set
(0845-0848) and 1018 pass (70 passed, 167 s), so nothing is broken — but
"15 of 450 000" reads as "no downstream consumer can notice", and a consumer
evaluating a discontinuous field at node/face coordinates can notice. Please
restate the identity result as what it is: the walk population moves, and it
moves by O(jump) for discontinuous fields; the sets measured were mostly not
the walk population.
For the record on the rest of the downstream question: swarm cell ownership
(swarm.py:5145, ddt.py:2401), adaptivity.py:752 and
gradient_evaluation.py:240 all use get_closest_cells, which the PR does not
touch; _dminterpolation_cache keys on an xxhash of the coordinate bytes plus
dofcount plus policy, so a cache hit implies identical coordinates and cannot
serve a stale hint; points_in_domain's returned mask is bit-identical by
construction (the new cells[...] assignments are the only additions, and
in_or_not[...] = cells[...] >= 0 is the same expression). We could not find a
consumer that depends on the old tie-break other than evaluation of
discontinuous fields.
M2. The dangerous direction of the rejection radius is untested, and item 2's contract assertion never fires in serial
tests/test_0761_point_locator.py pins "a point outside costs O(1) tests" and
"every returned cell contains its point". Neither is the failure mode that
would be silent: a point that IS in a local cell being rejected. Nothing in
the file would catch a wrong reach.
That gap is worth closing because the guard is cheap and the harness works. Our
negative control (brute-force oracle, 3-D box, 2000 interior points, sweeping
Mesh._LOCATOR_REACH_MARGIN):
| margin | false rejections |
|---|---|
| 2.0 (shipped) | 0 |
| 1.0 | 0 |
| 0.5 | 9 |
| 0.25 | 469 |
| 0.1 | 607 |
So the bound is tight at ~1.0 and the shipped 2.0 has a factor-2 margin over
the first observable failure — good, and exactly the kind of number that should
be in the test file rather than in a PR table.
Second, test_the_classifier_hands_over_the_cells_it_located ends with
offered = in_or_not & (cells >= 0)
if offered.any():
assert mesh._test_if_points_in_cells_internal(...).all()In serial _eval_use_robust_location() is False, so the classifier never
fills a hint and offered is always empty: measured 0 hints for 1465
interior points in serial, against 462/462 at np4. The strongest assertion in
item 2's only test is skipped in the default (serial) test run. The file
applies the "confirm the probe fires" discipline elsewhere
(assert per_owned_point >= 1.0, with a comment explaining why) — apply it
here: assert the serial expectation explicitly (offered.sum() == 0 in serial)
and require offered.any() under MPI.
Minor
m1. Reach invalidation is asserted in prose only
The comment says _local_cell_reach is "invalidated with _index by deform
and adapt". It is, and we checked it: a 2-D box deformed by ×50 gives
0.0706932 -> 3.53466 (ratio exactly 50.0), deformed back by /50 returns to
0.0706932, and 2000 queries in the expanded domain are located with 0 false
rejections. A stale small reach is the one thing that silently breaks the
rejection, so pin it with a two-line test.
m2. cell_hints is mutated in place
_function.pyx:1354, cells = np.ascontiguousarray(cell_hints, dtype=np.int64)
returns the caller's array when it is already int64 and contiguous (verified),
and cells[unhinted] = ... at :1362 then writes into it. Today the only
caller passes cell_hints[in_or_not], a fresh fancy-index copy, so nothing is
observable — but petsc_interpolate takes cell_hints as a documented
keyword. Use np.array(cell_hints, dtype=np.int64, copy=True).
m3. _build_kd_tree_index_PIC would leave a stale reach
It sets self._index (and a differently-built _centroid_index) but not
_local_cell_reach, so if it ever ran after a normal build it would leave the
old reach in place against new geometry — the one stale-value scenario the
getattr(..., None) guard does not cover. The PR already notes it has no
callers; we confirmed (one hit in the whole tree, the def itself). Delete it
or set the reach — do not leave a second index builder that skips the new
invariant.
m4. The face tolerance is absolute, the rejection radius is mesh-relative
_mark_faces_inside_and_out places the inner/outer control points at a fixed
±1e-3 in physical units, so the loose test admits a slab of
tol * 1e-3 ~ 1e-5 absolute, while the rejection radius is
2 * max cell reach. They cross over when the largest local cell has reach
below ~5e-6, at which point the radius can reject something the containment
test accepts. Not reachable in any realistic model — we built a mesh with the
whole domain 1e-4 across (reach 6.9e-6) and got 0 false rejections — and the
containment test is already meaningless at that scale for its own reasons. But
the two scales should be tied to each other rather than to each other's
accident.
Pre-existing defects surfaced by this review (separate issues, not blockers)
-
The 50-neighbour cap loses genuinely interior points. Attribution probe
(a Python transcription of the pre-Point locator: reject foreign points cheaply, locate once, and containment-check the cell hint (items 1-3 of #551, fixes the #432 class) #556 walk, run against the shipped one on
the same mesh and points, with the brute-force oracle as referee):
newly-lost = 0 in every case we ran — the PR loses nothing the old code
found. But both lose:mesh set tol oracle-contained lost (old) lost (new) graded 3-D ( x**4)uniform interior _EVAL_FACE_TOL3000 6 6 graded 3-D near-vertex 0.0 2615 29 29 sliver 2-D (y/200) uniform interior 0.0 3000 2 2 uniform 3-D box largest cell's vertices 0.0 24 2 2 This is the population that hits B1. Worth an issue of its own: the walk
should widenk(or walk the cell adjacency instead of the centroid kd-tree)
rather than cap at 50 nearest centroids. -
points_in_domainstill callsget_max_radius(), an@collective_operation
allgather, inside the routineevaluateuses to decide which ranks hold
points. The PR flags it and defers to item 4. Agreed.
What we attacked and could not break
- The rejection bound is sound by construction and holds in practice. Every
cell contributes its centroid to the control-point kd-tree and the reach is
the max over local cells of|vertex - centroid|, so for any pointpin
cellc,dist(p, nearest control point) <= dist(p, centroid_c) <= reach_c <= max reach. That holds for P1 simplices and for bi/trilinear quads and
hexes (their maps are convex combinations of the vertices), independently of
grading, aspect ratio or how farpis fromcentroid_c. Confirmed against
a brute-force oracle (every query point tested against every local cell with
the same containment test) with 0 newly-lost points on: uniform 2-D/3-D
simplex, quad 2-D, hex 3-D, annulus, spherical shell, graded 2-D/3-D
(x**4, grading 6.2-7.4), sliver 2-D (aspect 200), a 1e-4-wide mesh, a
6371-wide mesh, and all of the above at np1/np2/np4.sqr_dists=False
genuinely returns distances, not squares (checked), so the comparison is
dimensionally right. - Parallel, graded partition. 2-D np4 gives per-rank reaches
[0.1015, 0.0286, 0.1040, 0.0267](spread 3.9x); 3-D np4
[0.382, 0.385, 0.164, 0.387](2.4x). Across the global node set:
0 points claimed by 0 ranks in every configuration, mask identical to
points_in_domainon every rank, 0 hints that do not contain their point,
0 hints offered for exterior points. - Item 2 / the cache. Cache key is
(xxhash(coords), dofcount, policy), so
a hit implies byte-identical coordinates; the hint is only computed on a
miss; a stale hint cannot be reused. Splitting the location into two batches
(classifier + lazy fill) is only answer-safe because of item 1 — the old
walk's answer depended on the rest of the batch. Worth stating in the commit
message: commit 3 depends on commit 1 for correctness, not just for speed. - Item 3 / no performance cliff. Containment tests per point after the
collapse to one branch, serial: simplexcellSize=1/601.004 (interior) /
1.000 (nodes); quad 60x60 1.000 / 1.000. The claimed 1.0 holds on both the
paths that previously had their own fast branch. Outside points, measured
on our own meshes at PR head: 2-D 1.000, 3-D 1.000 tests per point, against
the 51.0 the merge-base build produces on the same query set. - The multi-mesh hint guard (
mesh is hinted_mesh) is unreachable today.
functions_unit_system._evaluate_implrejects an expression spanning two
meshes beforepetsc_interpolatesees it ("Expression contains MeshVariable
symbols from 2 different meshes"), so the guard is defensive-only and
untested. Correct to have; worth knowing it is not exercised. - Fail-before, independently reproduced by rebuilding
src/at the merge
base with the PR's tests in place: 7 failed, 11 passed, 2 skipped, including
3-D P1 evaluation at t=0.25 along cell edges is off by 3.036e-01 at 17 of 1318 points(uw.function.evaluate returns wrong values for P1 fields at points exactly on cell edges (3D) #432) anda point outside the 2-D mesh costs 51.0 containment tests per point.
Test results at PR head
tests/test_0761_point_locator.py: 18 passed 2 skipped (serial), 18/2 at
np2, 18/2 at np4.tests/test_0845*-0848*+tests/test_1018_rotated_freeslip.py:
70 passed, 167 s.- Full gate
pytest tests -m "level_1 and tier_a" -q -p no:cacheprovider --ignore=tests/test_0050_utils.py: 622 passed, 19 skipped, 1545
deselected, 1 xfailed, 0 failed in 396 s — reproduces the PR's claim
exactly. scripts/check_deprecated_patterns.py: clean, allowlist unchanged at 79.
Note what the gate does not cover, and why B1 got through it: every
level_1 / tier_a evaluation test runs on a quasi-uniform, undeformed mesh,
where _robust_owning_cells never returns -1 for an interior point, so the
dead RBF rung is never asked to run.
Artifacts
Probes written for this review (scratchpad):
oracle.py (brute-force containment oracle, 11 mesh cases),
attribute.py (pre-#556 walk transcribed and run side by side),
negctl_deform.py (reach-margin negative control + deform invalidation),
downstream.py (discontinuous-field consequence, serial fast-path counts),
nanprobe.py / nanprobe2.py / parnan.py (B1),
parallel_seam.py (graded-partition parallel contract).
The RBF fallback in petsc_interpolate has never run. It iterates mesh.vars.values(), and mesh.vars is a weakref.WeakValueDictionary whose .values() is a GENERATOR, not a view: the dofcount loop a few lines earlier consumes it, so the fallback loop iterates nothing and the NaN that DMInterpolationEvaluate_UW writes for an unplaced point is returned to the caller. That did not matter until the previous commit in this branch, which routed serial simplex evaluation through the containment-checked locator. That locator returns -1 - the older nearest-control-point hint never did - so the dead rung became load-bearing on the most common configuration in the library. Measured on a graded 3-D simplex box (cellSize 1/8, deformed x -> x**4), 3000 interior queries: the merge base returns 0 NaN, this branch returned 1. Materialising the list fixes it: same probe, 0 NaN, and the point takes the RBF value (0.2802) the rung was always supposed to write. The same exhausted generator was read again by the continuity gate, so all() over it was vacuously True and the gate has never bound. Un-breaking it is a behaviour change, so it is measured rather than assumed. The gate only does anything on meshes whose measured location capability is "continuous" - warped hexes, the cubed-sphere class; simplex, quad, annulus and rectilinear hex boxes are all "exact", where the gate is a no-op either way. On a warped hex box carrying a P1 and a P0, binding it over EVERY variable on the mesh would cost the CONTINUOUS field a factor of 15 in accuracy (linear field, max error 8.0e-3 -> 1.2e-1 at 62 of 1500 interior points) because one discontinuous variable elsewhere took every evaluation off the authoritative path. Scoped instead to the variables the call actually asks for - which is what the surrounding comment always said the policy was - the continuous field is bit-identical to before and only the P0 moves, by up to 1.71 on a field of range 2, which is the O(jump) correction the gate exists to make. Also here, from the same review: - petsc_interpolate copies cell_hints before filling its -1 entries. np.ascontiguousarray hands back the caller's own array when it is already int64 and contiguous, and cell_hints is a documented keyword, so the fill was writing through somebody else's array. - _build_kd_tree_index_PIC and _build_kd_tree_index_DS are deleted. Both set _index without the new _local_cell_reach, which is the one way a stale rejection radius could outlive the geometry it was measured on; both have had zero callers for a long time and are already on the readability review's delete list (READ-31). - The tie-break change is written down where it can be found: which containing cell a shared vertex/edge/face query returns has moved, and for a discontinuous field the cell is the answer. - The absolute face-tolerance slab and the mesh-relative rejection radius are noted as the different scales they are, with the domain size at which they would cross (about 5e-6 across, measured clean at 1e-4 and at 6371). Underworld development team with AI support from Claude Code
Nothing asserted that a point a local cell really does contain is
never rejected by the new radius, which is the failure mode that
would not announce itself: the point vanishes into the RBF fallback
in serial, or is claimed by nobody in parallel.
Four tests, all with their own negative control.
A brute-force oracle - every query against every local cell, with the
same containment predicate the locator uses - says which points are
genuinely owned, and the reach margin is swept to prove the probe
fires. Serial: 2-D {2.0: 0, 1.0: 0, 0.5: 5, 0.25: 16}, 3-D
{2.0: 0, 1.0: 0, 0.5: 4, 0.25: 168} out of 800 interior points. The
bound is tight at about 1.0, so the shipped 2.0 carries a factor of
two over the first observable loss. 0 at the shipped margin on every
rank at np2 and np4 as well.
A P0 discontinuous field evaluated at shared mesh vertices pins the
tie-break: the value must belong to one of the cells that contains
the query, not a specific cell (that was never the contract, and
pinning it would break on the next legitimate change). Every vertex
in the set is contained by 18 to 44 cells and their values differ by
more than 0.5, so a wrong cell is visible.
The RBF fallback rung gets a fault-injection test: five points are
forced to -1 and the result must be finite AND equal the RBF value
AND leave the other 395 points bit-identical. Waiting for a graded
mesh to lose a point naturally works, but how many it loses is a
property of whatever gmsh produced that day. Confirmed fail-before by
rebuilding with the one-line generator fix reverted: exactly the five
injected points come back NaN.
The reach invalidation was asserted in prose only. Deform a 2-D box by
50 and the stored reach must follow exactly, and every point in the
expanded domain must still be located.
Finally, test_the_classifier_hands_over_the_cells_it_located ended
with its real assertion inside `if offered.any():`, and offered is
ALWAYS empty in serial - 0 hints for 1465 interior points, against
462 of 462 at np4 - so the strongest claim in item 2's only test was
a no-op in the default test run. It now asserts the serial reality
(no hints, by design) and requires hints under MPI.
Underworld development team with AI support from Claude Code
PR #556 failed CI on test_rotated_workspace_deform_invalidates with "post-deform solve differs from fresh control by 1.16e-01", and the failure does not reproduce here. Not in the amr-dev toolchain (custom AMR PETSc), not in the dev toolchain CI itself uses (conda PETSc), not as a single test, not as the whole file, not in the exact CI batch shape (pytest tests/test_101*py tests/test_102*py: 137 passed in both environments), with or without CI's environment variables. Locally err is not small, it is exactly 0.0: one Krylov iteration, |r| = 8.267546030330813e-12 bit-identical between the post-deform solve and the fresh control. The locator work in this PR is measurably inert in this test: - the mesh is a StructuredQuadBox, and for a quad/hex mesh the merge base ALREADY located through _robust_owning_cells - item 3 changed the serial SIMPLEX path only, so the evaluate route here is unchanged; - there are 8 point-location calls during the test (from swarm.migrate -> points_in_domain and from global_evaluate, not from the assembly). Running every one of them twice, with the new rejection radius and with it disabled, gives identical cells: 0 differences out of 2270 points; - the deformed quad box measures location capability "exact" before and after the deform, so _hint_is_authoritative is True whatever the field continuity and the continuity gate cannot bind; - in serial the classifier hands over no cell hints at all, so item 2 is inert too; - the locator reach follows the deform (0.070710678 -> 0.073041055). Rather than guess at a fix for something that cannot be reproduced, the test now distinguishes the two things it could be. Both solvers assemble the same system on the same deformed mesh from a zero guess, so any deterministic solver owes them the same answer: if both report converged and the answers still differ, the OPERATORS differ and something survived the deform; if the convergence reports differ, the linear solves did. The failure message carries both solves' reason, iteration counts and residuals plus the locator reach before and after, and the reach invalidation is now asserted here as well as in test_0761 (the review's minor m1) because this is the only test in the suite that deforms a mesh between two solves. Underworld development team with AI support from Claude Code
|
Response commits B1 fixed — The
Binding it mesh-wide would cost the continuous field a factor 15 in accuracy for no correctness gain — one discontinuous variable elsewhere de-authorises everything. Scoped to the requested variables it is what the surrounding comment always claimed, leaves continuous evaluation bit-identical, and moves only the discontinuous field, which is the O(jump) correction the gate exists to make. M2 — M1 — documented rather than buried: the restated identity claim is in the PR body, a Minors: m1 reach-invalidation now asserted in tests (ratio exactly 50.0); m2 The CI failure on Gate 627 passed / 0 failed (622 before, +5 new). Fault set + 1018: 70 passed. Attribution oracle against the pre-#556 walk across 10 mesh cases: newly-lost 0 everywhere. Underworld development team with AI support from Claude Code |
… the difference CI's first diagnosis-grade report on test_rotated_workspace_deform_invalidates said both solves converged from an IDENTICAL initial residual (|r0| = 0.02785134866629431 to every digit, |r| ~ 1e-11 for both) and still disagreed by 7.02e-02. Two solves that both drive the residual to machine zero on the same system can only differ in the null space, so this round measures the null space instead of arguing about it. Solver side: the rigid-rotation gauge decision was invisible. Whether a mode is admitted decides whether a component is projected out of the answer, and _mode_satisfies_constraints made that call silently, per mode, from the boundary normals. It now optionally records the constraint violation, the operator violation and the verdict; _finalize_rotated_solution collects one record per offered mode and the solve result carries it as "rotation_gauge". Default arguments unchanged, so nothing else moves. Test side, all of it reported in the failure message and all of it computed EAGERLY, because an instrument that only runs when the test fails is an instrument that has never been run: - the difference field is decomposed onto the rigid-body span (the same modes the solver considers, plus translations, Gram-Schmidt'd in the same order, built from nodal coordinates so it does not inherit the machinery it is measuring). If the difference lives in the span, one solve admitted a mode the other rejected and this is a #543 gauge bug; if it lives off the span, the two operators differ and something survived the deform; - per-solver constrained-row count, distinct-row count and boundary list; - a point-location tally around the deform, the post-deform solve and the control, answering every call TWICE - as shipped and with the rejection radius set aside - and reporting how many answers the radius changed. That is the number that decides whether #556 is implicated, taken on whatever mesh the machine actually built rather than argued from a local run. Three negative controls, so none of those numbers is unchecked: a pure rigid rotation must decompose with 4e-18 of it off the span and a random field with 0.9985 of it off (both asserted on the mesh under test), and squeezing the reach margin to 0.05 must make the radius comparison see answers change. Measured here: 84 constrained rows for both solvers, one rotation mode offered and rejected by both with identical violation 3.045e-01, gauge_removed False for both, and during the deform 8 location calls over 2270 points with 38 returning -1 and radius_changed = 0. The locator changes nothing on this platform. If CI reports otherwise, that is the finding. Underworld development team with AI support from Claude Code
…es them The test made two claims in one assertion. One is well posed and is what #543 wrote it for: the deform invalidated the workspace. The other is not: that a post-deform solve matches a fresh control to 1e-6, on a system that does not determine its own answer in one direction (#560). They are now separate. (a) stays hard: rotation_reused and workspace_reused both False, the locator's rejection radius followed the deform, both solves converged. (b) is narrowed, not relaxed. Rotated free-slip on a curved boundary loses the constant-pressure gauge, and the solution acquires a component along one unpinned direction whose amplitude is round-off. Measured: a coordinate change of two machine epsilons (4.44e-16) moves the velocity by 1.33e-01, and the move does not scale with the perturbation - 4.4e-16, 2.2e-15, 1e-14, 1e-12 and 1e-9 all give between 4e-2 and 2e-1. That is why this assertion was green on macOS (err exactly 0.0 in 81 consecutive runs across two PETSc toolchains and nine PYTHONHASHSEEDs) and intermittently red on CI: it passes only where the two assemblies agree bitwise. The unpinned subspace is exactly one-dimensional - five different perturbations move the answer along the same direction to cosine 1.000000, the normalised difference set has singular values [2.236, 4.4e-9, 3.4e-9, 2.3e-9, 1.7e-9], and removing the leading direction leaves 2e-9 of each difference. So the test measures that direction with one extra perturbed solve and requires the two solutions to agree in every OTHER direction, at the same 1e-6 it always used. The tolerance is untouched; it is the claim that is made honest. Two things keep it from becoming a rubber stamp. A negative control injects a 1e-3 discrepancy orthogonal to the unpinned direction and asserts it survives the projection (measured 1.000e-03 against a 1e-4 floor), so the projection cannot absorb a real disagreement. And the test branches on what it measures: if a 2-eps perturbation stops moving the answer - i.e. when #560 is fixed - the projection becomes a no-op and the solutions are compared directly again, with the branch reported in the failure message. Fixing #560 strengthens this test instead of breaking it. All the instrumentation stays: the gauge decisions, the constrained-row counts, the locator tallies and the rigid-body decomposition are what turned an unreadable CI failure into a filed defect, and they are the diagnostic for the next one. Verified: this test 10/10 in amr-dev and 10/10 in dev, the CI batch shape (tests/test_101*py tests/test_102*py) 137 passed in dev, full level_1/tier_a gate 627 passed 0 failed. Underworld development team with AI support from Claude Code
The point locator: reject early, shrink the working set, and check the hint contains the point
Addresses items 1, 2 and 3 of #551. Items 4 (rank-independent collective structure),
5 (pointwise derivative evaluation) and 6 (the cheap follow-ups) are separate PRs.
Three changes, one file each side of the JIT boundary.
1. The walk gets a rejection path and a shrinking working set
Mesh._get_closest_local_cells_internalhandles a point the first containment testrejects by walking the 50 nearest cell centroids, and the only way out of that loop was
every lost point being found. A point no local cell could own therefore paid all 50
rounds, and every point already located was re-tested on each of them.
Two things now bound it.
A rejection radius. Every cell contributes its centroid to the control-point
kd-tree, so a point inside cell c is no further from its nearest control point than
from c's centroid, which for a convex cell is within that cell's vertex reach. The
largest local reach is recorded when the kd-tree is built (
Mesh._local_cell_reach,invalidated with
_indexbydeformandadapt) and a lost point beyond twice thatdistance is rejected before the walk starts. The factor of two is deliberate slack: the
in-cell test admits a thin slab outside each face, and a badly shaped cell expands
further under that slab than a well-shaped one. Two reaches is still about one cell.
A shrinking working set. A point leaves the walk as soon as a cell claims it, or as
soon as the (sorted) neighbour distances pass the rejection radius. The nearest
containing centroid wins.
2. Classification and interpolation share one location pass
points_in_domainlocated the near-boundary points, returned a boolean mask and threwthe owning cells away;
petsc_interpolatethen located every interior point again.Mesh._classify_points_in_domainreturns both — the mask and the cells theclassification actually looked up.
evaluatepasses those topetsc_interpolateascell_hints, which locates only the entries still marked-1, and only on theDMInterpolation cache miss that needs them.
points_in_domainkeeps its signature, itsanswer and its cost: it is now a one-line wrapper and does not search on the
interpolator's behalf.
3. The cell hint is containment-checked (#432)
Serial simplex meshes assert
"exact"location capability, which makes the UW3 hintauthoritative and bypasses
DMLocatePoints— while the hint they handed over wasget_closest_cells, a nearest-control-point kd-tree lookup with no containment testat all. On a tetrahedron nothing downstream can rescue that: the reference-coordinate
guard in
petsc_tools.cis a componentwise box clamp, and the reference tet is not thereference box. A query on a shared edge was answered by extrapolating the basis of a
cell that does not contain it.
Fixed by option (a) of the three the issue lists: containment-check the hint. The
serial-simplex branch now takes
_robust_owning_cells, the same containment-checkedlocator every other authoritative path already used, so the three branches collapse to
one call. Points it returns
-1for surface inunlocated_maskand take the RBFfallback that is already plumbed.
Why not the other two:
"exact"assertion would push serial simplex meshes back ontoDMLocatePoints, which is slower and re-opens the test_1052 VEP stability tests fail under PETSc 3.25 (variable-dt pure-VE blow-up), pass in level_1 gate #390 class of silent drops thebypass was added to close.
that cell's boundary. On a shared edge that returns the edge value of the wrong
cell, which is right only for continuous fields, and it does nothing when the
nominated cell is not adjacent to the point at all.
Option (a) costs one extra containment test per point on the fast path — measured at
1.0 tests per point for interior, on-face and on-edge queries — and it makes the word
"authoritative" true.
Measurements
2-D unstructured simplex box,
amr-dev, Darwin arm64, runs sequential. Probes are theones from the #551 evidence dossier, unchanged, plus a cell-identity probe that dumps
the located cell for every point of a battery.
The headline: does the rejection path recover the anti-scaling? No.
global_evaluateon the fixed 20 000-point global query set,cellSize = 1/100, bestof 3, max over ranks:
Flat. The issue's expectation that item 1 "alone should recover most of §1" is not
confirmed, and the reason is worth recording:
points_in_domainalready had its owncheap rejection — the boundary-sign /
far_from_domaintest — so the migrate claim loopwas never paying the 51× penalty in the first place. Phase attribution after the change
(3 calls, rank 0):
Swarm.migratepoints_in_domainagainst 31 / 64 / 69 % for migrate before. The residual is exactly the signature the
dossier named: the migrate claim loop classifies the whole global point set on every
rank on every round, and the number of rounds grows with rank count. That is item 4
(collective structure), not item 1 — the cost is the kd-tree query plus numpy over
20 000 points per rank per round, with the walk now contributing almost nothing.
Containment work inside those same calls did drop (rank 0, 3 calls, point-tests):
global_evaluatebeforeglobal_evaluateafterWhere the 51× actually was: a direct locator call
_robust_owning_cellson 4 000 points,cellSize = 1/100:Foreign / owned wall-clock ratio: 9.2× → 0.7× at np=1, 13.1–13.8× → 1.1–1.2× at np=4.
evaluateon the fixed global set (rank-local classification) improves in step, both incontainment work and in time:
evaluate(s), before → afterContainment tests per point, by query class
Cell-identity probe,
tol = _EVAL_FACE_TOL; 2-DcellSize = 1/40, 3-DcellSize = 1/12,4 000 points per set (before → after):
The 3-D residual (8.9–11.2 for random interior at np=4) is real work, not the old
pathology: those are points near the partition seam whose nearest control point is
inside the rejection radius, so the walk legitimately runs — and it stops as soon as the
sorted neighbour distances pass the radius, which is what turns 42 into 10.
Answers
The cell-identity probe dumps the located cell for every point of nine query sets
(vertices, edge midpoints, 3-D face centroids, cell centroids, random interior, on the
domain boundary, outside the domain, other ranks' nodes, and a batch poisoned with one
unfindable point) at two tolerances, in 2-D and 3-D, at np = 1, 2 and 4.
15 of roughly 450 000 located cells differ. Every one is 3-D, at the loose
evaluation tolerance, in the random-interior set, and every one is a point that both
the old and the new cell contain (verified point by point: 15 checked, 0 failures).
These are points sitting within the containment test's ~1e-5 face slab, i.e. effectively
on a shared face, where more than one cell qualifies. Which one was returned was
previously decided by whether some unrelated point in the same batch happened to be
findable — the walk kept going and a later round overwrote the earlier answer. It is now
the nearest containing centroid, deterministically. Nothing changed in 2-D, nothing at
tol = 0, and nothing for vertices, edges, faces, centroids, boundary or exteriorpoints.
Tests
tests/test_0761_point_locator.py, 2-D and 3-D,level_1/tier_a, 18 passed2 skipped.
from a high-frequency signal of position, so the reference
(1-t)u_a + t u_biscomputed without
uw.function.evaluate— the same discipline astest_0753_nested_mg_prolongation.py::test_reproduces_an_arbitrary_coarse_field, andfor the same reason. Covers
t = 0(shared vertex),t = 0.5(edge midpoint),t = 0.25(quarter point) and 3-D face centroids.3-D P1 evaluation at t=0.25 along cell edges is off by 3.036e-01 at 17 of 1318 points. That is uw.function.evaluate returns wrong values for P1 fields at points exactly on cell edges (3D) #432. Also failing before:a point outside the 2-D mesh costs 51.0 containment tests per point, andone unlocatable point added 991 (2-D) / 12251 (3-D) containment point-tests to a batch of 1000.between an edge and its midpoint, so a wrong-cell interpolant cannot pass by being
smooth. (ii)
get_closest_cellsis shown to nominate cells that do not contain thequery, so the edge test is pinning something real — 3-D quarter points, where a
midpoint is (measured) never misassigned and 2-D never is either, which is why uw.function.evaluate returns wrong values for P1 fields at points exactly on cell edges (3D) #432 is
a 3-D report; the 2-D case skips with that count in the message. (iii) The containment
counter is asserted to fire (an owned point must cost at least one test), so the
bounds above cannot pass by instrumenting nothing.
point rather than seconds — the count is what the algorithm does, a wall time is what
the machine was doing at the time. The parallel timings live in this PR, not in a
brittle assert.
tolerances, on random interior points. That is the property the shrinking working set
must preserve, and it is the one that is actually defined — which qualifying cell is
returned for a point on a shared face never was.
_classify_points_in_domainis asserted to agree withpoints_in_domainexactly, to offer no hint for exterior points, and to offer only hints that contain
their point.
Also run green:
tests/test_1018_rotated_freeslip.py(22), the evaluate/locator set(
0503,0503_evaluate2,0506,0507,0730,0753×2,0755–0757,0760,0820_in_cell_test_loose_semantics,0057_deformed_domain_membership— 116), and thefault set
0845–0848(48).Full gate
pytest tests -m "level_1 and tier_a" -q -p no:cacheprovider --ignore=tests/test_0050_utils.py: 622 passed, 19 skipped, 1 xfailed, 0 failed(1545 deselected), 409 s.
scripts/check_deprecated_patterns.py: clean, allowlistunchanged at 79.
What this PR does not do
global_evaluateis still there. It is the migrateclaim loop classifying the global point set on every rank on every round (evaluate / global_evaluate: one pass over the locator, not five patches (anti-scaling to np=8, a smoothed substitution, an uncontained cell hint) #551 item 4),
not the locator walk. This PR removes the amplifier and measures the residual; it does
not fix it.
evaluatesubstituting a smoothed degree-1 projection for any expressioncontaining a derivative, uw.function.evaluate returns negative values for a squared DERIVATIVE expression #491) is untouched.
UnitAwareArrayguard (Enforce the ND<->units boundary: evaluate/global_evaluate should reject or nondimensionalise UnitAwareArray input #279), and the O(cells)
_mark_faces_inside_and_outrebuild (0.73 s at 23 264cells, redone on every
deform) — is untouched. The face-marking rebuild is now thelargest single fixed cost in a cold evaluation after a mesh move.
Two things noticed in passing and deliberately not fixed here (Charter §9):
points_in_domaincallsget_max_radius(), which is@uw.collective_operationanddoes an allgather. Every current caller reaches it on all ranks, so it is not firing
today — but it is a collective inside the routine
evaluateuses to decide whichranks hold interior points, which is precisely the petsc_interpolate: the DMLocatePoints call count is rank-local (cache key + unreduced location capability) — hangs at np>=4 #314 shape. It belongs with item 4,
and the surrounding lines are being edited by the Empty-rank support, layer 2: evaluate/points_in_domain, radii/centroid reductions, gather_data NaN-stripping #405 work.
Mesh._build_kd_tree_index_PIChas no callers anywhere in the tree.Underworld development team with AI support from Claude Code