Adaptivity round 3: the 3D mover+adapt combination, spherical domains, and boundary snapping

docs/developer/design/ADAPTIVITY_3D_SPHERICAL_2026-07.md

@@ -626,6 +626,26 @@ worked example (3D box or annulus fault: metric-driven redistribution of the
 base + NVB band on the fault + Stokes FMG + advection-diffusion with field
 transfer across re-adaptation), plus a short how-to in `docs/advanced/`.
 
+**Round-3a status: 3D combination WORKS (2026-07-24, worktree
+`adaptivity-round3`).** First run stalled the refinement drain — and the
+cause was a PETSc label-semantics trap, not the driver: `DMLabelSetValue`
+does not remove a point from its previous stratum, so RE-seeding the
+tagged refinement state on the moved base left every cell in two strata
+and readers got the *old* (unmoved-order) state back — the driver was
+silently running the moved geometry with a stale seed. (The driver
+itself was first cleared by exhaustive oracle sweeps: all 72 single-tet
+states × 3 passes, and all 3,456 engine-valid two-tet state pairs with
+single-cell marking — native vs Python engine, all equal.)
+`write_tagged_state_label` now destroys and recreates the label;
+regression test in `test_0840`. Result on the dipping-fault box (ml=2):
+mover+adapt gives **33,570 cells vs 36,040 adapt-only**, slightly finer
+on-fault (med h 0.0282 vs 0.0292) with the over-refinement fraction
+down from 46% to 39% — the same budget-spends-better behaviour as 2D.
+MG on the combined child: 3 iterations, exact solution to 9e-10.
+Renders: `~/+Simulations/nvb_3d_adapt_evaluation/combo3d_*.png`. The
+worked example + `docs/advanced/` how-to remain the landing
+deliverable.
+
 ### Round 3b — spherical geometry + curved/internal-boundary clean-ups
 
 - **MMPDE on shells**: `SphericalShell` (solid, `dim==cdim==3`) is a valid
@@ -651,6 +671,44 @@ transfer across re-adaptation), plus a short how-to in `docs/advanced/`.
   path is boundary-label-driven). Validation: shell response benchmark
   (test_1064 pattern) on an NVB child.
 
+**Round-3b probe results (2026-07-24):**
+
+* *MMPDE on `SphericalShell`: validated, zero new code.* With the
+  registered radial `BoundingSurface`s, a boundary-layer metric toward
+  the inner sphere gives zero folds and holds BOTH spheres to machine
+  precision (radius drift ≤ 2e-16) while the radial cell-size profile
+  grades to ratio ~1.32 (better than the box's 1.12 — the boundary-layer
+  metric shape suits the mover well).
+* *The chord error is real and measured*: adapted annulus children carry
+  an outer-boundary radius error of 1.8e-3 at base h≈0.125
+  (the h²/8R chord sag), frozen at base resolution regardless of depth.
+* *Snap is safe on every gate*: projecting the child's boundary vertices
+  onto the analytic circles gives zero folds, radius error → 2e-16, and
+  the custom-P MG on the snapped child is untouched (2 iterations,
+  solution identical to the GAMG reference) — the coordinate-built
+  transfers tolerate the slightly non-nested snapped boundary exactly as
+  designed. **The chord-vs-snap ruling is now a pure decision** — no
+  technical blocker either way.
+* *A real MG gap found and fixed on the way*: the mesh-owned FMG
+  auto-pickup rejected any square finest transfer as "no coarsening" —
+  but a generation whose new vertices all land on a Dirichlet boundary
+  adds only constrained dofs, making the free-dof counts of the last two
+  levels coincide. Boundary-focused metrics (what curved domains invite)
+  hit this every time, so every annulus boundary-layer adapt was
+  silently falling back to the default preconditioner. Guard relaxed to
+  genuine row/operator mismatches; annulus children now solve with
+  custom MG in 2 iterations.
+* *Internal interfaces resolved (maintainer discussion, 2026-07-24):*
+  refinement always preserved embedded interfaces topologically
+  (conforming bisection + label carry), but they were chord-frozen and
+  invisible to the snap. The Internal circle/sphere of the
+  *InternalBoundary meshes now registers as a radial surface flagged
+  `interior=True`: adapt() snaps refinement onto the true interface
+  radius, while the movers keep interface nodes FULLY pinned (no normal
+  or tangential motion, even with slip_surfaces=True) — interface
+  motion is physics-owned. Tangential slide along an interface is a
+  one-flag change if ever wanted.
+
 ## Effort and risk, honestly
 
 | round | new-code size | risk | notes |

