diff --git a/docs/developer/guides/memory-diagnostics.md b/docs/developer/guides/memory-diagnostics.md
new file mode 100644
index 000000000..8dea38e70
--- /dev/null
+++ b/docs/developer/guides/memory-diagnostics.md
@@ -0,0 +1,153 @@
+# Memory diagnostics
+
+Long parallel runs occasionally OOM on HPC even when each step looks small.
+The `uw.utilities.memprobe` module gives you a way to "light up" memory
+tracking on demand, sample at regular intervals, and pin which subsystem is
+growing.
+
+## Quick start
+
+```bash
+UW_MEMPROBE=1 mpirun -n 16 python my_long_run.py
+```
+
+With this flag set, the `Stokes.solve()` and `NavierStokes.solve()` paths
+emit a one-line growth report each time they're called:
+
+```
+[memprobe] Stokes.solve:
+  RSS +0.42 MiB
+  kdtree: live +1, total_constructed +1
+```
+
+A clean run shows mostly zero deltas. A leak shows the same component
+growing on every step.
+
+## What's tracked
+
+| Signal | Source | Cost |
+|---|---|---|
+| Process RSS (MiB, current) | `psutil.Process().memory_info().rss`, fallback `/proc/self/statm` (Linux), last resort `resource.ru_maxrss` | free |
+| KDTree live count | `uw.kdtree.live_count()` | free |
+| KDTree total constructed | `uw.kdtree.total_constructed()` | free |
+| Per-class Python instance counts | `gc.get_objects()` walk | slow — gated behind `full=True` |
+
+The RSS source is **current** RSS where possible — it should drop when memory
+is freed, not just rise. The last-resort `resource.ru_maxrss` fallback returns
+the **peak** (high-water-mark) RSS instead and never decreases, so on systems
+where neither psutil nor `/proc/self/statm` is available you'll see growth but
+not recovery; install `psutil` to fix that.
+
+`KDTree` instances are tracked via Cython class counters in `__cinit__` and
+`__dealloc__`. CPython refcounting calls `__dealloc__` promptly when the
+refcount hits zero, so the count is accurate for typical use; it can lag if a
+KDTree ends up in a reference cycle that only the cyclic garbage collector
+can break. Call `gc.collect()` before reading if that matters.
+
+PETSc-side object and allocation tracking is **not** parsed from Python —
+PETSc's own `-log_view` and `-malloc_dump` runtime flags give the same
+information more reliably. To enable them from Python:
+
+```python
+from underworld3.utilities import memprobe
+memprobe.dump_petsc_leaks_at_finalize()  # equivalent to -malloc_dump -objects_dump
+```
+
+## API
+
+### Snapshots and diffs
+
+```python
+from underworld3.utilities import memprobe
+
+before = memprobe.snapshot()
+do_work()
+after = memprobe.snapshot()
+print(memprobe.format_diff("after-work", memprobe.diff(before, after)))
+```
+
+Add `full=True` to also walk Python-class counts:
+
+```python
+snap = memprobe.snapshot(full=True)
+# snap["py_classes"] = {"underworld3.swarm.Swarm": 3, ...}
+```
+
+### Probe context manager
+
+```python
+with memprobe.probe("step 42"):
+    advance_one_step()
+# On exit, the diff is emitted via `print` (configurable).
+```
+
+The `emit` keyword takes any callable: `with probe(..., emit=logger.info):`
+to route into the logging system, or a rank-aware writer for parallel runs:
+
+```python
+import underworld3 as uw
+emit = lambda s: uw.pprint(0, s)  # rank-0 only
+with memprobe.probe("step 42", emit=emit):
+    ...
+```
+
+### Decorator
+
+```python
+@memprobe.instrument("my-hot-loop")
+def step():
+    ...
+```
+
+When `memprobe.ENABLED` is `False` (the default) the decorator's wrapper is
+a single attribute lookup + branch — sub-microsecond — so it's safe to
+leave on hot paths permanently.
+
+`Stokes.solve()` and `NavierStokes.solve()` are pre-decorated. Add more if
+you want them.
+
+### Runtime toggles
+
+```python
+memprobe.enable()
+# ...probed region...
+memprobe.disable()
+```
+
+`UW_MEMPROBE=1` flips it on at import time.
+
+## Debugging recipes
+
+### "RSS grows X MiB per step — which component is it?"
+
+1. Set `UW_MEMPROBE=1` and run for ~20 steps to confirm the per-solve
+   growth pattern.
+2. Add `with memprobe.probe("step N", full=True):` around your step loop.
+   The `full=True` walks `gc.get_objects()` and lists Python class growth
+   sorted by absolute change — usually the dominant suspect is on top.
+3. If RSS grows but no Python class does, the leak is in PETSc memory or
+   C extensions. Re-run with `-log_view -malloc_dump` and inspect the PETSc
+   reports written at finalize.
+
+### "Are kd-trees being released properly?"
+
+Check `uw.kdtree.live_count()` directly, or look for the `kdtree: live +N`
+line in the diff. KDTrees should drop to zero when their owning object
+(typically a `Mesh` or `Swarm`) is destroyed.
+
+### Parallel runs
+
+`memprobe` runs on each rank independently. For meaningful aggregate
+output, pipe `emit` through a rank filter:
+
+```python
+import underworld3 as uw
+def root_only(s):
+    if uw.mpi.rank == 0:
+        print(s)
+
+with memprobe.probe("step", emit=root_only):
+    ...
+```
+
+Or compare per-rank snapshots manually for a load-imbalance view.
diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx
index 99a13087f..19acaedbb 100644
--- a/src/underworld3/ckdtree.pyx
+++ b/src/underworld3/ckdtree.pyx
@@ -15,6 +15,25 @@ cdef extern from "kdtree_interface.hpp" nogil:
         void find_closest_point( size_t  num_coords, const double* coords, long unsigned int* indices, double* out_dist_sqr, bool* found )
         size_t knnSearch(const double* query_point, const size_t num_closest, long unsigned int* indices, double* out_dist_sqr )
 
