Skip to content

Add dmfield_evaluate from PETSc - #436

Open
bknight1 wants to merge 7 commits into
developmentfrom
feat/dmfield_evaluate_function
Open

Add dmfield_evaluate from PETSc#436
bknight1 wants to merge 7 commits into
developmentfrom
feat/dmfield_evaluate_function

Conversation

@bknight1

Copy link
Copy Markdown
Member

DMFieldEvaluate computes the FE interpolant at each query point. For fields that are in the FE space (e.g., a quadratic velocity field on P2 elements), the gradient is exact to machine precision at any interior point.

https://petsc.org/main/manualpages/DM/DMFieldEvaluate/

DMFieldEvaluate computes the FE interpolant at each query point.  For fields that are in the FE space (e.g., a quadratic velocity field on P2 elements), the gradient is exact to machine precision at any interior point.

https://petsc.org/main/manualpages/DM/DMFieldEvaluate/
Copilot AI review requested due to automatic review settings July 27, 2026 04:41
@bknight1
bknight1 requested a review from lmoresi as a code owner July 27, 2026 04:41

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

Adds a new uw.function.dmfield_evaluate() API backed by PETSc DMFieldEvaluate to compute FE-interpolated values and derivatives (gradient/Hessian) at arbitrary coordinates, aiming to avoid projection artifacts for derivative-derived quantities (e.g. strain rate, viscosity).

Changes:

  • Introduces a Python API (dmfield_evaluate) with a per-(mesh,var) cache and a cache-clear helper.
  • Adds a new Cython extension (_dmfield_wrapper) to create/evaluate/destroy PETSc DMField objects.
  • Adds tests for FE-exact value/gradient/Hessian evaluation and for per-cell volumes computed via PETSc FVM geometry.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 8 comments.

Show a summary per file
File Description
tests/test_0580_dmfield_evaluate.py New tests covering DMField-based evaluation and mesh cell-volume totals.
src/underworld3/function/dmfield_evaluate.py New public API + cache management for DMField evaluation.
src/underworld3/function/_dmfield_wrapper.pyx New Cython/ctypes wrapper around PETSc DMFieldCreateDS / DMFieldEvaluate / DMFieldDestroy.
src/underworld3/function/init.py Exposes dmfield_evaluate* from underworld3.function.
src/underworld3/discretisation/discretisation_mesh.py Computes/stores _cell_volumes during nuke_coords_and_rebuild().
setup.py Registers the new _dmfield_wrapper extension for compilation.
Comments suppressed due to low confidence (3)

src/underworld3/function/_dmfield_wrapper.pyx:145

  • This call site also hard-codes libpetsc.dylib and assumes PETSC_DIR is set (os.environ["PETSC_DIR"]), which will raise KeyError in environments where petsc4py is available but PETSC_DIR/ARCH is not exported. Use the same cross-platform PETSc library resolution as in create() (and ideally cache the loaded CDLL).
        # Call via ctypes
        lib_path = os.path.join(
            os.environ["PETSC_DIR"],
            os.environ.get("PETSC_ARCH", ""),
            "lib", "libpetsc.dylib"
        )
        lib = ctypes.CDLL(lib_path)

src/underworld3/function/_dmfield_wrapper.pyx:180

  • Same portability issue in destroy(): libpetsc.dylib is macOS-specific and os.environ["PETSC_DIR"] can raise KeyError. Reuse the cross-platform PETSc library resolution logic (or a shared helper) here too.
            lib_path = os.path.join(
                os.environ["PETSC_DIR"],
                os.environ.get("PETSC_ARCH", ""),
                "lib", "libpetsc.dylib"
            )
            lib = ctypes.CDLL(lib_path)

tests/test_0580_dmfield_evaluate.py:82

  • Same MPI-safety issue here: v.sum() is rank-local. Use an allreduce before comparing to the global area (2.0).
        mesh = uw.meshing.UnstructuredSimplexBox(
            minCoords=(0, 0), maxCoords=(2, 1), cellSize=0.3
        )
        v = mesh._cell_volumes
        assert (v > 0).all()
        assert np.allclose(v.sum(), 2.0, rtol=1e-3)


💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment on lines +156 to +157
# Keep weakrefs so the cache is cleaned up when mesh/var die
_cache_owners[key] = (weakref.ref(mesh), weakref.ref(var))
Comment on lines +95 to +105
**Data freshness** — The variable's internal PETSc local vector
(``_lvec``) is read at the time of the *first* call. If the variable's
data changes later (e.g. after a Stokes solve), call
``mesh.update_lvec()`` before ``dmfield_evaluate()`` to sync the local
vector.

**Parallel** — This function is not yet collective across MPI ranks.
Each rank evaluates its own set of coordinates. For parallel-safe
evaluation across the whole domain, use ``uw.function.global_evaluate``
with an expression that has been pre-projected.

Comment on lines +56 to +66
# Locate libpetsc
petsc_dir = os.environ.get("PETSC_DIR", "")
petsc_arch = os.environ.get("PETSC_ARCH", "")
lib_dir = os.path.join(petsc_dir, petsc_arch, "lib")
if not os.path.isdir(lib_dir):
lib_dir = os.path.join(petsc_dir, "lib")
lib_path = os.path.join(lib_dir, "libpetsc.dylib")
if not os.path.exists(lib_path):
raise RuntimeError(f"Cannot find libpetsc at {lib_path}")

lib = ctypes.CDLL(lib_path)
Comment on lines +99 to +103
gradient : bool
If True, return first derivatives D.
hessian : bool
If False, return second derivatives H.

Comment on lines +149 to +156
ierr = lib.DMFieldEvaluate(
c_void_p(self._field_handle),
pts_py.handle,
0, # PETSC_REAL
B.ctypes.data_as(c_void_p) if B is not None and B.size else ctypes.c_void_p(),
D.ctypes.data_as(c_void_p) if D is not None and D.size else ctypes.c_void_p(),
H.ctypes.data_as(c_void_p) if H is not None and H.size else ctypes.c_void_p(),
)
Comment thread tests/test_0580_dmfield_evaluate.py Outdated

@pytest.fixture(scope="module")
def tri_mesh():
"""Unstructured simplex box — non-affine elements."""
Comment thread tests/test_0580_dmfield_evaluate.py Outdated
Comment on lines +68 to +72
mesh = uw.meshing.StructuredQuadBox(elementRes=(8, 8))
v = mesh._cell_volumes
assert (v > 0).all()
assert np.allclose(v.sum(), 1.0, rtol=1e-3)

Comment on lines +2709 to +2713
# Cell volumes via PETSc's FVM geometry (one call per cell).
cStart, cEnd = self.dm.getHeightStratum(0)
self._cell_volumes = numpy.abs(numpy.array(
[self.dm.computeCellGeometryFVM(c)[0]
for c in range(cStart, cEnd)]))
@lmoresi

lmoresi commented Jul 27, 2026

Copy link
Copy Markdown
Member

Adversarial review — live probes on the PR head (5f2f907)

Reviewed both for inherent validity and for fit with this repo's evaluation stack (uw.function.evaluate / the Clement derivative path), per project review policy. All probes ran in an isolated worktree on macOS/arm64, PETSc 3.25.0 (real, double, 32-bit indices); parallel probes at np=2.

Verdict: 10 findings — 3 CRITICAL, 3 MAJOR, 4 MINOR. The primitive is valid in concept — petsc4py has no DMField binding (confirmed), and FE-exact one-sided gradients plus Hessians are a real capability UW3 lacks — but this implementation needs a substantial rework before it can merge, and both CI jobs on the PR are currently red.

Findings (ranked)