src/underworld3/discretisation/discretisation_mesh.py

@@ -3136,7 +3136,14 @@ def boundary_slip(self, slip_spec=True, reference_coords=None,
         # labels (dict ``False``) still slide-without-restore regardless of
         # kind (handled in ``project`` below). A label with no boundary facets
         # at all stays unusable → its vertices pin (the safe default).
-        surf = dict(self.bounding_surfaces)
+        # INTERIOR surfaces (embedded interfaces, e.g. the Internal circle/
+        # sphere of the *InternalBoundary meshes) are never slip-eligible:
+        # interface motion is physics-owned (deform / free-surface
+        # machinery), so their nodes fall through to the pinned default —
+        # exactly as when no surface was registered. Their registration
+        # exists for adapt()'s refinement snapping.
+        surf = {lab: s for lab, s in dict(self.bounding_surfaces).items()
+                if not getattr(s, "interior", False)}
         unreg = [lab for lab in slip_labels if lab not in surf]
         if unreg:
             facets, _opp = _boundary_facets(self, cdim)
@@ -6516,6 +6523,94 @@ def cell_geometry(dm):
                 vol_scaled = numpy.sqrt(numpy.abs(numpy.linalg.det(G)))
             return X.mean(axis=1), vol_scaled ** (1.0 / dim), cs
 
+        # Analytic boundary surfaces to snap each generation onto. New
+        # boundary vertices are CHORD midpoints, so without snapping the
+        # boundary geometry of a curved domain stays frozen at base
+        # resolution no matter how deep the refinement. Per the 2026-07
+        # round-3b ruling, EVERY generation snaps (not just the returned
+        # child): each intermediate level is a valid mesh in its own
+        # right — extractable for solvers — and the metric marks on true
+        # geometry. On plane surfaces (boxes) the chord midpoint already
+        # lies in the plane, so the snap is exactly a no-op and the flat
+        # confluence gates are untouched.
+        # NOTE: surfaces must be disjoint or intersect only where each
+        # projection is an exact no-op (concentric radial pairs; orthogonal
+        # planes, whose corner vertices are fixed points of both). Sequential
+        # restore does not converge to the intersection of two CURVED
+        # surfaces — junction handling as in boundary_slip is a follow-up
+        # for the day such a mesh registers surfaces.
+        _snap_surfs = [s for s in dict(self.bounding_surfaces).values()
+                       if getattr(s, "kind", None) in ("radial", "plane")
+                       and not getattr(s, "is_free", False)]
+
+        def snap_level_boundaries(dm):
+            if not _snap_surfs:
+                return
+            from underworld3.meshing.smoothing import _pinned_mask
+
+            def sync_mask(mask):
+                # In 3D a boundary-edge midpoint can live on a rank whose
+                # only cells containing it are INTERIOR members of the
+                # edge's star — that rank sees no labelled face and would
+                # skip the snap the face-owning rank applies, leaving one
+                # global vertex with two coordinates. Reduce the mask over
+                # the point SF (ADD ≡ logical-or here) so every rank
+                # holding the vertex agrees; pre-snap coordinates are
+                # rank-consistent and restore is a pure function of them,
+                # so a consistent mask gives a consistent snap.
+                if uw.mpi.size == 1:
+                    return mask
+                from underworld3.meshing.smoothing.graph import (
+                    _build_scalar_dm)
+                dm_s = _build_scalar_dm(dm)
+                lvec = dm_s.createLocalVector()
+                gvec = dm_s.createGlobalVector()
+                lvec.array[:] = mask.astype(float)
+                dm_s.localToGlobal(lvec, gvec, addv=True)
+                dm_s.globalToLocal(gvec, lvec)
+                out = numpy.asarray(lvec.array) > 0.5
+                lvec.destroy(); gvec.destroy(); dm_s.destroy()
+                return out
+
+            vec = dm.getCoordinatesLocal()
+            arr = vec.array.reshape(-1, self.cdim)
+            pre = arr.copy()
+            snapped_any = numpy.zeros(arr.shape[0], dtype=bool)
+            for s in _snap_surfs:
+                mask = sync_mask(_pinned_mask(dm, (s.label,)))
+                if mask.any():
+                    arr[mask] = s.restore(arr[mask])
+                    snapped_any |= mask
+            dm.setCoordinatesLocal(vec)
+            if snapped_any.any():
+                # A snap moves boundary vertices by the chord sagitta
+                # (~h²/8R); on a base coarse enough that h ≈ R this could
+                # invert a sliver silently. Fail loudly instead: no cell
+                # incident to a snapped vertex may flip orientation.
+                cs_, ce_ = dm.getHeightStratum(0)
+                vS_, vE_ = dm.getDepthStratum(0)
+                flipped = 0
+                for c in range(cs_, ce_):
+                    vs = [p - vS_ for p in dm.getTransitiveClosure(c)[0]
+                          if vS_ <= p < vE_]
+                    if not any(snapped_any[v] for v in vs):
+                        continue
+                    e0 = pre[vs[1:]] - pre[vs[0]]
+                    e1 = arr[vs[1:]] - arr[vs[0]]
+                    if (numpy.sign(numpy.linalg.det(e0))
+                            != numpy.sign(numpy.linalg.det(e1))):
+                        flipped += 1
+                if uw.mpi.size > 1:
+                    from mpi4py import MPI as _MPI
+                    flipped = uw.mpi.comm.allreduce(flipped, op=_MPI.SUM)
+                if flipped:
+                    raise RuntimeError(
+                        f"adapt: snapping boundary vertices to the analytic "
+                        f"surfaces inverted {flipped} cell(s) — the base mesh "
+                        "is too coarse for the boundary curvature (chord "
+                        "sagitta ~ cell size). Refine the base mesh or adapt "
+                        "without registered bounding surfaces.")
+
         # Refine from the mesh's CURRENT geometry. Node redistribution
         # (redistribute_nodes) moves mesh.dm's coordinates while the static
         # base hierarchy keeps the originals — without this carry, an adapt
@@ -6588,6 +6683,7 @@ def cell_geometry(dm):
                 for cidx in marked:
                     lab.setValue(cidx, DM_ADAPT_REFINE)
                 current_dm = _nvbx.refine(d, "adapt")
+                snap_level_boundaries(current_dm)
                 level_dms.append(current_dm)
                 if verbose:
                     fs, fe = current_dm.getHeightStratum(0)
@@ -6628,6 +6724,24 @@ def cell_geometry(dm):
                 marked = [int(cids[j]) for j in sel]
                 markers_per_level.append(marked)
                 nvb.refine(set(marked))
+                # Snap INSIDE the engine: its own coordinates feed the next
+                # generation's marking AND the next midpoints, so snapping
+                # only the exported DM would leave the engine's geometry on
+                # the chords.
+                if _snap_surfs:
+                    _val2surf = {v: s for s in _snap_surfs
+                                 for (nm, v) in carry if nm == s.label}
+                    _sv = {}
+                    for _fkey, _val in nvb.facet_label.items():
+                        _s = _val2surf.get(_val)
+                        if _s is not None:
+                            _sv.setdefault(id(_s), (_s, set()))[1].update(_fkey)
+                    for _s, _vs in _sv.values():
+                        _idx = numpy.fromiter(_vs, dtype=numpy.int64)
+                        _snapped = _s.restore(
+                            numpy.array([nvb.coords[i] for i in _idx]))
+                        for _k, _i in enumerate(_idx):
+                            nvb.coords[_i] = _snapped[_k]
                 level_dms.append(nvb.to_dm(boundaries=carry, regions=rcarry,
                                            comm=self.dm.comm))
                 if verbose:
@@ -6664,6 +6778,7 @@ def cell_geometry(dm):
                     uw.pprint(0, f"[adapt] level {level}: refining {len(cell_ids)} "
                                  f"of {ncells} cells")
                 current_dm = custom_mg.sbr_refine(current_dm, cell_ids)
+                snap_level_boundaries(current_dm)
                 level_dms.append(current_dm)
 
         if verbose:

src/underworld3/meshing/annulus.py

@@ -1141,6 +1141,7 @@ def AnnulusInternalBoundary(
     degree: int = 1,
     qdegree: int = 2,
     filename=None,
+    refinement=None,
     gmsh_verbosity=0,
     verbose=False,
 ):
@@ -1403,6 +1404,7 @@ def annulus_internal_mesh_refinement_callback(dm):
         boundaries=boundaries,
         boundary_normals=None,
         coordinate_system_type=CoordinateSystemType.CYLINDRICAL2D,
+        refinement=refinement,
         refinement_callback=annulus_internal_mesh_refinement_callback,
         return_coords_to_bounds=annulus_internal_return_coords_to_bounds,
         verbose=verbose,
@@ -1421,6 +1423,17 @@ class boundary_normals(Enum):
     x, y = new_mesh.X
     new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])]
 
+    # Radial bounding surfaces: Upper/Lower slip+snap as usual; the embedded
+    # Internal circle is INTERIOR — adapt() snaps refinement onto its true
+    # radius, but the movers keep its nodes pinned (interface motion is
+    # physics-owned, 2026-07 round-3b ruling).
+    from underworld3.meshing.bounding_surface import register_radial_surfaces
+    register_radial_surfaces(
+        new_mesh, centre=(0.0, 0.0),
+        label_radius={"Upper": radiusOuter, "Lower": radiusInner,
+                      "Internal": radiusInternal},
+        interior=("Internal",))
+
     return new_mesh
 
 
@@ -1687,6 +1700,7 @@ def annulus_internal_mesh_refinement_callback(dm):
         boundaries=boundaries,
         boundary_normals=None,
         coordinate_system_type=CoordinateSystemType.CYLINDRICAL2D,
+        refinement=refinement,
         refinement_callback=annulus_internal_mesh_refinement_callback,
         return_coords_to_bounds=annulus_internal_return_coords_to_bounds,
         verbose=verbose,

src/underworld3/meshing/bounding_surface.py

@@ -61,14 +61,21 @@ class BoundingSurface:
     """
 
     def __init__(self, mesh, label, kind, *, centre=None, radius=None,
-                 point=None, normal=None, reference_facets=None, is_free=False):
+                 point=None, normal=None, reference_facets=None, is_free=False,
+                 interior=False):
         if kind not in _VALID_KINDS:
             raise ValueError(
                 f"BoundingSurface kind must be one of {_VALID_KINDS}; got {kind!r}")
         self._mesh = mesh
         self.label = str(label)
         self.kind = kind
         self.is_free = bool(is_free) or kind == "free"
+        # An INTERIOR surface (an embedded interface such as the Internal
+        # circle/sphere of the *InternalBoundary meshes) is snapped-to by
+        # adapt() so refinement follows its true geometry, but is NOT
+        # slip-eligible: interface motion is physics-owned (deform /
+        # free-surface machinery), so the movers keep its nodes pinned.
+        self.interior = bool(interior)
         self.centre = None if centre is None else np.asarray(centre, dtype=float).ravel()
         self.radius = None if radius is None else _as_float(radius)
         self.point = None if point is None else np.asarray(point, dtype=float).ravel()
@@ -210,15 +217,20 @@ def __repr__(self):
 
 
 # -- constructor-side registration helpers ---------------------------------
-def register_radial_surfaces(mesh, centre, label_radius):
+def register_radial_surfaces(mesh, centre, label_radius, interior=()):
     """Register ``radial`` surfaces: ``label_radius = {label_name: radius}``.
 
     Called by analytic radial-boundary constructors (Annulus, SphericalShell,
-    CubedSphere). ``centre`` is the common centre (e.g. the origin)."""
+    CubedSphere). ``centre`` is the common centre (e.g. the origin). Labels
+    named in ``interior`` are embedded interfaces: adapt() snaps refinement
+    onto them, but they stay pinned under the movers (interface motion is
+    physics-owned)."""
     centre = np.asarray(centre, dtype=float).ravel()
+    interior = set(interior)
     for lab, r in label_radius.items():
         mesh.register_tangent_slip_provider(
-            lab, BoundingSurface(mesh, lab, "radial", centre=centre, radius=r))
+            lab, BoundingSurface(mesh, lab, "radial", centre=centre, radius=r,
+                                 interior=lab in interior))
 
 
 def register_plane_surfaces(mesh, label_plane):

src/underworld3/meshing/spherical.py

@@ -780,6 +780,17 @@ def spherical_mesh_refinement_callback(dm):
         sympy.Matrix([-y, x, 0]),
     ]
 
+    # Radial bounding surfaces: Upper/Lower slip+snap as usual; the embedded
+    # Internal sphere is INTERIOR — adapt() snaps refinement onto its true
+    # radius, but the movers keep its nodes pinned (interface motion is
+    # physics-owned, 2026-07 round-3b ruling).
+    from underworld3.meshing.bounding_surface import register_radial_surfaces
+    register_radial_surfaces(
+        new_mesh, centre=(0.0, 0.0, 0.0),
+        label_radius={"Upper": radiusOuter, "Lower": radiusInner,
+                      "Internal": radiusInternal},
+        interior=("Internal",))
+
     return new_mesh
 
 

src/underworld3/utilities/custom_mg.py