+# Module-level live-instance counter for memory introspection.
+# Incremented in __cinit__, decremented in __dealloc__. CPython refcounting
+# calls __dealloc__ promptly when the refcount hits zero, so the count is
+# accurate for typical use; it can lag if a KDTree ends up in a reference
+# cycle that only the cyclic garbage collector can break — call
+# gc.collect() before reading if that matters. Read via
+# uw.utilities.memprobe.snapshot() or directly via uw.kdtree.live_count().
+cdef long _live_instances = 0
+cdef long _total_constructed = 0
+
+def live_count():
+    """Number of KDTree instances currently alive on this rank."""
+    return _live_instances
+
+def total_constructed():
+    """Total KDTree instances ever constructed on this rank."""
+    return _total_constructed
+
+
 cdef class KDTree:
     """
     Unit-aware KD-Tree for spatial indexing and queries.
@@ -84,10 +103,16 @@ cdef class KDTree:
         self.points = points
         self.index = new KDTree_Interface(<const double *> &points[0][0], points.shape[0], points.shape[1])
 
+        global _live_instances, _total_constructed
+        _live_instances += 1
+        _total_constructed += 1
+
         super().__init__()
 
     def __dealloc__(self):
         del self.index
+        global _live_instances
+        _live_instances -= 1
 
     @property
     def n(self):
diff --git a/src/underworld3/cython/petsc_types.pyx b/src/underworld3/cython/petsc_types.pyx
index 609dad2a6..247401832 100644
--- a/src/underworld3/cython/petsc_types.pyx
+++ b/src/underworld3/cython/petsc_types.pyx
@@ -1,9 +1,36 @@
-from libc.stdlib cimport malloc
+from libc.stdlib cimport malloc, free
 
 cdef class PtrContainer:
 
+    def __cinit__(self):
+        self.fns_residual = NULL
+        self.fns_bcs = NULL
+        self.fns_jacobian = NULL
+        self.fns_bd_residual = NULL
+        self.fns_bd_jacobian = NULL
+
+    def __dealloc__(self):
+        if self.fns_residual != NULL:
+            free(self.fns_residual)
+        if self.fns_bcs != NULL:
+            free(self.fns_bcs)
+        if self.fns_jacobian != NULL:
+            free(self.fns_jacobian)
+        if self.fns_bd_residual != NULL:
+            free(self.fns_bd_residual)
+        if self.fns_bd_jacobian != NULL:
+            free(self.fns_bd_jacobian)
+
     cpdef allocate(self, int n_res, int n_bcs, int n_jac, int n_bd_res, int n_bd_jac):
         """Allocate function pointer arrays of the given sizes."""
+
+        # Free existing memory if already allocated
+        if self.fns_residual != NULL: free(self.fns_residual)
+        if self.fns_bcs != NULL:      free(self.fns_bcs)
+        if self.fns_jacobian != NULL: free(self.fns_jacobian)
+        if self.fns_bd_residual != NULL: free(self.fns_bd_residual)
+        if self.fns_bd_jacobian != NULL: free(self.fns_bd_jacobian)
+
         self.fns_residual    = <PetscDSResidualFn*>   malloc(n_res    * sizeof(PetscDSResidualFn))
         self.fns_bcs         = <PetscDSResidualFn*>   malloc(n_bcs    * sizeof(PetscDSResidualFn))
         self.fns_jacobian    = <PetscDSJacobianFn*>   malloc(n_jac    * sizeof(PetscDSJacobianFn))
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index bf3e40752..35e321906 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -285,6 +285,11 @@ def __init__(
         self._registered_swarms = weakref.WeakSet()
         self._registered_surfaces = weakref.WeakSet()  # Surfaces using this mesh
         self._registered_submeshes = weakref.WeakSet()  # Submeshes from extract_region
+
+        # _mesh_update_lock: Re-entrant lock to coordinate mesh deformation.
+        # Held by mesh_update_callback during _deform_mesh(). Checked by
+        # MeshVariable callbacks (blocking=False) to skip PETSc sync during
+        # sensitive coordinate changes.
         self._mesh_update_lock = threading.RLock()
 
         comm = PETSc.COMM_WORLD
@@ -594,14 +599,25 @@ class replacement_boundaries(Enum):
         # to handle that so we just wrap it here.
 
         def mesh_update_callback(array, change_context):
-            print(f"Mesh update callback - mesh deform")
-            coords = array.reshape(-1, array.owner.cdim)
-            self._deform_mesh(coords, verbose=True)
+            mesh = array.owner
+            if mesh is None:
+                # This guard handles cases where the array is accessed during
+                # object teardown (e.g. at application exit or during mesh
+                # replacement), where the owning Python mesh object has already
+                # been garbage collected but the NDArray proxy still exists.
+                return
+
+            if verbose:
+                uw.pprint(0, f"Mesh update callback - mesh deform")
+
+            coords = array.reshape(-1, mesh.cdim)
+            mesh._deform_mesh(coords, verbose=verbose)
 
             # Increment mesh version to notify registered swarms of coordinate changes
-            with self._mesh_update_lock:
-                self._mesh_version += 1
-                print(f"Mesh version incremented to {self._mesh_version}")
+            with mesh._mesh_update_lock:
+                mesh._mesh_version += 1
+                if verbose:
+                    uw.pprint(0, f"Mesh version incremented to {mesh._mesh_version}")
 
             return
 
@@ -1337,10 +1353,14 @@ def _re_extract_from_parent(self, verbose=False):
             )
 
             def mesh_update_callback(array, change_context):
-                coords = array.reshape(-1, array.owner.cdim)
-                self._deform_mesh(coords, verbose=False)
-                with self._mesh_update_lock:
-                    self._mesh_version += 1
+                mesh = array.owner
+                if mesh is None:
+                    return
+
+                coords = array.reshape(-1, mesh.cdim)
+                mesh._deform_mesh(coords, verbose=False)
+                with mesh._mesh_update_lock:
+                    mesh._mesh_version += 1
                 return
 
             self._coords.add_callback(mesh_update_callback)
@@ -1777,42 +1797,43 @@ def _deform_mesh(self, new_coords: numpy.ndarray, verbose=False):
         The coord array that is passed in should match the shape of self.data
         """
 
-        coord_vec = self.dm.getCoordinatesLocal()
-        coords = coord_vec.array.reshape(-1, self.cdim)
-        coords[...] = new_coords[...]
-
-        self.dm.setCoordinatesLocal(coord_vec)
-        self.nuke_coords_and_rebuild()
-
-        # Rebuild the _coords array view.  nuke_coords_and_rebuild may
-        # replace the coordinate vector internally (createCoordinateSpace),
-        # leaving self._coords as a stale numpy view of the old buffer.
-        import underworld3.utilities
-        old_callbacks = getattr(self._coords, "_callbacks", [])
-        self._coords = underworld3.utilities.NDArray_With_Callback(
-            numpy.ndarray.view(
-                self.dm.getCoordinatesLocal().array.reshape(-1, self.cdim)
-            ),
-            owner=self,
-        )
-        for cb in old_callbacks:
-            self._coords.add_callback(cb)
-
-        # BUGFIX(#122): mark registered solvers for rebuild. Since PR #127
-        # ("Trust JIT cache: skip DM rebuild on constant-only parameter
-        # changes") a solver with is_setup=True trusts its cached PETSc DM
-        # / SNES assembly and skips rebuild on the next solve(). After a
-        # coordinate change the cached DM still carries pre-deform
-        # coordinates, so F(v_prev) ≈ 0 and the solver converges in 0
-        # iterations without updating the solution. mesh.adapt() already
-        # does this; _deform_mesh must match.
-        for solver in self._equation_systems_register:
-            if solver is not None and hasattr(solver, "is_setup"):
-                solver.is_setup = False
-
-        # Propagate coordinate changes to registered submeshes
-        for submesh in self._registered_submeshes:
-            submesh.sync_coordinates_from_parent()
+        with self._mesh_update_lock:
+            coord_vec = self.dm.getCoordinatesLocal()
+            coords = coord_vec.array.reshape(-1, self.cdim)
+            coords[...] = new_coords[...]
+
+            self.dm.setCoordinatesLocal(coord_vec)
+            self.nuke_coords_and_rebuild()
+
+            # Rebuild the _coords array view.  nuke_coords_and_rebuild may
+            # replace the coordinate vector internally (createCoordinateSpace),
+            # leaving self._coords as a stale numpy view of the old buffer.
+            import underworld3.utilities
+            old_callbacks = getattr(self._coords, "_callbacks", [])
+            self._coords = underworld3.utilities.NDArray_With_Callback(
+                numpy.ndarray.view(
+                    self.dm.getCoordinatesLocal().array.reshape(-1, self.cdim)
+                ),
+                owner=self,
+            )
+            for cb in old_callbacks:
+                self._coords.add_callback(cb)
+
+            # BUGFIX(#122): mark registered solvers for rebuild. Since PR #127
+            # ("Trust JIT cache: skip DM rebuild on constant-only parameter
+            # changes") a solver with is_setup=True trusts its cached PETSc DM
+            # / SNES assembly and skips rebuild on the next solve(). After a
+            # coordinate change the cached DM still carries pre-deform
+            # coordinates, so F(v_prev) ≈ 0 and the solver converges in 0
+            # iterations without updating the solution. mesh.adapt() already
+            # does this; _deform_mesh must match.
+            for solver in self._equation_systems_register:
+                if solver is not None and hasattr(solver, "is_setup"):
+                    solver.is_setup = False
+
+            # Propagate coordinate changes to registered submeshes
+            for submesh in self._registered_submeshes:
+                submesh.sync_coordinates_from_parent()
 
         return
 
@@ -3777,12 +3798,15 @@ def adapt(self, metric_field, verbose=False):
 
             # Rebuild the callback for mesh deformation
             def mesh_update_callback(array, change_context):
-                print(f"Mesh update callback - mesh deform")
+                if verbose:
+                    uw.pprint(0, f"Mesh update callback - mesh deform")
+
                 coords = array.reshape(-1, array.owner.cdim)
-                self._deform_mesh(coords, verbose=True)
+                self._deform_mesh(coords, verbose=verbose)
                 with self._mesh_update_lock:
                     self._mesh_version += 1
-                    print(f"Mesh version incremented to {self._mesh_version}")
+                    if verbose:
+                        uw.pprint(0, f"Mesh version incremented to {self._mesh_version}")
                 return
 
             self._coords.add_callback(mesh_update_callback)
diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py
index 510787459..9cd5943d3 100644
--- a/src/underworld3/discretisation/discretisation_mesh_variables.py
+++ b/src/underworld3/discretisation/discretisation_mesh_variables.py
@@ -492,36 +492,44 @@ def _create_variable_array(self, initial_data=None):
         # Single callback function (following mesh_update_callback pattern)
         def variable_update_callback(array, change_context):
             """Callback to sync variable changes back to PETSc (like mesh.points)"""
+            var = array.owner
+            if var is None:
+                # This guard handles cases where the array is accessed during
+                # object teardown (e.g. at application exit or mesh rebuilds),
+                # where the owning Python variable has already been garbage
+                # collected but the NDArray proxy still exists.
+                return
+
             # Only act on data-changing operations (following mesh.points pattern)
             data_changed = change_context.get("data_has_changed", True)
             if not data_changed:
                 return
 
             # Prevent recursion by checking if we're already in a callback
-            if hasattr(self, "_in_callback") and self._in_callback:
+            if hasattr(var, "_in_callback") and var._in_callback:
                 return
 
             # Set recursion guard
-            self._in_callback = True
+            var._in_callback = True
 
             try:
                 # Skip updates during mesh coordinate changes to prevent corruption
                 # Check if mesh is currently being updated
-                if hasattr(self.mesh, "_mesh_update_lock"):
+                if hasattr(var.mesh, "_mesh_update_lock"):
                     # Try to acquire lock without blocking - if we can't, skip update
-                    if not self.mesh._mesh_update_lock.acquire(blocking=False):
+                    if not var.mesh._mesh_update_lock.acquire(blocking=False):
                         return
                     try:
                         # Persist changes to PETSc (like mesh callback updates coordinates)
-                        self.pack_uw_data_to_petsc(array, sync=True)
+                        var.pack_uw_data_to_petsc(array, sync=True)
                     finally:
-                        self.mesh._mesh_update_lock.release()
+                        var.mesh._mesh_update_lock.release()
                 else:
                     # Fallback if no lock exists
-                    self.pack_uw_data_to_petsc(array, sync=True)
+                    var.pack_uw_data_to_petsc(array, sync=True)
             finally:
                 # Clear recursion guard
-                self._in_callback = False
+                var._in_callback = False
 
         # Register the callback (following mesh.points pattern)
         array_obj.add_callback(variable_update_callback)
@@ -556,34 +564,38 @@ def _create_flat_data_array(self, initial_data=None):
         # Callback for flat data format
         def flat_data_update_callback(array, change_context):
             """Callback to sync flat data changes back to PETSc"""
+            var = array.owner
+            if var is None:
+                return
+
             # Only act on data-changing operations
             data_changed = change_context.get("data_has_changed", True)
             if not data_changed:
                 return
 
             # Prevent recursion by checking if we're already in a callback
-            if hasattr(self, "_in_flat_callback") and self._in_flat_callback:
+            if hasattr(var, "_in_flat_callback") and var._in_flat_callback:
                 return
 
             # Set recursion guard
-            self._in_flat_callback = True
+            var._in_flat_callback = True
 
             try:
                 # Skip updates during mesh coordinate changes to prevent corruption
-                if hasattr(self.mesh, "_mesh_update_lock"):
-                    if not self.mesh._mesh_update_lock.acquire(blocking=False):
+                if hasattr(var.mesh, "_mesh_update_lock"):
+                    if not var.mesh._mesh_update_lock.acquire(blocking=False):
                         return
                     try:
                         # Use pack_raw for flat data format
-                        self.pack_raw_data_to_petsc(array, sync=True)
+                        var.pack_raw_data_to_petsc(array, sync=True)
                     finally:
-                        self.mesh._mesh_update_lock.release()
+                        var.mesh._mesh_update_lock.release()
                 else:
                     # Fallback if no lock exists
-                    self.pack_raw_data_to_petsc(array, sync=True)
+                    var.pack_raw_data_to_petsc(array, sync=True)
             finally:
                 # Clear recursion guard
-                self._in_flat_callback = False
+                var._in_flat_callback = False
 
         # Register the callback
         array_obj.add_callback(flat_data_update_callback)
@@ -2535,12 +2547,16 @@ def _create_canonical_data_array(self):
 
         # Create NDArray_With_Callback with proper shape and data
         from underworld3.utilities import NDArray_With_Callback
-
-        array_obj = NDArray_With_Callback(flat_petsc_data)
+        array_obj = NDArray_With_Callback(flat_petsc_data, owner=self)
 
         # Single canonical callback for PETSc synchronization
+
         def canonical_data_callback(array, change_context):
             """ONLY callback that handles PETSc synchronization - prevents conflicts"""
+            var = array.owner
+            if var is None:
+                return
+
             # Only act on data-changing operations
             data_changed = change_context.get("data_has_changed", True)
             if not data_changed:
@@ -2556,26 +2572,26 @@ def canonical_data_callback(array, change_context):
 
             canonical_array = np.atleast_2d(array)
 
-            if canonical_array.shape != (canonical_array.shape[0], self.num_components):
+            if canonical_array.shape != (canonical_array.shape[0], var.num_components):
                 # Only reshape if we actually need to
-                canonical_array = canonical_array.reshape(-1, self.num_components)
+                canonical_array = canonical_array.reshape(-1, var.num_components)
 
             # Skip updates during mesh coordinate changes to prevent corruption
-            if hasattr(self.mesh, "_mesh_update_lock"):
-                if not self.mesh._mesh_update_lock.acquire(blocking=False):
+            if hasattr(var.mesh, "_mesh_update_lock"):
+                if not var.mesh._mesh_update_lock.acquire(blocking=False):
                     return
                 try:
                     # STEP 1: Sync to PETSc using established method with correct shape
-                    self.pack_raw_data_to_petsc(canonical_array, sync=True)
+                    var.pack_raw_data_to_petsc(canonical_array, sync=True)
                 finally:
-                    self.mesh._mesh_update_lock.release()
+                    var.mesh._mesh_update_lock.release()
             else:
                 # Fallback if no lock exists
-                self.pack_raw_data_to_petsc(canonical_array, sync=True)
+                var.pack_raw_data_to_petsc(canonical_array, sync=True)
 
             # STEP 2: Handle variable-specific updates (extensible like SwarmVariable)
-            if hasattr(self, "_on_data_changed"):
-                self._on_data_changed()
+            if hasattr(var, "_on_data_changed"):
+                var._on_data_changed()
 
         array_obj.add_callback(canonical_data_callback)
         return array_obj
