Foolproof mesh-coordinate mutation: capability gate + deform() + SL-field CARRY transfer

src/underworld3/discretisation/discretisation_mesh.py

@@ -952,6 +952,18 @@ def mesh_update_callback(array, change_context):
         # discretisation/remesh.py.
         self._remesh_hooks = []
 
+        # Capability gate for coordinate mutation. `_deform_mesh` (the raw
+        # primitive that moves nodes WITHOUT the field/SL-history transfer)
+        # is only permitted inside a sanctioned context: a remesh transaction
+        # (`_in_remesh_transfer`, set by remesh_with_field_transfer) or a
+        # `_coord_mutation()` scope (depth>0, opened by `deform`,
+        # `ephemeral_coords`, and trusted internal movers). Outside those,
+        # a bare call on a mesh that already carries variables/history raises
+        # with a pointer to the public methods — so the field/history transfer
+        # can never be silently skipped. See `deform`/`ephemeral_coords`.
+        self._in_remesh_transfer = getattr(self, "_in_remesh_transfer", False)
+        self._coord_mutation_depth = 0
+
         self._evaluation_hash = None
         self._evaluation_interpolated_results = None
         self._dm_initialized = False
@@ -1034,7 +1046,11 @@ class ElementInfo:
                         f"{float(numpy.sqrt(_d2.max())):.2e}); the reloaded "
                         "reference geometry does not match the saved fine."
                     )
