evaluate / global_evaluate: one pass over the locator, not five patches (anti-scaling to np=8, a smoothed substitution, an uncontained cell hint)

[lmoresi]

uw.function.evaluate / global_evaluate: the subsystem needs one pass, not five patches

We have five open issues against the point-evaluation entry points. They are not five
independent bugs. Reading the code and measuring at development @ e475246f, the
subsystem has one architectural problem and four defects that sit on top of it: point
location is done by a hand-rolled kd-tree + half-space walk that has no cheap rejection
path, is re-run several times per call, is not shared between the classifier and the
interpolator, and is bypassed entirely on the one path where its answer is trusted.

Everything below is measured on a 2-D unstructured simplex box, cellSize=1/100
(23 264 cells, 46 929 P2 DOF) unless stated, amr-dev, Darwin arm64, runs sequential
(concurrent PETSc contends ~26%).

---

1. Parallel anti-scaling: global_evaluate gets slower as ranks are added

Not previously filed. This is the finding with the largest practical cost.

Fixed global query set — 20 000 points, identical on every rank, so the answer is
rank-count independent and wall time should be flat or better as np grows. Best of 3
calls, max over ranks:

np	global_evaluate (s)	vs np=1	efficiency T₁/(np·T_np)
1	0.0790	1.00×	1.00
2	0.1325	0.60×	0.30
4	0.1789	0.44×	0.11
8	0.2232	0.35×	0.044

np=8 is 2.83× slower than serial for identical work.

evaluate on the same fixed global set (timing only — in parallel it classifies against
the rank-local mesh, so the answer is rank-local): 0.0518 / 0.0436 / 0.0382 / 0.0359 s at
np = 1/2/4/8. Efficiency 1.00 / 0.59 / 0.34 / 0.18. It is nearly flat because
points_in_domain sees all 20 000 points on every rank at every np (measured: 60 000
points classified per 3 calls, unchanged from np=1 to np=8) — per-rank work proportional
to the GLOBAL point count, the classic signature.

The benign case, for contrast: rank-local nodal points (each rank asks only for its own
DOF coordinates) scale acceptably. At cellSize=1/150 (105 185 P2 DOF), evaluate repeat
cost 0.2623 / 0.1568 / 0.0730 s at np = 1/2/4 → efficiency 1.00 / 0.84 / 0.90. Answers are
correct: max|evaluate − nodal data| = 5.55e-16 at every np.

Where the time goes

Phase attribution inside global_evaluate, 3 calls, 20 000-point global set:

phase	np=1	np=4	np=8
Swarm.migrate (forward)	0.0743 s (31%)	0.3486 s (64%)	0.4656 s (69%)
…of which mesh.points_in_domain	0.1104 s	0.3551 s	0.4359 s
points_in_domain calls per global_evaluate	2	5	6
Swarm + 4 SwarmVariable construction	0.0021 s (0.9%)	0.0023 s (0.4%)	0.0023 s (0.3%)

Constructing a throwaway swarm and four swarm variables per call (_function.pyx:447-492)
is not the problem. The migration claim loop is: it calls the full point classification once
per round, and the number of rounds grows with rank count (2 → 6).

The amplifier: a point you do not own costs 51× a point you do

_get_closest_local_cells_internal (discretisation/discretisation_mesh.py:6073-6101):
points the first containment test rejects enter a k = min(n_cells, 50) nearest-centroid
walk whose only exit is every lost point found.

query set (4 000 points)	µs/point	in-cell test calls	points tested / points asked
owned by this rank (its own nodal coords)	0.67 – 0.77	2	1.0
not owned by this rank	7.45 – 8.91	51	50.7 – 51.0

10–13× wall clock, exactly 51× the containment work. Two compounding faults in that loop:

• no cheap rejection. A point beyond max_radius of the nearest local centroid cannot
  be in any local cell, but nothing tests that — a genuinely foreign point runs all 50
  iterations.