diff --git a/src/underworld3/function/expressions.py b/src/underworld3/function/expressions.py
index 2f24914b3..2a35b6312 100644
--- a/src/underworld3/function/expressions.py
+++ b/src/underworld3/function/expressions.py
@@ -659,8 +659,13 @@ def __new__(
 
         # Check both dicts for name collisions
         name_exists_persistent = name in UWexpression._expr_names
-        # Check ephemeral dict - need to look for any key starting with this name
-        name_exists_ephemeral = any(k[0] == name for k in UWexpression._ephemeral_expr_names)
+        # Check ephemeral dict - need to look for any key starting with this name.
+        # Snapshot keys via list(...) before iterating: the underlying dict can
+        # be mutated mid-iteration by weakref finalizers running asynchronously
+        # during cyclic GC, raising "dictionary changed size during iteration".
+        name_exists_ephemeral = any(
+            k[0] == name for k in list(UWexpression._ephemeral_expr_names)
+        )
 
         # Determine unique ID for disambiguation
         # When _unique_name_generation=True, ALWAYS use instance_no as _uw_id
diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py
index ff8dc593d..41da07950 100644
--- a/src/underworld3/swarm.py
+++ b/src/underworld3/swarm.py
@@ -393,17 +393,25 @@ def _create_variable_array(self, initial_data=None):
         # Single callback function (following swarm_update_callback pattern)
         def variable_update_callback(array, change_context):
             """Callback to sync variable changes back to PETSc (like swarm.points)"""
+            var = array.owner
+            if var is None:
+                # This guard handles cases where the array is accessed during
+                # object teardown (e.g. at application exit or mesh rebuilds),
+                # where the owning Python variable has already been garbage
+                # collected but the NDArray proxy still exists.
+                return
+
             # Only act on data-changing operations (following swarm.points pattern)
             data_changed = change_context.get("data_has_changed", True)
             if not data_changed:
                 return
 
             # Skip updates during coordinate changes to prevent corruption
-            if hasattr(self.swarm, "_migration_disabled") and self.swarm._migration_disabled:
+            if hasattr(var.swarm, "_migration_disabled") and var.swarm._migration_disabled:
                 return
 
             # Persist changes to PETSc (like swarm callback updates coordinates)
-            self.pack_uw_data_to_petsc(array, sync=True)
+            var.pack_uw_data_to_petsc(array, sync=True)
 
         # Register the callback (following swarm.points pattern)
         array_obj.add_callback(variable_update_callback)
@@ -4537,10 +4545,14 @@ def advection(
             #
             #
 
-        # Remove points no longer in the domain
+        # Re-route particles to their owning ranks and remove any that
+        # have genuinely left the domain. Use the default max_its so that
+        # boundary particles whose owner is the 2nd/3rd-closest centroid
+        # get reclaimed via the kdtree retry — max_its=1 here was an
+        # accidental regression that deleted boundary particles (issue #175,
+        # reported by @bknight1).
         self.migrate(
             delete_lost_points=True,
-            max_its=1,
         )
 
         return
diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index 8b2af5d88..5e26d4cc2 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -62,6 +62,7 @@
 from underworld3.cython.generic_solvers import SNES_MultiComponent
 from underworld3 import VarType
 import underworld3.timing as timing
+from underworld3.utilities import memprobe
 from underworld3.utilities._api_tools import (
     uw_object,
     SymbolicProperty,
@@ -1232,6 +1233,7 @@ def _create_stress_history_ddt(self, order=2):
         self.Unknowns.DFDt.enable_source_snapshot()
 
     @timing.routine_timer_decorator
+    @memprobe.instrument("Stokes.solve")
     def solve(
         self,
         zero_init_guess: bool = True,
@@ -3724,6 +3726,7 @@ def penalty(self, value):
         self._penalty.sym = value
 
     @timing.routine_timer_decorator
+    @memprobe.instrument("NavierStokes.solve")
     def solve(
         self,
         zero_init_guess: bool = True,
diff --git a/src/underworld3/utilities/__init__.py b/src/underworld3/utilities/__init__.py
index 495301e31..1731805b3 100644
--- a/src/underworld3/utilities/__init__.py
+++ b/src/underworld3/utilities/__init__.py
@@ -83,3 +83,4 @@ def _append_petsc_path():
 )
 
 from . import retention_curves
+from . import memprobe
diff --git a/src/underworld3/utilities/memprobe.py b/src/underworld3/utilities/memprobe.py
new file mode 100644
index 000000000..52e22860c
--- /dev/null
+++ b/src/underworld3/utilities/memprobe.py
@@ -0,0 +1,298 @@
+"""
+Lightweight memory-growth diagnostics.
+
+Designed to surface slow leaks in long parallel runs (HPC OOM-class bugs)
+without slowing normal use. Three signal sources:
+
+* **Process RSS** — current resident-set size via :mod:`psutil` if available,
+  falling back to ``/proc/self/statm`` on Linux, then to :mod:`resource`'s
+  high-water-mark RSS as a last resort. Current RSS is preferred so that
+  freed memory shows up as a negative delta.
+* **Live ``uw.kdtree.KDTree`` count** — UW's nanoflann wrapper isn't visible
+  to PETSc's object tracker, so it's a common silent-leak suspect.
+* **Per-class Python instance counts** for the ``underworld3`` package, gated
+  behind ``full=True`` because :func:`gc.get_objects` is slow at scale.
+
+For PETSc-side object/allocation tracking, use the runtime flags this module
+documents (``-log_view``, ``-malloc_dump``) — parsing PETSc's viewer output
+from Python is brittle and adds no signal beyond what those flags give.
+
+Activation
+----------
+- ``UW_MEMPROBE=1`` env var → enables instrumentation hooks at import time.
+- :func:`enable` / :func:`disable` for runtime toggling.
+- :func:`probe` context manager for ad-hoc snapshot+diff blocks.
+- :func:`instrument` decorator for per-method instrumentation; a no-op when
+  disabled (one attribute lookup + branch, sub-microsecond).
+
+Examples
+--------
+>>> import underworld3 as uw
+>>> from underworld3.utilities import memprobe
+>>>
+>>> with memprobe.probe("step 42"):
+...     do_one_step()           # diff is logged on exit
+>>>
+>>> @memprobe.instrument("stokes-solve")
+... def solve(...):
+...     ...
+
+For a long run, set ``UW_MEMPROBE=1`` and add a probe block per step. To pin
+PETSc-side leaks, run with ``-log_view -malloc_dump`` (or call
+:func:`dump_petsc_leaks_at_finalize`).
+"""
+from __future__ import annotations
+
+import functools
+import gc
+import os
+import sys
+from collections import Counter
+from contextlib import contextmanager
+from typing import Any
+
+__all__ = [
+    "ENABLED",
+    "enable",
+    "disable",
+    "snapshot",
+    "diff",
+    "probe",
+    "instrument",
+    "dump_petsc_leaks_at_finalize",
+    "format_diff",
+]
+
+ENABLED: bool = bool(int(os.environ.get("UW_MEMPROBE", "0")))
+
+
+def enable() -> None:
+    """Turn instrumentation on at runtime."""
+    global ENABLED
+    ENABLED = True
+
+
+def disable() -> None:
+    """Turn instrumentation off at runtime."""
+    global ENABLED
+    ENABLED = False
+
+
+def _rss_mb() -> float:
+    """Current resident-set size in MiB.
+
+    Prefers ``psutil.Process().memory_info().rss`` (true current RSS, drops
+    when memory is freed). Falls back to ``/proc/self/statm`` on Linux when
+    psutil isn't available, then to ``resource.ru_maxrss`` (peak/high-water
+    RSS — does not decrease) as a last resort. The fallback is good enough
+    for spotting growth, less good for spotting recovery.
+    """
+    try:
+        import psutil
+
+        return psutil.Process().memory_info().rss / (1024 * 1024)
+    except ImportError:
+        pass
+
+    if sys.platform.startswith("linux"):
+        try:
+            with open("/proc/self/statm", "r") as f:
+                # statm: size resident shared text lib data dt (in pages)
+                resident_pages = int(f.read().split()[1])
+            page_size = os.sysconf("SC_PAGESIZE")
+            return (resident_pages * page_size) / (1024 * 1024)
+        except (OSError, ValueError):
+            pass
+
+    # Last resort: peak RSS — does NOT decrease, so growth shows but
+    # recovery does not. ru_maxrss is KiB on Linux, bytes on macOS.
+    import resource
+
+    rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
+    if sys.platform == "darwin":
+        return rss / (1024 * 1024)
+    return rss / 1024
+
+
+def _kdtree_counts() -> dict[str, int]:
+    try:
+        from underworld3 import kdtree
+
+        return {
+            "live": kdtree.live_count(),
+            "total_constructed": kdtree.total_constructed(),
+        }
+    except (ImportError, AttributeError):
+        return {"live": -1, "total_constructed": -1}
+
+
+def _python_class_counts(prefix: str = "underworld3") -> Counter[str]:
+    """Count live Python objects whose class lives under ``prefix``.
+
+    Walks ``gc.get_objects()`` once. Slow at large counts (~10–100 ms for
+    typical UW runs). Only called when ``full=True`` is requested.
+    """
+    counts: Counter[str] = Counter()
+    for obj in gc.get_objects():
+        cls = type(obj)
+        mod = getattr(cls, "__module__", "") or ""
+        # Some objects expose a non-string descriptor as ``__module__``;
+        # skip those rather than crashing the walk.
+        if isinstance(mod, str) and mod.startswith(prefix):
+            counts[f"{mod}.{cls.__name__}"] += 1
+    return counts
+
+
+def snapshot(full: bool = False) -> dict[str, Any]:
+    """Capture a memory snapshot.
+
+    Parameters
+    ----------
+    full
+        If ``True``, also walk ``gc.get_objects()`` to count live Python
+        instances by class (slow). Defaults to ``False``.
+
+    Returns
+    -------
+    dict
+        Keys: ``rss_mb``, ``kdtree`` (mapping with ``live`` and
+        ``total_constructed``), ``py_classes`` (mapping; only present when
+        ``full=True``).
+    """
+    snap: dict[str, Any] = {
+        "rss_mb": _rss_mb(),
+        "kdtree": _kdtree_counts(),
+    }
+    if full:
+        snap["py_classes"] = dict(_python_class_counts())
+    return snap
+
+
+def diff(before: dict[str, Any], after: dict[str, Any]) -> dict[str, Any]:
+    """Compute a structured diff between two snapshots.
+
+    Returns deltas only (excludes anything unchanged or absent). Python
+    class entries are sorted by absolute change so the dominant suspects
+    appear first.
+    """
+    delta: dict[str, Any] = {}
+
+    drss = after["rss_mb"] - before["rss_mb"]
+    # Threshold below the format precision (0.01 MiB ≈ 10 KiB): smaller
+    # noise — inevitable when the source is RSS in KiB pages — would print
+    # as "+0.00 MiB" and defeat the "no change" fast path.
+    if abs(drss) >= 0.01:
+        delta["rss_mb"] = drss
+
+    kd_delta: dict[str, int] = {}
+    for key in ("live", "total_constructed"):
+        d = after["kdtree"].get(key, 0) - before["kdtree"].get(key, 0)
+        if d != 0:
+            kd_delta[key] = d
+    if kd_delta:
+        delta["kdtree"] = kd_delta
+
+    if "py_classes" in before and "py_classes" in after:
+        b = before["py_classes"]
+        a = after["py_classes"]
+        keys = set(b) | set(a)
+        py_delta = {k: a.get(k, 0) - b.get(k, 0) for k in keys}
+        py_delta = {k: v for k, v in py_delta.items() if v != 0}
+        if py_delta:
+            delta["py_classes"] = dict(
+                sorted(py_delta.items(), key=lambda kv: abs(kv[1]), reverse=True)
+            )
+
+    return delta
+
+
+def format_diff(label: str, delta: dict[str, Any]) -> str:
+    """Render a diff for stdout / logs. One line for trivial deltas, more if needed."""
+    if not delta:
+        return f"[memprobe] {label}: no change"
+
+    parts = [f"[memprobe] {label}:"]
+    if "rss_mb" in delta:
+        parts.append(f"  RSS {delta['rss_mb']:+.2f} MiB")
+    if "kdtree" in delta:
+        kd = delta["kdtree"]
+        chunk = ", ".join(f"{k} {v:+d}" for k, v in kd.items())
+        parts.append(f"  kdtree: {chunk}")
+    if "py_classes" in delta:
+        parts.append("  py_classes (top 10 by |Δ|):")
+        for cls, d in list(delta["py_classes"].items())[:10]:
+            parts.append(f"    {cls}: {d:+d}")
+    return "\n".join(parts)
+
+
+@contextmanager
+def probe(label: str, full: bool = False, *, emit=print):
+    """Context manager: snapshot before, snapshot after, log the diff.
+
+    Parameters
+    ----------
+    label
+        Short identifier for the block, included in the emitted diff.
+    full
+        Walk ``gc.get_objects()`` for per-class deltas (slow).
+    emit
+        Callable receiving the formatted diff string. Defaults to
+        :func:`print`. Use a logger or rank-aware writer for parallel runs.
+    """
+    if not ENABLED:
+        yield
+        return
+    before = snapshot(full=full)
+    try:
+        yield
+    finally:
+        after = snapshot(full=full)
+        emit(format_diff(label, diff(before, after)))
+
+
+def instrument(label: str, full: bool = False, *, emit=print):
+    """Decorator wrapping a method/function with :func:`probe`.
+
+    Fast-returns when :data:`ENABLED` is ``False`` (sub-microsecond overhead),
+    so it's safe to leave permanent decorations on hot paths.
+
+    Examples
+    --------
+    >>> @instrument("stokes-solve")
+    ... def solve(self, ...):
+    ...     ...
+    """
+    def decorator(func):
+        @functools.wraps(func)
+        def wrapper(*args, **kw):
+            if not ENABLED:
+                return func(*args, **kw)
+            with probe(label, full=full, emit=emit):
+                return func(*args, **kw)
+        return wrapper
+    return decorator
+
+
+def dump_petsc_leaks_at_finalize(filename: str | None = None) -> None:
+    """Ask PETSc to dump unfreed objects at finalize.
+
+    Sets the runtime options ``malloc_dump`` and ``objects_dump`` so PETSc
+    writes its leak report when the process tears down. When ``filename``
+    is given, also sets ``malloc_view`` to the same path for a more verbose
+    allocation dump.
+
+    Note that ``petsc4py``'s :class:`Options` API expects keys *without* the
+    leading ``-`` that you'd type on the command line — same convention used
+    elsewhere in this codebase.
+
+    Equivalent to running with ``-malloc_dump -objects_dump`` on the command
+    line. Useful when you can't change the launch command but can call this
+    from Python early in the run.
+    """
+    from petsc4py import PETSc
+
+    opts = PETSc.Options()
+    opts.setValue("malloc_dump", "")
+    opts.setValue("objects_dump", "")
+    if filename:
+        opts.setValue("malloc_view", filename)
diff --git a/tests/parallel/test_0765_swarm_advection_no_loss.py b/tests/parallel/test_0765_swarm_advection_no_loss.py
new file mode 100644
index 000000000..6a68d3f52
--- /dev/null
+++ b/tests/parallel/test_0765_swarm_advection_no_loss.py
@@ -0,0 +1,96 @@
+"""
+Regression test for swarm particle loss during advection across processor boundaries.
+
+Bug: ``Swarm.advection()`` finished with ``self.migrate(delete_lost_points=True,
+max_its=1)``. ``max_its=1`` only tries the *closest* domain centroid for each
+unclaimed particle, so any particle whose owning rank happened to be the
+2nd-or-3rd closest centroid (typical near a process boundary) failed
+``points_in_domain`` on its sole try and got deleted.
+
+Reproducer adapted from @bknight1's report on issue #175.
+
+See: https://github.com/underworldcode/underworld3/issues/175
+
+Run with:
+    mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_0765_swarm_advection_no_loss.py
+    mpirun -n 4 python -m pytest --with-mpi tests/parallel/test_0765_swarm_advection_no_loss.py
+"""
+
+import pytest
+import numpy as np
+import sympy
+import underworld3 as uw
+from mpi4py import MPI
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(60)]
+
+
+@pytest.mark.mpi(min_size=2)
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_advection_preserves_particle_count_order1():
+    """Order-1 advection across rank boundaries must not lose particles."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0),
+        maxCoords=(1.0, 1.0),
+        cellSize=0.1,
+    )
+
+    swarm = uw.swarm.Swarm(mesh)
+    var = uw.swarm.SwarmVariable("test_var", swarm, 1)
+    swarm.populate(fill_param=1)
+    var.data[...] = uw.mpi.rank
+
+    comm = uw.mpi.comm
+    initial_count = comm.allreduce(swarm.dm.getLocalSize(), op=MPI.SUM)
+
+    # Constant rightward velocity; dt=0.6 sweeps particles ~6 cells across the
+    # domain in one step. With max_its=1 in the post-advection migrate this
+    # used to drop boundary particles.
+    v_fn = sympy.Matrix([1.0, 0.0])
+    swarm.advection(v_fn, 0.6, order=1)
+
+    final_count = comm.allreduce(swarm.dm.getLocalSize(), op=MPI.SUM)
+
+    assert final_count == initial_count, (
+        f"Lost {initial_count - final_count} particles "
+        f"(initial={initial_count}, final={final_count})"
+    )
+
+
+@pytest.mark.mpi(min_size=2)
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_advection_preserves_particle_count_order2():
+    """Order-2 mid-point advection across rank boundaries must not lose particles.
+
+    Uses solid-body rotation about the domain centre so particles stay inside
+    the unit box for the whole rotation. This isolates rank-boundary loss
+    from genuine domain-exit, while still forcing every particle across
+    multiple processor cuts over many substeps.
+    """
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0),
+        maxCoords=(1.0, 1.0),
+        cellSize=0.1,
+    )
+
+    swarm = uw.swarm.Swarm(mesh)
+    swarm.populate(fill_param=1)
+
+    comm = uw.mpi.comm
+    initial_count = comm.allreduce(swarm.dm.getLocalSize(), op=MPI.SUM)
+
+    # Solid-body rotation about (0.5, 0.5) — radius <= sqrt(0.5) so nothing
+    # can escape the unit box. omega=2*pi makes one revolution in dt=1.
+    x, y = mesh.X
+    omega = 2 * sympy.pi
+    v_fn = sympy.Matrix([-omega * (y - 0.5), omega * (x - 0.5)])
+    swarm.advection(v_fn, 0.5, order=2, step_limit=True)
+
+    final_count = comm.allreduce(swarm.dm.getLocalSize(), op=MPI.SUM)
+
+    assert final_count == initial_count, (
+        f"Lost {initial_count - final_count} particles "
+        f"(initial={initial_count}, final={final_count})"
+    )
diff --git a/tests/test_0780_memprobe.py b/tests/test_0780_memprobe.py
new file mode 100644
index 000000000..d1b4b5414
--- /dev/null
+++ b/tests/test_0780_memprobe.py
@@ -0,0 +1,153 @@
+"""
+Smoke tests for the memprobe diagnostic module.
+
+These tests check the *plumbing*: instrumentation enable/disable, snapshot
+shape, KDTree live-count tracking, decorator no-op-when-disabled, and
+context-manager diff emission. They do not validate that any real leak is
+or isn't present — that's what users do with the tool, not what the tool
+itself can verify.
+"""
+import gc
+import numpy as np
+import pytest
+
+import underworld3 as uw
+from underworld3.utilities import memprobe
+
+
+@pytest.fixture(autouse=True)
+def _restore_enabled_flag():
+    """Each test starts with memprobe disabled regardless of env."""
+    saved = memprobe.ENABLED
+    memprobe.disable()
+    yield
+    memprobe.ENABLED = saved
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_snapshot_shape():
+    """Snapshot must contain rss_mb (float) and kdtree (dict with two keys)."""
+    snap = memprobe.snapshot()
+
+    assert isinstance(snap["rss_mb"], float)
+    assert snap["rss_mb"] > 0
+    assert set(snap["kdtree"].keys()) == {"live", "total_constructed"}
+
+    # full=False should NOT walk gc
+    assert "py_classes" not in snap
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_snapshot_full_includes_py_classes():
+    snap = memprobe.snapshot(full=True)
+    assert "py_classes" in snap
+    # Sanity: gc walked something
+    assert isinstance(snap["py_classes"], dict)
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_kdtree_live_count_tracks_construction_destruction():
+    """KDTree __cinit__/__dealloc__ must update the live-instance counter.
+
+    The test asserts only on deltas around the operation under test, never
+    on absolute counts, because earlier tests in the same pytest session
+    can leave KDTree references alive that the cyclic GC may collect at
+    any time, shifting the baseline.
+    """
+    pts = np.random.random((20, 2))
+
+    # Flush any pending cyclic-GC clean-ups so the baseline doesn't shift
+    # under us between snapshots.
+    gc.collect()
+    before_live = uw.kdtree.live_count()
+    before_total = uw.kdtree.total_constructed()
+
+    tree = uw.kdtree.KDTree(pts)
+    assert uw.kdtree.live_count() - before_live == 1
+    assert uw.kdtree.total_constructed() - before_total == 1
+
+    after_create_live = uw.kdtree.live_count()
+
+    del tree
+    gc.collect()
+
+    # The only thing we care about: this construction/destruction pair
+    # produced a +1/-1 swing relative to its immediate before/after, and
+    # total_constructed never went backwards.
+    assert uw.kdtree.live_count() - after_create_live == -1
+    assert uw.kdtree.total_constructed() >= before_total + 1
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_diff_reports_kdtree_growth():
+    pts = np.random.random((20, 2))
+
+    before = memprobe.snapshot()
+    tree = uw.kdtree.KDTree(pts)  # noqa: F841 — kept alive across the diff
+    after = memprobe.snapshot()
+
+    delta = memprobe.diff(before, after)
+    assert delta["kdtree"]["live"] == 1
+    assert delta["kdtree"]["total_constructed"] == 1
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_diff_drops_unchanged_keys():
+    snap = memprobe.snapshot()
+    delta = memprobe.diff(snap, snap)
+    # Identical snapshots produce no deltas
+    assert "kdtree" not in delta
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_probe_context_emits_when_enabled(capsys):
+    memprobe.enable()
+    with memprobe.probe("test-block"):
+        _ = uw.kdtree.KDTree(np.random.random((10, 2)))
+
+    captured = capsys.readouterr()
+    assert "[memprobe] test-block" in captured.out
+    assert "kdtree" in captured.out
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_probe_context_silent_when_disabled(capsys):
+    assert memprobe.ENABLED is False
+    with memprobe.probe("silent"):
+        _ = uw.kdtree.KDTree(np.random.random((10, 2)))
+
+    captured = capsys.readouterr()
+    assert captured.out == ""
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_instrument_decorator_no_op_when_disabled(capsys):
+    """Decorator must not emit (or noticeably slow) when ENABLED is False."""
+    @memprobe.instrument("test-fn")
+    def f(x):
+        return x * 2
+
+    assert f(3) == 6
+    captured = capsys.readouterr()
+    assert captured.out == ""
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_instrument_decorator_emits_when_enabled(capsys):
+    @memprobe.instrument("test-fn")
+    def f(x):
+        return x * 2
+
+    memprobe.enable()
+    assert f(3) == 6
+    captured = capsys.readouterr()
+    assert "[memprobe] test-fn" in captured.out
diff --git a/tests/test_0852_swarm_integration_statistics.py b/tests/test_0852_swarm_integration_statistics.py
index f2ef4aef1..2606c583f 100644
--- a/tests/test_0852_swarm_integration_statistics.py
+++ b/tests/test_0852_swarm_integration_statistics.py
@@ -122,7 +122,11 @@ def test_clustered_swarm_shows_difference(self):
         # Arithmetic mean: biased toward left half values
         # (75% particles at x ≈ 0.25) + (25% particles at x ≈ 0.75)
         # ≈ 0.75 * (1 + 0.25) + 0.25 * (1 + 0.75) = 1.375
-        arithmetic_mean = s_var.array.mean()
+        # 
+        # CRITICAL MPI FIX: Use global swarm statistics for arithmetic mean
+        # to ensure comparison with integration mean (which is always global)
+        # works correctly on all ranks.
+        arithmetic_mean = s_var.global_sum() / s_var.global_size()
 
         # Integration-based mean: weights by volume equally
         # ∫∫(1 + x) dA = [x + 0.5x²]₀¹ × 1 = (1 + 0.5) = 1.5