-                self._deform_mesh(_def_canon[_idx].reshape(-1, _cdim))
+                # Restoring saved deformed geometry; fields are reloaded
+                # separately from the checkpoint, so this is a sanctioned
+                # internal coordinate move (no live-state transfer needed).
+                with self._coord_mutation():
+                    self._deform_mesh(_def_canon[_idx].reshape(-1, _cdim))
 
         if verbose and uw.mpi.rank == 0:
             print(
@@ -1930,7 +1946,10 @@ def sync_coordinates_from_parent(self):
         new_sub_coords = numpy.array(self.X.coords)
         new_sub_coords[sub_rows] = self.parent.X.coords[parent_rows]
 
-        self._deform_mesh(new_sub_coords)
+        # Submesh follows its parent's geometry; this is a sanctioned
+        # internal coordinate move (the parent owns the transfer).
+        with self._coord_mutation():
+            self._deform_mesh(new_sub_coords)
         self._parent_mesh_version = self.parent._mesh_version
 
     def _re_extract_from_parent(self, verbose=False):
@@ -2759,6 +2778,116 @@ def unregister_remesh_hook(self, op):
             kept.append(r)
         self._remesh_hooks = kept
 
+    # ------------------------------------------------------------------
+    # Coordinate-mutation capability gate + sanctioned entry points
+    # ------------------------------------------------------------------
+    def _assert_coord_mutation_allowed(self):
+        """Guard for :meth:`_deform_mesh`.
+
+        Moving coordinates with the raw primitive skips the field /
+        SL-DDt-history transfer. That is only safe (a) before the mesh
+        carries any variables/history (construction, restart-before-solvers)
+        or (b) inside a sanctioned scope — a remesh transaction
+        (``_in_remesh_transfer``) or a ``_coord_mutation()`` scope opened by
+        :meth:`deform`, :meth:`ephemeral_coords`, or a trusted internal mover.
+        Outside those, on a live mesh, raise with a pointer to the public API.
+        """
+        allowed = (getattr(self, "_in_remesh_transfer", False)
+                   or getattr(self, "_coord_mutation_depth", 0) > 0)
+        if allowed:
+            return
+        has_state = bool(getattr(self, "vars", None)) or bool(
+            getattr(self, "_remesh_hooks", None))
+        if not has_state:
+            return
+        raise RuntimeError(
+            "Mesh._deform_mesh() is an internal primitive — it moves nodes "
+            "WITHOUT transferring fields or solver/DDt history onto the new "
+            "layout, which corrupts the solution. This mesh already carries "
+            "variables, so a direct call (or an in-place write to "
+            "`mesh.X.coords`) is rejected.\n"
+            "  Use instead:\n"
+            "    • mesh.deform(new_coords, dt=…)        — impose an arbitrary "
+            "node displacement (free surface / prescribed motion)\n"
+            "    • mesh.adapt(metric) / mesh.OT_adapt(field)\n"
+            "    • uw.meshing.smooth_mesh_interior(…) / uw.meshing.follow_metric(…)\n"
+            "  These route the field + SL/DDt-history transfer "
+            "(remesh_with_field_transfer). For trusted scheme-internal trial "
+            "meshes that will be discarded, use `with mesh.ephemeral_coords(): …`."
+        )
+
+    @contextmanager
+    def _coord_mutation(self):
+        """Internal: sanction direct ``_deform_mesh`` calls within this scope.
+
+        Re-entrant. Does NOT itself transfer fields — callers either are the
+        transfer transaction (REMAP + ``on_remesh`` already run by
+        ``remesh_with_field_transfer``), restore saved state separately, or
+        intend an ephemeral trial (see :meth:`ephemeral_coords`).
+        """
+        self._coord_mutation_depth += 1
+        try:
+            yield
+        finally:
+            self._coord_mutation_depth -= 1
+
+    @contextmanager
+    def ephemeral_coords(self):
+        """Trusted scheme-internal trial coordinate moves, restored on exit.
+
+        For schemes (e.g. RK4 surface stages) that probe trial meshes to
+        compute velocities/increments and must NOT commit a transfer. The
+        coordinates are snapshotted on enter and restored on exit, so the
+        intermediate meshes are genuinely ephemeral — only the final
+        committed move (via :meth:`deform`) updates fields + history.
+        """
+        saved = numpy.asarray(self.X.coords).copy()
+        self._coord_mutation_depth += 1
+        try:
+            yield
+        finally:
+            try:
+                self._deform_mesh(saved)
+            finally:
+                self._coord_mutation_depth -= 1
+
+    def deform(self, new_coords, *, dt=None, zero=None, verbose=False):
+        """Move the mesh to ``new_coords``, transferring all fields + history.
+
+        The public, foolproof way to impose an arbitrary node displacement
+        (free surface, prescribed mesh motion). Wraps the remesh transaction
+        so that REMAP variables are re-interpolated onto the new layout and
+        every registered ``on_remesh`` hook fires — in particular the
+        SemiLagrangian DDt's coherent ALE (carry its history stack + a
+        one-step ``v_mesh = Δx/dt`` correction consumed by the next solve).
+
+        Parameters
+        ----------
+        new_coords : ndarray
+            Target vertex coordinates (shape of ``mesh.X.coords``).
+        dt : float, optional
+            Timestep for the ALE ``v_mesh = Δx/dt`` correction. Required when
+            SemiLagrangian history is present and the move is advective (a
+            free-surface step); omit for a pure geometric re-mesh.
+        zero : list of MeshVariable, optional
+            Variables to zero after the move (e.g. V, P for a cold restart).
+        verbose : bool
+
+        Returns
+        -------
+        bool
+            True if the mesh moved (geometry changed), False for a no-op.
+        """
+        from underworld3.discretisation.remesh import remesh_with_field_transfer
+
+        _nc = numpy.asarray(new_coords)
+
+        def _do_move():
+            self._deform_mesh(_nc)
+
+        return remesh_with_field_transfer(
+            self, _do_move, dt=dt, extra_zero=zero, verbose=verbose)
+
     def _deform_mesh(self, new_coords: numpy.ndarray, verbose=False,
                      active_vars=None):
         """
@@ -2777,8 +2906,19 @@ def _deform_mesh(self, new_coords: numpy.ndarray, verbose=False,
         recompute n_outer×); the wrapper does a full recompute once at
         sweep exit by calling ``_deform_mesh`` again with
         ``active_vars=None``.
+
+        .. warning::
+
+           This is the **internal** coordinate-mutation primitive. It moves
+           nodes WITHOUT transferring fields or solver/DDt history onto the
+           new layout. It may only be called inside a sanctioned coordinate-
+           mutation scope (see :meth:`deform`, :meth:`ephemeral_coords`, and
+           ``remesh_with_field_transfer``). A bare call on a mesh that already
+           carries variables raises — use the public methods instead.
         """
 
+        self._assert_coord_mutation_allowed()
+
         with self._mesh_update_lock:
             coord_vec = self.dm.getCoordinatesLocal()
             coords = coord_vec.array.reshape(-1, self.cdim)
@@ -3712,7 +3852,10 @@ def apply_snapshot_payload(self, payload: dict) -> None:
                 f"({coords.shape} vs current {expected_shape}); programming "
                 f"error since _mesh_version matched"
             )
-        self._deform_mesh(coords)
+        # Snapshot restore: variables are reloaded just below, so this is a
+        # sanctioned internal coordinate move (no live-state transfer needed).
+        with self._coord_mutation():
+            self._deform_mesh(coords)
 
         current_vars = {var.clean_name: var for var in self.vars.values()}
         for var_clean_name, saved_array in payload["vars"].items():

