Refined-submesh (coarse/fine companion) investigation

docs/developer/design/submesh-solver-architecture.md

@@ -21,6 +21,93 @@ Underworld3 needs to support solving different equations on different subsets of
 
 5. **Normalised `Gamma_N`** (merged) — `mesh.Gamma_N` now returns a unit normal. Penalty and Nitsche BCs are mesh-independent.
 
+## Two Submesh Flavours: Subdomain and Resolution Level
+
+A *submesh* in UW3 is any mesh pulled out of a parent mesh that retains a
+lineage link (`parent`, registration in `parent._registered_submeshes`)
+and supports explicit field transfer back and forth. There are two
+flavours, and they share one usage pattern:
+
+> **get a submesh → build a solver on it → map fields back and forth**
+
+| | Subdomain | Resolution level |
+|---|---|---|
+| Constructor | `mesh.extract_region("Inner")` | `coarsened_companion(fine, levels=1)` |
+| PETSc mechanism | `DMPlexFilter` (cell subset) | `dm.refine()` nested hierarchy |
+| Parent ↔ child map | `getSubpointIS()` (point-level) | nested `createInjection` / `createInterpolation` (DOF-level) |
+| Transfer fidelity | exact (shared nodes) | exact (nested FE) |
+| Example | `docs/examples/submesh_investigation/test_region_ds_submesh.py` | `docs/examples/submesh_investigation/example_refined_companion.py` |
+
+Both examples solve the **same** annulus + radial-buoyancy Stokes problem
+so the flavours are directly comparable: one solves on a *subdomain* of
+the annulus, the other on a *coarser resolution* of the whole annulus,
+and both map the solution back to the parent.
+
+### Design contract for the resolution-level flavour: refine-DM mode only
+
+A resolution level is extractable **only when a genuine nested
+refinement hierarchy exists** (the mesh was built with
+`refinement >= 1`). The hierarchy is the source of truth; any level can
+be pulled out as a standalone solver-ready `uw.Mesh`, exactly as a
+subdomain is pulled out with `extract_region`.
+
+**If there is no refinement relationship the operation is unavailable.**
+There is deliberately:
+
+- **no geometric / `DMPlexComputeInterpolatorGeneral` fallback**, and
+- **no KDTree coordinate-matching fallback**.
+
+`coarsened_companion` raises a clear error on a non-refined mesh rather
+than silently degrading to an approximate transfer.
+
+Transfer between levels uses PETSc's *nested* interpolator/injector
+(`plex.c:10328`, `DMPlexComputeInterpolatorNested` /
+`DMPlexComputeInjectorFEM`), which is:
+
+- **exact** — the prolongation is the FE embedding; injection is a pure
+  scatter of coincident DOFs;
+- **parallel-local** — the injector builds per-rank `COMM_SELF` scatters
+  (`plexfem.c:3739-3741`); `DMRefine` preserves the coarse partition, so
+  no MPI communication is needed for restrict/prolongate;
+- **not dependent on point location** — no grid hashing, no geometry
+  search.
+
+```{note}
+The nested path triggers only when the fine DM's `getCoarseDM()` *is*
+the coarse DM and the regular-refinement flag is set. Independently
+cloning two hierarchy levels breaks this and petsc4py exposes no setter
+to restore it. The working construction is to build the transfer pair by
+**refining a single-field clone of the coarse level**
+(`dm_f = dm_c.refine()`): `refine()` itself establishes the linkage and
+the flag. See `refined_pair_prototype.py` (`_linked_pair`).
+```
+
+Empirically (box and annulus, P2 velocity), prolongating a coarse Stokes
+solution to the fine mesh and sampling it back recovers the coarse
+solution to **O(1e-15)** — machine precision — with the geometric escape
+hatch deliberately removed, proving the transfer is genuinely the nested
+operator.
+
+```python
+fine   = uw.meshing.Annulus(radiusOuter=1.5, radiusInner=0.5,
+                            cellSize=1/16, refinement=2)
+coarse = coarsened_companion(fine, levels=1)   # parent = fine
+
+v_c = uw.discretisation.MeshVariable("V", coarse, coarse.dim, degree=2)
+stokes = Stokes(coarse, velocityField=v_c, ...)   # solve cheaply, coarse
+stokes.solve()
+
+v_f = uw.discretisation.MeshVariable("Vf", fine, fine.dim, degree=2)
+prolongate(coarse, v_c, v_f)                      # exact, fills all fine DOFs
+```
+
+Status: prototype + both parallel examples + gating tests
+(`test_refined_pair_solver.py`) + contract test
+(`test_refined_pair_contract.py`) all passing. The investigation under
+`docs/examples/submesh_investigation/` is the sign-off artifact;
+promotion of `coarsened_companion` into the `Mesh` API proper is a
+follow-up.
+
 ## Design Principles
 
 ### 1. Separate meshes, separate variables, explicit copies

docs/examples/submesh_investigation/example_refined_companion.py

@@ -0,0 +1,148 @@
+"""
+Refined-companion approach: pull a COARSE level out of a refined mesh's
+nested hierarchy, solve Stokes on it, map the solution back to the fine
+mesh exactly.
+
+This is the sibling of ``test_region_ds_submesh.py``. Both follow the
+same pattern -- *get a submesh, build a solver, map back and forth* --
+but the submesh here is a different *resolution* of the whole domain
+rather than a *subdomain*:
+
+    test_region_ds_submesh.py : extract_region("Inner")   -> subdomain
+    example_refined_companion.py : coarsened_companion(...) -> coarse level
+
+The same annulus + radial-buoyancy problem is used so the two examples
+are directly comparable.
+
+Design contract (refine-DM mode only): the coarse companion is available
+ONLY because a genuine nested refinement hierarchy exists. Transfer
+between levels uses PETSc's *nested* interpolator/injector -- exact,
+parallel-local, no geometric point location, no KDTree. On a mesh with
+no refinement relationship the companion is simply not offered (raises).
+
+Usage:
+    pixi run -e amr-dev python -u \
+        docs/examples/submesh_investigation/example_refined_companion.py
+"""
+
+import numpy as np
+import sympy
+
+import underworld3 as uw
+from underworld3.systems import Stokes
+
+import refined_pair_prototype as rpp
+
+# --- Parameters ---
+
+r_outer = 1.5
+r_inner = 0.5
+cellsize = 1 / 16
+refine_levels = 2          # fine mesh = 2 uniform refinements of the base
+companion_levels = 1       # solve one level coarser than the finest
+n = 2
+k = 1
+stokes_tol = 1.0e-6
+vel_penalty = 1.0e6
+
+# --- Fine mesh (the full-resolution domain) ---
+
+uw.pprint(0, "Creating fine annulus mesh (with refinement hierarchy)...")
+fine = uw.meshing.Annulus(
+    radiusOuter=r_outer,
+    radiusInner=r_inner,
+    cellSize=cellsize,
+    refinement=refine_levels,
+)
+fS, fE = fine.dm.getHeightStratum(0)
+uw.pprint(
+    0,
+    f"Fine mesh: {fE - fS} cells, "
+    f"hierarchy depth {len(fine.dm_hierarchy)}",
+)
+
+# --- Pull out the coarse companion (the "submesh") ---
+
+uw.pprint(0, "Pulling coarse level out of the nested hierarchy...")
+coarse = rpp.coarsened_companion(fine, levels=companion_levels)
+cS, cE = coarse.dm.getHeightStratum(0)
+uw.pprint(0, f"Coarse companion: {cE - cS} cells")
+uw.pprint(0, f"  parent is fine mesh: {coarse.parent is fine}")
+uw.pprint(0, f"  registered with parent: {coarse in fine._registered_submeshes}")
+uw.pprint(0, f"  boundaries: {[b.name for b in coarse.boundaries]}")
+
+# --- Variables on the coarse companion ---
+
+v = uw.discretisation.MeshVariable("V", coarse, coarse.dim, degree=2)
+p = uw.discretisation.MeshVariable("P", coarse, 1, degree=1, continuous=True)
+
+# --- Coordinate system (same radial-buoyancy problem as the rock/air ex) ---
+
+unit_rvec = coarse.CoordinateSystem.unit_e_0
+r, th = coarse.CoordinateSystem.xR
+Gamma = coarse.Gamma
+
+# --- Stokes solver on the coarse companion ---
+
+stokes = Stokes(coarse, velocityField=v, pressureField=p)
+stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+stokes.saddle_preconditioner = 1.0
+
+rho = ((r / r_outer) ** k) * sympy.cos(n * th)
+stokes.bodyforce = rho * (-1.0 * unit_rvec)
+
+# Free-slip on both annulus boundaries
+stokes.add_natural_bc(vel_penalty * Gamma.dot(v.sym) * Gamma, "Upper")
+stokes.add_natural_bc(vel_penalty * Gamma.dot(v.sym) * Gamma, "Lower")
+
+stokes.tolerance = stokes_tol
+
+uw.pprint(0, "Solving Stokes on the coarse companion...")
+stokes.solve(verbose=False)
+
+v_mag = np.sqrt(v.data[:, 0] ** 2 + v.data[:, 1] ** 2)
+uw.pprint(
+    0,
+    f"Coarse solve: |v| max={v_mag.max():.6e} mean={v_mag.mean():.6e}",
+)
+
+# --- Map the coarse solution back to the fine mesh ---
+
+uw.pprint(0, "Mapping coarse velocity back to the fine mesh (nested FE)...")
+v_fine = uw.discretisation.MeshVariable("Vf", fine, fine.dim, degree=2)
+rpp.prolongate(coarse, v, v_fine)
+
+vf_mag = np.sqrt(v_fine.data[:, 0] ** 2 + v_fine.data[:, 1] ** 2)
+uw.pprint(
+    0,
+    f"Prolongated to fine: |v| max={vf_mag.max():.6e} "
+    f"mean={vf_mag.mean():.6e}",
+)
+
+# Round-trip: sampling the prolongated fine field back to the coarse
+# companion must recover the coarse solution exactly (it lives in the
+# fine FE space by construction).
+v_back = uw.discretisation.MeshVariable("Vb", coarse, coarse.dim, degree=2)
+rpp.sample(coarse, v_fine, v_back)
+rt_err = np.linalg.norm(v_back.data - v.data) / np.linalg.norm(v.data)
+
+# --- Report ---
+
+uw.pprint(0, "=" * 60)
+uw.pprint(0, "Refined-companion approach (coarse level of refined mesh)")
+uw.pprint(0, f"  Coarse cells:        {cE - cS}")
+uw.pprint(0, f"  Fine cells:          {fE - fS}")
+uw.pprint(0, f"  Max |v| (coarse):    {v_mag.max():.10e}")
+uw.pprint(0, f"  Max |v| (-> fine):   {vf_mag.max():.10e}")
+uw.pprint(
+    0,
+    f"  Prolongation max|v| rel diff: "
+    f"{abs(vf_mag.max() - v_mag.max()) / v_mag.max():.3e}",
+)
+uw.pprint(0, f"  Round-trip rel error (sample o prolongate): {rt_err:.3e}")
+uw.pprint(0, "  -> transfer is exact (nested FE): differences are O(1e-15)")
+uw.pprint(0, "=" * 60)
+
+assert rt_err < 1.0e-8, f"nested round-trip not exact: {rt_err:.3e}"
+uw.pprint(0, "example_refined_companion: OK")