• the working set never shrinks. coords[lost_points] is re-tested in full every
  iteration, including points already located; and one unfindable point forces all 50
  iterations for the whole batch.

In parallel, the fraction of points a rank does not own is what grows with np. That is
the whole mechanism.

Secondary: a hidden O(cells) rebuild inside the first evaluate after every mesh move

_mark_faces_inside_and_out (discretisation_mesh.py:5505-…) is a per-cell Python loop
that builds the face control points, triggered lazily by the first containment test:

local cells	first call	cached
23 264	0.709 s	2 µs
5 726	0.177 s	2 µs

30.5 / 30.9 µs per local cell, linear. It is invalidated by nuke_coords_and_rebuild
(discretisation_mesh.py:2781-2782), i.e. by every deform() and every adapt() — so a
moving-mesh or adaptive run pays it once per mesh update, buried inside whatever evaluate
happens to be first.

This is also the whole of the apparent "nodal points are 6.7× slower than scattered points"
effect: cold, at cellSize=1/150, nodal 2.0347 s vs scattered-interior 0.3037 s — but
1.596 s of the 1.73 s gap is this one-time build, charged to whichever point set ran first.
Warm, the nodal/scattered gap is 5% (0.2607 vs 0.2476 s). Nodal points do take ~2× the
in-cell tests (they sit on cell faces, so more of them miss the first test), but at these
sizes that is ~4% of the call.

Is it a regression? No — we tested and rejected the three plausible recent suspects

1. Not the DMInterpolation cache. Hit/miss was 1 miss + 2 hits in every configuration;
   the coords hash is stable across repeat calls, and cache-hit cost is flat with np.
2. Not the out-of-domain allgather block (_function.pyx:588-640, added
   f6f0b963, 2026-06-06). local_fallback=False changes nothing: np=4 0.1789 → 0.1802 s,
   np=8 0.2232 → 0.2278 s. The extrapolated-flag count is 0 for in-domain nodal queries, so
   the block does not fire.
3. Not c52201b8 ("parallel locator hardening", 2026-05-28), which introduced the
   parallel-only _robust_owning_cells / _eval_use_robust_location. Forcing
   _eval_use_robust_location() False at run time — the pre-c52201b8 routing — moves
   nothing: np=8 global_evaluate 0.2249 → 0.2217 s, evaluate/local-nodal 0.0177 →
   0.0177 s. Both routes go through the same 50-neighbour walk, which predates it (6cb735b7).

The anti-scaling is structural, not a recent regression. It has simply never been
measured, because the rank-local pattern (which does scale) is what the tests exercise.

---

2. #491 — evaluate silently substitutes a smoothed degree-1 projection for the expression

Reproduced verbatim at head:

evaluate(d)            min -5.6261  max +6.0884
evaluate(d)**2         min +0.0000   negatives 0            <- correct
evaluate(d**2)         min -2.6386   negatives 614 of 4000
evaluate(sqrt(d**2))   min -0.0622   negatives 25
max|evaluate(d**2) - evaluate(d)**2| = 3.3521
CONTROL (no derivative): max|evaluate(f**2) - evaluate(f)**2| = 8.7e-19

Mechanism (located by the parallel triage session, confirmed here): when any derivative
appears anywhere in the expression, evaluate_nd (_function.pyx:964-985) replaces the
whole composite with an L2 Projection onto a hard-coded degree-1 continuous work
variable plus a gradient penalty (_function.pyx:748-767). The user's expression is never
evaluated pointwise; a smoothed continuous interpolant of it is.

Verdict: BUG, not a contract issue. The decisive number is not the square — it is the
exponential:

expression	evaluate min
d**2	−2.63858
d*d	−2.63858
sqrt(d**2)	−0.06218
f*d (mixed)	−3.22276
exp(d)	−7.25198

exp is strictly positive on the reals. No interpolation-undershoot argument survives that.
"Compose invariants in numpy" is a workaround for a substitution the caller was never told
about, not a documented contract. Unsafe shapes: every expression whose tree contains a
Derivative and is not itself linear in that derivative — squares, products of
derivatives, sqrt, exp, any yield/invariant law. Linear-in-derivative expressions
survive only because projection commutes with them.