CRITICAL-1 — Silent wrong values for points raw DMLocatePoints cannot locate (the #390 bug class, un-mitigated) — PROBED

DMFieldEvaluate_DS calls raw DMLocatePoints (petsc dmfieldds.c:259) with none of the fallback ladder UW3 added in PR #395. Two failure modes observed, both bad:

  • Silent garbage: dropped points are never written, and the wrapper allocates outputs with np.empty — callers receive uninitialized heap contents with no error and no NaN. On StructuredQuadBox(8×8), P1 field s = x+2y: top-face interior points 7/7 silently wrong (max err 0.91; uw.function.evaluate on the same points: 4.4e-16). The full P2 node set — the nodal-fill use case — 22/289 (7.6%) silently wrong, including 14 strictly interior edge midpoints (returned values include stale previous results and raw 10000.0 garbage).
  • Unpredictable hard error: the same 289-point set evaluated as one batch raises PETSc error 62 ("Point could not be located") — whether you get the error or the silent corruption depends on batching.
  • Points at y = 1+1e-9: silently wrong on both quad and simplex meshes (evaluate degrades gracefully via its fallback).
  • Simplex meshes are clean (0/433 P2-node failures), but quad boxes are a supported mesh class and this PR's own tests use them.

Minimum mitigation: allocate outputs with np.full(..., np.nan) so unlocated points are detectable. Real fix: route location through the same ladder evaluate uses.

CRITICAL-2 — Documented D (and H) layout is transposed for vector fields; the docstring's own examples read the wrong quantity — PROBED

PETSc's layout is point-major, then component, then direction (rD[(i*Nc + j)*dimC + l], dmfieldds.c:568). The PR reshapes to (n_points, dim, nc) and documents D[k,i,j] = ∂(component j)/∂xᵢ. Probe with the asymmetric field v = (y, 0), where ∂vx/∂y = 1:

D[:,1,0] (docstring slot for dvx/dy): [0. 0. 0.]
D[:,0,1] (transposed slot):           [1. 1. 1.]

Hessian probe (w = (0, xy)) confirms the same transpose for H. Every shipped test passes because v=(x², y²) has a diagonal gradient and the strain-rate check sums the symmetric pair — insensitive to the transpose by construction; scalar fields coincidentally read correctly in both conventions. Anyone computing vorticity or any non-symmetric velocity-gradient quantity from the documented slots gets silently wrong physics. Correct reshape is (n_points, nc, dim) (or transpose and fix the docs/examples).

CRITICAL-3 — Parallel is completely broken, contradicting the docstring — PROBED (np=2)

Docstring: "Each rank evaluates its own set of coordinates." Reality: every call fails with PETSc error 56 ("Trying parallel point location: only local point location supported") — including each rank evaluating its own var.coords, and including a 0-point rank. Cause: the points Vec is created on the mesh communicator instead of COMM_SELF, so DMLocatePoints_Plex refuses. (At least it errors cleanly on all ranks rather than hanging.)

MAJOR-4 — CI is red on both jobs; the feature is macOS-only as written

  • test: all 11 dmfield tests fail on Linux — hardcoded libpetsc.dylib (Linux needs .so) and PETSC_DIR unset in CI degenerates the path. Also inconsistent internally: create() uses os.environ.get(...), evaluate()/destroy() use os.environ[...] (KeyError).
  • Deprecated-pattern scan: fails on dmfield_evaluate.py:176 [except-pass] (Charter §4: bare swallow with no sanctioned-failure comment).

MAJOR-5 — The cache leaks every mesh and variable it ever touches; the cleanup machinery is dead code — PROBED

CachedDMField holds strong refs to mesh and var; _cache_cleanup() is never registered as a finalizer; the _cache_owners weakrefs are never consulted. After del var; del mesh; gc.collect() the cache entry survives and the whole mesh (DM, lvec, DMField) is immortal. In adaptive/remeshing workflows this accumulates entire mesh hierarchies. (The strong refs incidentally moot the id()-reuse key collision — objects never die — which is not the defense intended.)

MAJOR-6 — No accuracy win over the existing Clement derivative path for smooth fields — PROBED

Head-to-head the PR never ran: P2 field sin(πx)sin(πy), ds/dx at 200 interior points vs analytic, three resolutions:

h dmfield L2 / max Clement path L2 / max
0.1 9.9e-3 / 4.3e-2 8.9e-3 / 2.4e-2
0.05 2.2e-3 / 7.8e-3 1.9e-3 / 8.4e-3
0.025 5.8e-4 / 2.3e-3 4.7e-4 / 1.3e-3

Recovery averaging is superconvergent for smooth fields; "FE-exact" is exact only w.r.t. the interpolant. DMField's genuine wins: machine precision for fields in the FE space (confirmed at 1e-16 — and a solved P2 velocity is in the FE space, so the strain-rate pitch is sound there), one-sided semantics at discontinuities, and Hessians (which UW3 currently has no direct source for).

MINOR-7 — The datatype argument 0 is PETSC_DATATYPE_UNKNOWN, not PETSC_REAL — wrong-but-lucky (petsc only branches on == PETSC_SCALAR, so 0 falls into the real branch by accident). Should be the named constant.

MINOR-8 — The Cython file buys nothing: it is ctypes at runtime, re-CDLL'd on every call, with a c_int-for-PetscInt hazard on 64-bit-index builds. The repo's established pattern for exactly this is a Cython extern — see DMPlexInsertBoundaryValues in src/underworld3/cython/petsc_extras.pxi (PR #412). The three signatures needed are trivial (petscdmfield.h:39,45,56); an opaque DMField typedef plus petsc4py's PetscDM/PetscVec eliminates the dylib lookup, the datatype literal, and the int-width assumption.

MINOR-9 — mesh._cell_volumes is an unrelated drive-by: an eager O(cells) Python loop on every mesh build and every deform(), with no consumer in the codebase (only the PR's own tests). Measured: ~1 s per million cells, paid per free-surface/MMPDE step forever. Should be a lazy cached property, in its own PR.

MINOR-10 — Docstring defects: hessian parameter description inverted; Notes tell the user to call mesh.update_lvec() first while the function already does it unconditionally; the module example teaches the transposed D indexing (CRITICAL-2); one dead branch (nc > 0 is always true).

Against the UW3 use case (candidate replacement for the Clement derivative sub-path)

  • Nodal fill: at gradient discontinuities dmfield returns the one-sided element gradient of whichever cell wins point location — arbitrary-sided, not a controlled choice (probe: T=|x−0.5| kink nodes, dmfield all −1, Clement averages to 0; real lid-driven P2 velocity, ∂vx/∂x at shared vertices: median difference 3.8e-2, max 20.0 at the lid corner). A semantic change, not a drop-in.
  • Boundary points: disqualifying as-is on quad/hex meshes (CRITICAL-1) — this primitive must adopt the feat(evaluate): measured point-location capability + NaN/RBF fallback ladder (#392) #395 location ladder or hard-fail before it can sit under uw.function.evaluate.
  • Where it genuinely earns its place: exact strain rates of solved FE velocity fields, Hessians for adaptivity metrics, one-sided evaluation where averaging is wrong.

Attacks that failed (confirmations)

Data freshness after var.data writes: fresh values (the unconditional update_lvec works). mesh.deform() with a cached DMField: exact values and gradients on the deformed mesh (DM/lvec mutated in place — the cache stays coherent). Values at true mesh nodes: exact to 4.4e-16 on both mesh classes; at arbitrary interior points dmfield's values (4.3e-4) beat evaluate's (2.7e-3) on the P2 sin field. The 13 shipped tests pass locally on macOS in ~3 s including the create/destroy loop. Zero-point parallel input errors cleanly rather than hanging. A batch location error does not poison the cache.

UNPROBED: 3D, np>2, periodic meshes, annulus/spherical, 64-bit-index and complex builds, hessian=True, gradient=False.

Recommended path

Worth landing eventually — UW3 has no Hessian source and no FE-exact gradient path, and this fills both — but as a rework: DMFieldCreateDS/Evaluate/Destroy externs in petsc_extras.pxi (the #412 pattern), COMM_SELF point Vecs, NaN-filled outputs with a hard check, the corrected (n, nc, dim) layout with fixed examples, a real cache lifecycle (finalizers, no strong refs), the location fallback (or an explicit simplex-only guard), and the _cell_volumes loop split out as a lazy property in its own PR.

Adversarial review per project policy — Underworld development team with AI support from Claude Code

@lmoresi lmoresi left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Requesting changes formally, per the adversarial review posted below (probe evidence and file/line specifics there). The capability is wanted — UW3 has no Hessian source and no FE-exact gradient path, and petsc4py has no DMField binding, so a wrapper is the right idea — but this implementation isn't at the bar we hold our own work to, and both CI jobs are red.

Blocking items before this can merge:

  1. No silent garbage (CRITICAL-1): outputs must be NaN-initialized with a hard failure (or the PR #395 location ladder) for unlocatable points. 7.6% of P2 nodes on a plain quad box currently come back as uninitialized memory with no error.
  2. Fix the D/H array layout and every example that teaches it (CRITICAL-2): PETSc is point-component-direction; the documented (n, dim, nc) indexing is transposed for vector fields. Add an asymmetric-field test (e.g. v=(y,0)) that actually pins the layout — the current symmetric tests cannot.
  3. Parallel must work or refuse loudly (CRITICAL-3): points Vec on COMM_SELF, plus an explicit statement of the rank-local contract.
  4. Replace the ctypes/dylib mechanism with Cython externs in petsc_extras.pxi — the established pattern from PR #412 (DMPlexInsertBoundaryValues). That removes the macOS-only library path, the per-call CDLL, the PETSC_DATATYPE_UNKNOWN literal, and the PetscInt-width hazard in one move, and it's what makes CI green on Linux.
  5. Real cache lifecycle (MAJOR-5): no strong mesh/var refs, registered finalizers — the current cache immortalizes every mesh it touches, which is disqualifying for adaptive/remeshing workflows.
  6. Split mesh._cell_volumes into its own PR as a lazy cached property — it's an eager per-deform() O(cells) Python loop with no consumer in this PR.
  7. Honest framing: the smooth-field head-to-head shows the Clement path is superconvergent and slightly better there; this primitive's real value is exact strain rates of solved FE fields, one-sided evaluation at discontinuities, and Hessians. Docs should say that, not 'avoids boundary artifacts' generally.

Happy to pair on the extern wrapper — the three signatures are small and the #412 pattern maps directly.

Underworld development team with AI support from Claude Code

Replace ctypes/DMField wrapper with Cython externs (petsc_extras.pxi
pattern from PR #412). Remove per-(mesh,var) cache in favour of
create-evaluate-destroy per call (~0.5 us overhead vs ~50-200 us
evaluate). Correct D/H array layout from (n, dim, nc) to (n, nc, dim)
matching PETSc point-component-direction storage. Initialise outputs
with NaN so unlocated points are detectable. Document the all-ranks
collective contract for DMFieldEvaluate. Remove _cell_volumes.

Addresses all 10 findings from the adversarial review:
- CRITICAL-1 (silent garbage): NaN init + test_unlocated_points_nan
- CRITICAL-2 (transposed D/H): correct (n, nc, dim) + asymmetric tests
- CRITICAL-3 (broken parallel): COMM_SELF Vec, all-ranks-must-call, MPI test
- MAJOR-4 (CI red / macOS-only): Cython externs replace ctypes/.dylib
- MAJOR-5 (cache leak): no cache — fresh DMField per call
- MAJOR-6 (accuracy docs): honest DMField vs Clement guidance
- MINOR-7 (datatype 0): PETSC_REAL named constant
- MINOR-8 (ctypes buys nothing): Cython externs earn their keep
- MINOR-9 (cell_volumes): removed, deferred to separate PR
- MINOR-10 (docstring defects): all fixed

Tests: 17 serial + 1 MPI (18 total) — passing on StructuredQuadBox,
UnstructuredSimplexBox, and mpiexec -n 2.
@bknight1

Copy link
Copy Markdown
Member Author

Response to adversarial review — all 10 findings fixed

This rework addresses every finding from the adversarial review. The core changes are summarised below, followed by a finding-by-finding response.

Design decisions

  • No cache: DMFieldCreateDS (~0.3 µs) + DMFieldDestroy (~0.25 µs) overhead is negligible vs DMFieldEvaluate (~50–200 µs). Create → evaluate → destroy per call. Eliminates MAJOR-5 (cache leaks) and 6 associated failure modes (weakref finalizers, use-after-free, immortal meshes, global mutable state, is_valid lifecycle, clear_cache API).
  • No RBF fallback: Unlocated points return NaN. DMField's value proposition is FE-exact evaluation — silently mixing RBF-interpolated values would undermine that contract.
  • No C-side helper: DMFieldEvaluate_DS uses continue on negative cell indices, so NaN initialisers survive. Verified by test_unlocated_points_nan.
  • COMM_SELF point Vecs: Each rank provides its own local coordinates. DMFieldEvaluate (and DMLocatePoints internally) is collective on the DM — all ranks must call, zero-point ranks pass an empty Vec.
  • Cython externs: Follow the petsc_extras.pxi pattern from PR CBF flux on driven boundaries, kdtree module identity, global velocity_error #412 (see DMPlexInsertBoundaryValues). No ctypes, no .dylib, no os.environ.

Finding-by-finding

Finding Status Fix
CRITICAL-1 — Silent garbage for unlocated points ✅ Fixed Outputs initialised with np.full(..., np.nan). PETSc's DMFieldEvaluate_DS skips unlocated cells — NaN survives. test_unlocated_points_nan verifies outside-domain points return NaN in all arrays. Upper-boundary face tested (test_upper_face_no_crash).
CRITICAL-2 — D/H layout transposed for vector fields ✅ Fixed Corrected reshape to (n_points, nc, dim) matching PETSc's point→component→direction C order. test_gradient_asymmetric_layout (v=(y,0)) pins D[k, j, i] = ∂ⱼ/∂xᵢ. test_hessian_asymmetric_layout (w=(x², xy)) pins H[k, j, i, l]. All existing gradient tests updated.
CRITICAL-3 — Parallel completely broken ✅ Fixed Points Vec created on COMM_SELF (rank-local data). DMFieldEvaluate is collective on the DM — all ranks call, zero-point ranks pass an empty Vec. test_parallel_rank_local verified with mpiexec -n 2 (rank 0: real points, rank 1: empty).
MAJOR-4 — CI red; macOS-only ✅ Fixed Hardcoded libpetsc.dylib replaced by Cython externs in petsc_extras.pxi linked at build time. No os.environ["PETSC_DIR"] calls. except: pass replaced with except Exception: + sanctioned-failure comment (passes deprecated-pattern scanner).
MAJOR-5 — Cache leaks every mesh ✅ Fixed No cache. Fresh DMFieldEvaluator created and destroyed per call. No _field_cache, no weakref finalizers, no dmfield_evaluate_clear_cache() API. test_repeated_create_destroy stress-tests 10 create→evaluate→destroy cycles with no PETSc leak or crash.
MAJOR-6 — No accuracy win over Clement path for smooth fields ✅ Fixed Docstrings now state honestly: DMField for FE-space fields (machine-precision), one-sided element-boundary evaluation, and Hessians. Clement path for smooth non-FE fields (superconvergent). Users directed to uw.function.evaluate for RBF-filled boundary values.
MINOR-7 — Datatype 0 = PETSC_DATATYPE_UNKNOWN ✅ Fixed PETSC_REAL named constant via anonymous cdef enum in petsc_extras.pxi.
MINOR-8 — Cython file buys nothing ✅ Fixed _dmfield_wrapper.pyx uses Cython externs — portable, type-safe, no per-call CDLL, no c_int-for-PetscInt hazard. The include "../cython/petsc_extras.pxi" path is the same pattern used by _function.pyx.
MINOR-9_cell_volumes unrelated drive-by ✅ Fixed Removed from discretisation_mesh.py. Deferred to a separate PR as a lazy @cached_property.
MINOR-10 — Docstring defects ✅ Fixed hessian parameter corrected to "If True". update_lvec note corrected (no manual sync needed — function calls it internally). Layout examples corrected to (n, nc, dim). Dead nc > 0 branch removed.

Test results

18 tests, 0 failures (17 serial + 1 MPI)

$ pytest tests/test_0580_dmfield_evaluate.py -v
17 passed, 1 skipped (MPI needs --with-mpi)
$ mpiexec -n 2 pytest ... --with-mpi
1 passed (test_parallel_rank_local)

Coverage includes: quadratic velocity fields (P2), scalar fields (P1), unstructured simplex and structured quad meshes, asymmetric layout guards for both gradient and Hessian, unlocated-point NaN semantics, empty-coordinates edge case, non-square component count (nc=3, dim=2), NULL-pointer PETSc paths (D=NULL, H=NULL), create/destroy stress test (10 cycles), and MPI collective contract (rank 0 with points + rank 1 with zero points).

Key API (corrected layout)

B, D, H = dmfield_evaluate(velocity, mesh.X.coords)
# B.shape == (n_points, nc)
# D.shape == (n_points, nc, dim)  — D[k, j, i] = ∂(component j)/∂xᵢ
# H.shape == (n_points, nc, dim, dim)

bknight1 added 3 commits July 28, 2026 09:33
…Wincompatible-pointer-types

PETSc's DMFieldCreateDS/Destroy expect DMField* (struct _p_DMField*),
not void*. macOS clang warns but Linux gcc errors with
-Wincompatible-pointer-types. Declare the opaque type via Cython's
ctypedef struct with name-mapping to "DMField".
…type

PETSc's DMField is typedef struct _p_DMField *DMField — a pointer
typedef. Cython's "DMField" name mapping generated 'struct DMField',
which does not exist in PETSc. This compiled with a warning on macOS
clang but errored on Linux gcc (-Wincompatible-pointer-types).

Mapping to 'struct _p_DMField' generates code matching PETSc's actual
struct tag, so Cython's _p_DMField* becomes struct _p_DMField * which
is exactly DMField.
…rrors on macOS

macOS clang treats -Wincompatible-pointer-types as a warning (not an
error) by default, while Linux gcc errors. Add these flags so pointer
type mismatches (like void* vs DMField*) fail the build on macOS too,
allowing CI-compatible testing before pushing.
@bknight1
bknight1 requested a review from lmoresi July 28, 2026 03:11
@lmoresi

lmoresi commented Jul 28, 2026

Copy link
Copy Markdown
Member

One structural requirement before merge: this must not become a second evaluation front door

The rework looks strong (verification round in progress), but there is one repositioning we need: dmfield_evaluate must not sit in the public uw.function namespace as a peer of evaluate.

The concern is concrete: uw.function.evaluate carries four things user code silently depends on — the units/non-dimensionalisation frame, the point-location fallback ladder (PR #395; raw DMLocatePoints misses are exactly what your NaN contract now surfaces rather than repairs), sympy expression composition, and the collective global_evaluate semantics. A visible uw.function.dmfield_evaluate WILL leak into user scripts as a faster "evaluate", and every one of those four contracts silently disappears. We don't want a second front door that diverges from the first.

Asks:

  1. Remove the re-export from uw.function.__init__ (line ~98). The module stays importable at its full path for internal/expert use — no public alias.
  2. Lead the module and function docstrings with the positioning: internal primitive — user code should call uw.function.evaluate, followed by the one-line list of what this bypasses (units frame, location fallback, expressions, collective evaluation) and the sanctioned uses (internal derivative/Hessian machinery, expert diagnostics, adaptivity metrics).
  3. Optional but preferred: underscore-prefix the module (_dmfield.py) so the contract is unmissable at the import site.

The intended end-state (from the round-1 review) is this primitive slotted UNDER uw.function.evaluate as the exact-derivative/Hessian provider — at that point users get its benefits through the front door that already handles units and location. Everything else about the rework can proceed to merge once the verification round confirms it.

Underworld development team with AI support from Claude Code

- Remove re-export from uw.function namespace (no public front door)
- Underscore-prefix module name (_dmfield_evaluate.py) per UW3 convention
- Lead docstrings with 'internal primitive — use uw.function.evaluate'
- List bypassed contracts (units, location fallback, expressions,
  collective semantics) and sanctioned uses

This prevents dmfield_evaluate from becoming a second evaluation front
door that silently diverges from uw.function.evaluate.
@bknight1

Copy link
Copy Markdown
Member Author

The structural requirement has been addressed in commit 95883b3.

Three changes:

  1. Re-export removed from uw.function.__init__. The module is importable
    only at its full path: from underworld3.function._dmfield_evaluate import dmfield_evaluate.
    import underworld3 as uw; uw.function.dmfield_evaluate now raises AttributeError.

  2. Module and function docstrings lead with "Internal primitive — user code should
    call uw.function.evaluate" and enumerate the four bypassed contracts (units frame,
    point-location fallback, sympy expression composition, collective global_evaluate
    semantics). Sanctioned uses are listed: exact derivatives of solved FE fields, Hessians
    for adaptivity, one-sided element-boundary probes, and internal expression-stack bypass.

  3. Module underscore-prefixed to _dmfield_evaluate.py, so the contract is
    unmissable at the import site.

The end-state integration (slotting this under uw.function.evaluate as the
derivative/Hessian engine) remains as follow-up work for a later PR.

@lmoresi

lmoresi commented Jul 28, 2026

Copy link
Copy Markdown
Member

Round-2 adversarial verification — the rework holds

Verified at head 6a94ba79 with live probes in an isolated worktree (macOS/arm64, PETSc 3.25.0), rerunning every round-1 attack. All three round-1 CRITICALs are demonstrably dead:

Round-1 CRITICAL Status Evidence
Silent np.empty garbage for unlocatable points DEAD The round-1 case (full P2 node set, quad box): batched → loud PETSc-62 RuntimeError; one-by-one → exactly 22 NaN, 0 garbage values across all probes. Far-outside and y=1+1e-9 → all-NaN. Never a stale number anywhere.
Transposed D/H layout DEAD v=(y,0): ∂vx/∂y lands in the documented D[k,0,1] = 1.0 exactly, off-slots ≤3e-15; Hessian w=(x²,xy) all slots correct ≤5.7e-14. Confirmed against dmfieldds.c's cD[(nc*p+g)*dimR+d] — the reshape is right, not lucky. Docs teach the correct indexing.
Parallel dead DEAD COMM_SELF point Vecs: np2, each rank its own partition coords — values 3.3e-16, gradients 1.4e-14, zero-point rank fine.

Also confirmed: no leak (1000 create-evaluate-destroy cycles: ΔRSS 0.02 MB, 135 µs/call — inside the claimed envelope); the _cell_volumes drive-by is gone (diff purely additive); Charter-clean (the one shutdown swallow carries its sanctioned-failure comment); ctypes/dylib fully replaced by the petsc_extras.pxi externs; Linux CI genuinely runs and passes the 17 serial tests (round-1's all-red is fixed).

Remaining conditions before merge

  1. (Required — docs only) The unlocated-point contract as documented is false: "outside points are returned as NaN" is not guaranteed. Actual behavior is NaN or PETSc error 62, batch- and partition-dependent (a point in a location-hash candidate cell that fails the interior test → error; a point absent from the SF graph → NaN). The code fails LOUD, which is safe — but the docstrings in dmfield_evaluate.py and _dmfield_wrapper.pyx must say "may raise for near-boundary/off-rank unlocatable points; pre-filter or catch", and the wrapper's NaN-mechanism comment ("continue on negative cell index") mis-describes the real path.
  2. (Required — one sentence) State the full parallel contract: all ranks must call (skipping → hard deadlock, probed: hung in update_lvec until MPI timeout), and every queried point must lie in the calling rank's local subdomain — a valid interior point owned by another rank gives NaN or error 62.
  3. (Required — from the maintainer comment above) The namespace repositioning: no public uw.function.dmfield_evaluate re-export; internal-primitive docstring positioning.
  4. (Non-blocking follow-up issue) Pre-locate points to restore uniform NaN semantics; note test_parallel_rank_local is skipped in CI (--with-mpi local-only), so CRITICAL-3's guard is not CI-exercised.

With 1–3 in, this merges. The rework is genuinely strong — the no-cache design eliminated an entire class of lifecycle failures, and the layout/NaN/parallel behavior now matches what the tests claim.

Adversarial review per project policy — Underworld development team with AI support from Claude Code

…eview

- _dmfield_wrapper.pyx: document dual failure mode (error 62 vs NaN)
- _dmfield_evaluate.py: add rank-local subdomain constraint to parallel note
- _dmfield_evaluate.py: document error 62/NaN outcomes + pre-filter guidance

Addresses remaining conditions from round-2 adversarial review (PR #436).
@bknight1
bknight1 dismissed lmoresi’s stale review July 28, 2026 14:01

Addressed per commit fb673db — docstring fixes for all three review conditions (unlocated-point docs, parallel contract, namespace)

@bknight1

Copy link
Copy Markdown
Member Author

All conditions verified and pushed to the PR head.

Condition verification:

  • ✅ Condition 1 (unlocated-point docs): wrapper and function docstrings
    now document the dual error-62/NaN outcome with pre-filter guidance
    (mesh.points_in_domain() or catch).
  • ✅ Condition 2 (parallel contract): rank-local subdomain constraint
    added to the parallel note — each point must lie in the calling rank's
    local subdomain.
  • ✅ Condition 3 (namespace): no public re-export from uw.function,
    module underscored (_dmfield_evaluate.py), all docstrings lead with
    "internal primitive — call uw.function.evaluate" and list bypassed
    contracts and sanctioned uses.

Ready to merge when you are.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants