Add dmfield_evaluate from PETSc

[bknight1]

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]

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.

[lmoresi]

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]

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

[bknight1]

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)