---

3. #432 — P1 fields wrong at points exactly on 3-D cell edges

Mechanism (located by the parallel triage session): serial simplex meshes take an
unconditionally "exact" location capability (discretisation_mesh.py:6346-6347), which
makes the cell hint authoritative, so petsc_interpolate uses get_closest_cells
(:5941-5985) — a kd-tree nearest-control-point lookup with no containment test at
all — and DMLocatePoints is bypassed (petsc_tools.c:72-99). The only remaining guard
is a componentwise ξ∈[−1,1] box clamp (petsc_tools.c:285-299), which is a no-op on a
tetrahedron (the reference tet is not the reference box). A point on a shared edge gets
whichever cell owns the nearest control point; if that is not a cell containing the point,
the basis is evaluated at reference coordinates outside the tet and extrapolates.

Our instrumentation corroborates the bypass independently: on the serial simplex cache-miss
path we measured get_closest_cells called with all 20 000 query points and zero
containment tests performed on them.

This is the same class as #390 (quad TOP-face silent zeros) — a UW3 hint-policy problem, not
a PETSc one — and it is a recurrence: it has been fixed and returned. There is no regression
test pinning on-edge/on-face queries.

---

4. #314 / #513 — the np≥4 hangs are one defect, two triggers (#405 is separate)

Mechanism (located by the parallel triage session; #513 now closed into #314): the number
of DMLocatePoints calls on the world mesh DM inside petsc_interpolate is a rank-local
function of that rank's coordinate array. Two rank-local switches produce it:

• the DMInterpolationCache is keyed on a per-rank coords hash (_function.pyx:1293,
  dminterpolation_cache.py:114-135); a hit skips create_structure → DMLocatePoints
  entirely, and hit-vs-miss is decided per rank;
• _location_capability() is deliberately not reduced (discretisation_mesh.py:6334-6339)
  yet is part of that cache key, so a capability split forces a hit on some ranks and a miss
  on others.

The entry guards (_function.pyx:1166-1170, :1012-1016) correctly get every rank into
petsc_interpolate; nothing keeps them in lockstep once inside. This resolves #513's own
caveat — its non-empty control also hung because emptiness was never the mechanism;
per-rank divergence of coords_array[in_or_not] (_function.pyx:1027) is, and rank-distinct
partially-off-domain points produce that just as reliably as a zero-interior-point rank.

#405 is a different defect and fires earlier: a rank-asymmetric ValueError from
unguarded _radii.min()/.max() (discretisation_mesh.py:6512, 6527) and four sibling sites,
before this path is reached. get_min_radius feeds estimate_dt and penalty scaling at 14
call sites and sits under Swarm.migrate, so a zero-cell rank takes down ordinary
time-stepping, not just evaluation. Fix order: #405, then #314.

---

5. #279 — the ND↔units boundary is still unenforced

Confirmed still reproducing at head by the parallel triage session: dimensional
UnitAwareArray coordinates (250 km, 750 km) against a unit-box ND DM return 0.800 and
0.933 instead of 0.25 and 0.75 — no exception, no warning. evaluate_nd strips the
subclass with np.array(...).view(np.ndarray) (_function.pyx:949) and keeps the numeric
value. Small, mechanical, but it is the exact mechanism by which #267 slipped in.

---

What a fix has to address

One root cause, in the locator. Items 1, 3 and (via the cache key) 4 are all
consequences of point location being a bespoke kd-tree + half-space walk with no clean
contract. The work is:

1. Give the walk a rejection path and a shrinking working set
   (_get_closest_local_cells_internal:6073-6101). A bounding-box / max-radius pre-test
   turns the 51× foreign-point penalty into ~1×; dropping located points from the working
   set each iteration removes the rest. This alone should recover most of §1.
2. Locate once per call. points_in_domain locates the near-boundary points and throws
   the owning cells away; petsc_interpolate then re-locates every interior point
   (_function.pyx:1026-1032). One classification should produce both the in/out mask and
   the cell hints.
3. Make the cell hint containment-checked, or stop calling it authoritative (§3). Simplex
   meshes assert "exact" unconditionally (:6346-6347) while the hint they hand over is a
   nearest-control-point lookup. Either check containment before trusting it, or restrict
   the bypass to meshes where nearest-control-point provably equals owning-cell. A ξ-clamp
   in barycentric coordinates would also close it for simplices; the current box clamp
   cannot.
4. Make the collective structure rank-independent (§4). The DMLocatePoints call count
   must not depend on rank-local coords. Either reduce the cache decision and the location
   policy across ranks, or move the collective out of the cached region entirely.
5. Evaluate expressions pointwise (§2). Project the derivative leaves to nodal values,
   then compose the expression pointwise at the query points — do not project the composite.
   Failing that, evaluate must refuse nonlinear-in-derivative expressions with an
   actionable error rather than return a smoothed surrogate. The degree-1 hard-coding is
   independently wrong for a P2 field.
6. Cheap, separable follow-ups: guard the empty-rank reductions (Empty-rank support, layer 2: evaluate/points_in_domain, radii/centroid reductions, gather_data NaN-stripping #405, 5 named sites);
   reject or non-dimensionalise UnitAwareArray coordinates at the ND entry points (Enforce the ND<->units boundary: evaluate/global_evaluate should reject or nondimensionalise UnitAwareArray input #279);
   rebuild _mark_faces_inside_and_out incrementally or vectorise it (30 µs/cell, redone on
   every deform).

Genuinely separate work: #405 (rank-asymmetric raises — different failure mode, fires
first, blocks the rest) and #279 (an API guard, not a locator problem). Everything else is
one pass over the locator plus the derivative-path fix.

Regression tests the fix must land with, since #432/#390 have both recurred:

• on-edge / on-face / shared-vertex queries in 2-D and 3-D, against the closed-form P1
  edge value (1−t)·u_a + t·u_b (tests/test_0753_nested_mg_prolongation.py::test_reproduces_an_arbitrary_coarse_field
  already computes exactly this and deliberately avoids evaluate);
• evaluate(exp(d)) > 0 and evaluate(d**2) ≥ 0 for a P2 field;
• a np=1/2/4/8 wall-time ceiling on a fixed global query set — the anti-scaling in §1 is
  invisible to every existing test because they all use the rank-local pattern.

---

Issues superseded

issue	why it folds in
#491	The mechanism (whole-composite → degree-1 continuous projection) is item 5 above; it is the derivative half of the same "evaluate does not do what it says" problem.
#432	The mechanism (uncontained nearest-control-point hint + a box clamp that is a no-op on a tet) is item 3; it is the same locator. Recurrence of #390 — closing it again without a pinned regression test just restarts the clock.
#314	The mechanism (rank-local DMLocatePoints call count via the coords-hash cache key and the unreduced location policy) is item 4; it is the same locator, seen through the cache. #513 already closed into this.
#405	Keep open and fix first. Listed here for ordering only — it is a rank-asymmetric raise, not a locator defect, and it fires before any of the above.
#279	Fold in as item 6 — it is an entry-point guard on the same two functions and is a 10-line change; keeping it separate costs more than carrying it.

Underworld development team with AI support from Claude Code

[lmoresi]

Items 1-3 are in PR #556. One correction to this issue's own analysis, from the measurement:

Item 1 does NOT recover the anti-scaling, and §1's attribution was wrong about which walk pays it. points_in_domain already had a cheap rejection of its own (the boundary-sign / far_from_domain test), so the migrate claim loop was never paying the 51x. After the fix, migrate is still 68.8% of np=8 and global_evaluate on the fixed global set is essentially unchanged (0.2564 -> 0.2319 s at np=8). The residual is the claim loop classifying the whole global point set on every rank every round — that is item 4, not item 1.

So the anti-scaling headline moves to item 4 in full, and item 4 is now the highest-value remaining work in this issue.

What items 1-3 did deliver (all measured, in the PR):

• the 51x foreign-point penalty is gone where it actually lived, in _robust_owning_cells: foreign points 7.99 -> 0.53 us/pt at np=1 and 9.55 -> 0.77 us/pt at np=4; containment tests per foreign point 51.0 -> 1.0. Foreign/owned cost ratio 13.1-13.8x -> 1.1-1.2x. evaluate point-tests per located point at np=8: 4.19 -> 1.27.
• item 2 landed — one classification now produces both the in/out mask and the cell hints, with a lazy fill (an eager fill made serial evaluate 29% slower by searching on cache hits).
• item 3 / uw.function.evaluate returns wrong values for P1 fields at points exactly on cell edges (3D) #432 fixed by containment-checking the hint (option (a)); three location branches collapse to one. Fail-before validated: 3-D P1 evaluation at t=0.25 along cell edges was off by 3.04e-01 at 17 of 1318 points.
• Answer identity: 15 of ~450 000 located cells differ, all 3-D at the loose tolerance, and every one verified to be a point both cells contain — the winner used to be the last of up to 50 rounds and depended on whether an unrelated point in the same batch was findable; it is now the nearest containing centroid.

Two findings for the follow-on work, noticed and deliberately not fixed:

• points_in_domain calls get_max_radius(), which is @uw.collective_operation and allgathers. It does not fire today, but it is a collective inside the routine evaluate uses to decide which ranks hold points — exactly the item-4 shape, and its lines sit in the region the Empty-rank support, layer 2: evaluate/points_in_domain, radii/centroid reductions, gather_data NaN-stripping #405 fix is touching.
• Mesh._build_kd_tree_index_PIC has no callers anywhere in the tree.

Underworld development team with AI support from Claude Code

[lmoresi]

Items 1-3 landed in PR #556 (merged 3f54218e). Also closed on the way: #558 — the dead RBF rung was the cause of the parallel NaNs too, 0/18/120 at np1/2/4 to 0/0/0.

Remaining work, and a correction to the numbering. What this issue calls "item 4" is two separable defects, and they should not be confused because fixing the first does nothing for the second:

4a — the hang (correctness). The DMLocatePoints call count inside petsc_interpolate is a rank-local function of the rank's coordinate array: the DMInterpolationCache is keyed on a per-rank coords hash (_function.pyx:1293, dminterpolation_cache.py:114-135) and _location_capability() is deliberately unreduced (discretisation_mesh.py:6334-6339) yet is part of that key. Two triggers reach it — a zero-interior-point rank, and rank-distinct partially-off-domain points.

4b — the anti-scaling (performance). Swarm.migrate's claim loop classifies the whole global point set on every rank, every round, and the rounds grow 2 -> 5 -> 6 with rank count. That is 69% of np=8 time and it is the reason global_evaluate is 2.83x slower at np=8 than serial. This is the item with the practical cost, and #556 confirmed it is untouched by the locator work.

5 — pointwise composition (#491's mechanism): project the derivative leaves to nodal values and compose the expression pointwise, rather than L2-projecting the whole composite onto a hard-coded degree-1 continuous work variable. The degree-1 hard-coding is independently wrong for a P2 field.

6 — the cheap ones: the 50-neighbour cap still losing genuinely interior points (6 of 3000 on a graded 3-D mesh — they now take an RBF value rather than a NaN, but a fallback is not a location); _mark_faces_inside_and_out at 30 us/cell, rebuilt on every deform()/adapt(); the ND<->units guard (formerly #279).

Note for whoever takes 4a: points_in_domain calls get_max_radius(), which is @uw.collective_operation and allgathers — a collective inside the routine evaluate uses to decide which ranks hold points. It does not fire today but it is the same shape as the defect.

Underworld development team with AI support from Claude Code