docs/examples/submesh_investigation/probe_hierarchy_labels.py

@@ -0,0 +1,80 @@
+"""Probe: do boundary labels survive on the coarse DMs in dm_hierarchy?
+
+Gating question for the refined-submesh-pair investigation. If the coarse
+levels of the refinement hierarchy don't carry the named boundary labels
+(with non-empty strata), we cannot build a solver-ready coarse companion
+and the whole approach needs rethinking.
+
+Run: pixi run -e amr-dev python -u docs/examples/submesh_investigation/probe_hierarchy_labels.py
+"""
+
+import underworld3 as uw
+
+
+def dump_hierarchy_labels(mesh, label):
+    print(f"\n=== {label} ===")
+    hier = mesh.dm_hierarchy
+    print(f"dm_hierarchy length: {len(hier)}  (index 0 = coarsest, -1 = finest)")
+    bnames = [b.name for b in mesh.boundaries] if mesh.boundaries is not None else []
+    print(f"boundaries enum: {bnames}")
+
+    for lvl, dm in enumerate(hier):
+        cStart, cEnd = dm.getHeightStratum(0)
+        ncells = cEnd - cStart
+        print(f"\n  level {lvl}: {ncells} cells")
+        for b in mesh.boundaries:
+            if b.name in ("Null_Boundary", "All_Boundaries"):
+                continue
+            lab = dm.getLabel(b.name)
+            if not lab:
+                print(f"    {b.name:<20} value={b.value:<5} LABEL ABSENT")
+                continue
+            sis = lab.getStratumIS(b.value)
+            size = sis.getSize() if sis else 0
+            flag = "" if size > 0 else "  <-- EMPTY"
+            print(f"    {b.name:<20} value={b.value:<5} stratum size={size}{flag}")
+
+
+def main():
+    # 2D box (simplex, gmsh-generated boundary labels)
+    box = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0),
+        maxCoords=(1.0, 1.0),
+        cellSize=0.25,
+        refinement=2,
+    )
+    dump_hierarchy_labels(box, "UnstructuredSimplexBox 2D, refinement=2")
+
+    # Annulus with an internal boundary (different boundary-construction path).
+    # AnnulusInternalBoundary has no `refinement` kwarg, so build it and
+    # re-wrap its DM through the Mesh constructor with refinement=2 — this is
+    # exactly the DM-passthrough path coarsened_companion will rely on.
+    annulus_base = uw.meshing.AnnulusInternalBoundary(
+        radiusOuter=1.0,
+        radiusInternal=0.7,
+        radiusInner=0.5,
+        cellSize=0.2,
+    )
+    annulus = uw.discretisation.Mesh(
+        annulus_base.dm,
+        boundaries=annulus_base.boundaries,
+        coordinate_system_type=annulus_base.CoordinateSystemType,
+        refinement=2,
+        qdegree=annulus_base.qdegree,
+    )
+    dump_hierarchy_labels(
+        annulus, "AnnulusInternalBoundary (re-wrapped DM), refinement=2"
+    )
+
+    # Spherical shell (3D, yet another path)
+    shell = uw.meshing.SphericalShell(
+        radiusOuter=1.0,
+        radiusInner=0.5,
+        cellSize=0.4,
+        refinement=1,
+    )
+    dump_hierarchy_labels(shell, "SphericalShell 3D, refinement=1")
+
+
+if __name__ == "__main__":
+    main()