src/underworld3/systems/solvers.py

@@ -3365,6 +3365,22 @@ def __init__(
             theta=theta,
         )
 
+        # Phase-2 remesh: the advected field must ride the mesh (CARRY)
+        # coherently with its own semi-Lagrangian history (psi_star,
+        # already CARRY + managed). Under the default REMAP it would be
+        # geometrically re-interpolated onto the new node positions AND
+        # then have v_mesh = Δx/dt subtracted in the next trace-back —
+        # a DOUBLE compensation for the mesh motion, inconsistent with
+        # the once-compensated CARRY'd history. Stamping it CARRY +
+        # managed-by-DuDt makes deform()/adapt transfer it on the single
+        # ALE trace-back path (and REMAP it correctly on an OT opt-out
+        # reset). Only meaningful when DuDt traces back (SemiLagrangian);
+        # Eulerian/Lagrangian fields keep the default policy.
+        if isinstance(self.Unknowns.DuDt, SemiLagrangian_DDt):
+            from underworld3.discretisation.remesh import RemeshPolicy
+            self.u.remesh_policy = RemeshPolicy.CARRY
+            self.u._remesh_managed_by = self.Unknowns.DuDt
+
         return
 
     @property

tests/test_0820_deform_mesh_solver_rebuild_regression.py

@@ -48,7 +48,10 @@ def deform_with_amplitude(amp, var_name):
         diff.solve()
         disp = np.zeros((mesh.X.coords.shape[0], mesh.dim))
         disp[:, -1] = uw.function.evaluate(Dz.sym[0], mesh.X.coords)[:, 0, 0]
-        mesh._deform_mesh(mesh.X.coords + disp)
+        # This regression deliberately tests the gated internal primitive's
+        # solver-rebuild (is_setup=False) contract — open the sanctioned scope.
+        with mesh._coord_mutation():
+            mesh._deform_mesh(mesh.X.coords + disp)
 
     stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
     stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel

tests/test_0825_deform_mesh_cache_invalidation.py

@@ -57,7 +57,11 @@ def test_topology_version_bumps_on_deform():
     coords = np.asarray(mesh.X.coords)
     perturbed = coords.copy()
     perturbed[:, 0] += 1.0e-3 * coords[:, 1]
-    mesh._deform_mesh(perturbed)
+    # _deform_mesh is the internal primitive (gated against bare use on a
+    # live mesh); this test deliberately exercises it, so open the
+    # sanctioned _coord_mutation scope.
+    with mesh._coord_mutation():
+        mesh._deform_mesh(perturbed)
     assert mesh._topology_version > before, (
         "_topology_version should be incremented by _deform_mesh")
 
@@ -69,7 +73,8 @@ def test_evaluation_hash_invalidated_on_deform():
     coords = np.asarray(mesh.X.coords)
     perturbed = coords.copy()
     perturbed[:, 1] += 5.0e-4 * coords[:, 0]
-    mesh._deform_mesh(perturbed)
+    with mesh._coord_mutation():
+        mesh._deform_mesh(perturbed)
     assert mesh._evaluation_hash is None
     assert mesh._evaluation_interpolated_results is None
 
@@ -83,7 +88,8 @@ def test_dminterpolation_cache_invalidated_on_deform():
     coords = np.asarray(mesh.X.coords)
     perturbed = coords.copy()
     perturbed[:, 1] += 1.0e-3 * coords[:, 0]
-    mesh._deform_mesh(perturbed)
+    with mesh._coord_mutation():
+        mesh._deform_mesh(perturbed)
     n_after = len(getattr(
         mesh._dminterpolation_cache, "_cache", {}))
     assert n_after == 0, (

tests/test_0850_mesh_smoothing.py

@@ -143,7 +143,9 @@ def test_aspect_ratio_improves_on_perturbed_mesh(self):
         is_interior = ~is_bnd
         coords[is_interior] += 0.04 * rng.standard_normal(
             (int(is_interior.sum()), coords.shape[1]))
-        mesh._deform_mesh(coords)
+        # deliberate direct use of the gated internal primitive
+        with mesh._coord_mutation():
+            mesh._deform_mesh(coords)
         ar_before = _min_aspect_ratio(mesh)
         smooth_mesh_interior(mesh, n_iters=10, alpha=0.5)
         ar_after = _min_aspect_ratio(mesh)