@@ -888,7 +888,15 @@ def auto_inject_custom_mg(solver, field_id=None):
             solver.snes.setUp()
             op_n = int(solver.snes.getJacobian()[0].getSize()[0])
             pr, pc = (int(v) for v in Ps[-1].getSize())
-            if op_n > 0 and (pr != op_n or pc >= pr):   # rows!=op or no coarsening
+            # Only a genuine size mismatch disqualifies the transfer. A
+            # SQUARE finest transfer (pc == pr) is legitimate and common on
+            # boundary-focused refinement: a generation whose new vertices
+            # all land on a Dirichlet boundary adds only CONSTRAINED dofs,
+            # so the free-dof counts of the last two levels coincide — the
+            # level is redundant but correct, and rejecting it silently
+            # downgraded every curved-domain boundary-layer adapt to the
+            # default preconditioner (round-3b annulus finding, 2026-07).
+            if op_n > 0 and pr != op_n:
                 import warnings
                 warnings.warn(
                     "custom_mg: mesh-owned adapt-mesh FMG transfer is incompatible "

src/underworld3/utilities/nvb.py

@@ -70,8 +70,15 @@ def write_tagged_state_label(dm):
     perms = {p: i for i, p in enumerate(itertools.permutations(range(4)))}
     cS, cE = dm.getHeightStratum(0)
     vS, vE = dm.getDepthStratum(0)
-    if not dm.hasLabel(TAGGED_STATE_LABEL):
-        dm.createLabel(TAGGED_STATE_LABEL)
+    # RE-seeding must start from a FRESH label: DMLabelSetValue does not
+    # remove a point from its previous stratum, so overwriting an existing
+    # seed leaves every cell in TWO strata and readers get the OLD value
+    # back — the driver then runs the new geometry with the stale seed and
+    # the conforming drain deadlocks (found by the moved-base composition,
+    # 2026-07 round 3a).
+    if dm.hasLabel(TAGGED_STATE_LABEL):
+        dm.removeLabel(TAGGED_STATE_LABEL)
+    dm.createLabel(TAGGED_STATE_LABEL)
     label = dm.getLabel(TAGGED_STATE_LABEL)
     comm = dm.comm.tompi4py()